How to implement periodic boundary conditions with callbacks?

Specific Domains
,
1

I am trying to use 2 callbacks when solving a differential equation.
I am trying to solve for a pendulum, so I have a cyclic coordinate that goes from (-\pi, \pi). So when the coordinate goes above \pi I want it to shift by -2\pi, and when it goes below -\pi I want it to shift by +2\pi.

I tried to have two callbacks:

function condition_pi(u, t, integrator) 
    u[1] - π
end

function condition_min_pi(u, t, integrator) 
    u[1] + π
end

function affect_pi!(integrator)
    integrator.u[1] = integrator.u[1] - 2π
end

function affect_min_pi!(integrator)
    integrator.u[1] = integrator.u[1] + 2π
end

cb_pi = ContinuousCallback(condition_pi, affect_pi!)
cb_min_pi = ContinuousCallback(condition_min_pi, affect_min_pi!)
cb_set = CallbackSet(cb_pi, cb_min_pi)

This does not work. I think the reason is that when the solution hits x=\pi, it sends it back by -2\pi, which then makes it less than -\pi, and then it wants to send it forward by +2\pi, and there is an infinite loop.

So the question is, how can I implement these boundary conditions?

2

Seems a lot easier to allow your angles \theta to be arbitrary, write your equations in terms of functions like \sin(\theta) that are insensitive to shifts of 2\pi, and then transform the coordinates modulo 2\pi afterwards if you like (e.g. for plotting solutions).

1 Like
3

AFAIK, you need to set affect_neg argument to nothing.

4

I know. But I have other reasons that I want the coordinate to be between (-\pi, \pi). I have some functions that cannot be defined beyond this range.

5

What is affect_neg? I couldn’t find it in the documentation.

6

If you can change your calculation domain to go from [0, 2π) instead, then you can simply do

u -> mod(u[1], 2π)
7

This would be a pain for me since I have some other complicated functions that are defined for the range (-\pi, \pi).

8

Callbacks have both increasing and decreasing crossing affects. You don’t want negative crossing. By default, they are both set to the same affect.

9

Would this be numerically accurate enough for your usecase?

u -> mod(u[1] + π, 2π) - π
10

Could you elaborate? I don’t understand what you mean.

11

Nice idea!
I will try, though I need very good numerical accuracy.

12

Check the affect! and affect_neg! arguments here:
https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#ContinuousCallback

1 Like
13

Thanks!
That actually works:

Here is the code:

function affect_pi!(integrator)
    integrator.u[1] = integrator.u[1] - 2π
end

function affect_pi_neg!(integrator)
    integrator.u[1] = integrator.u[1]
end

function affect_min_pi!(integrator)
    integrator.u[1] = integrator.u[1]
end


function affect_min_pi_neg!(integrator)
    integrator.u[1] = integrator.u[1] + 2π
end

cb_pi = ContinuousCallback(condition_pi, affect_pi!, affect_pi_neg!)
cb_min_pi = ContinuousCallback(condition_min_pi, affect_min_pi!, affect_min_pi_neg!)
14

Why not just call ũ = mod2pi(u[1] + π) - π before calling those functions? Why does the ODE solver have to know that you are interpreting a coordinate cyclically?

2 Likes
15

That’s also a nice idea. I tested this, but it works less well than using the callback. Perhaps it’s an issue of the numerical accuracy of mod?
I am using ArbFloat to solve to try to achieve good accuracy.