DifferentialEquations and a phase angle variable

Is there a way to use the solver from DifferentialEquations with a phase variable (where only the value within [0, 2pi] matters)?

The equation system generates a continuously increasing phi, so unless special precautions are taken, at some point phi becomes so large that the floating point resolution does not resolve values on the order of 1 anymore. And then the whole calculation is screwed.

Ideally a phase phi reaching 2pi would wrap around to 0 somehow. All parameters that are derived from phi only enter the equation system as argument of sin or cos.

Julia newbie here, so please bear with me if questions seem clueless…

You can define a callback so that when it goes over 2pi you just set it back to zero.

Or you can make a fun number type.

Sounds promising - do you have any pointers to documentation for me please?

https://diffeq.sciml.ai/stable/features/callback_functions/#Example-1:-Bouncing-Ball-1 is a good example

Ah great, thanks! -A

OK now I’m running into an “interesting” problem. Our code looks like this:

function condition(u, t, integrator)
     u[1] > 2*π

function affect!(integrator)
    integrator.u[1] = integrator.u[1] % (2*π)

and later

sol = solve(prob, maxiters=1e10, alg_hints=[:nonstiff], callback=cb, adaptive=false, dt=1.0, saveat=maxtime - maxtime/1000:maxtime)

Curiously, after setting the callback we 1) do not get dt=1 in the output anymore, but much larger time steps (only points where the callback is triggered???) and 2) the result covers the entire time from 0 to maxtime, while we’re only interested in the last part.
2) is not a huge problem, but 1) is… We’re interested in the region between the callbacks too… Any ideas on how to fix this?

You’re using a discrete callback? I assumed you’d want a continuous callback here.

Note that adaptive time stepping and saving are different ideas. dt and save points have no relation. I would be curious why adaptivity is off here.

Set save_positions = (false,false) to turn off callback saving.

1 Like

I’ve modeled sine wave sources (VCOs) by using equations of the form du[1] = f and V = A*sinpi(2*u[1]). I let the integrator turn frequency into phase (in units of cycles). Of course if your simulation runs long enough you will run into Patriot missile type bugs, but in my case the loss of resolution didn’t get anywhere near enough that I needed to do what Chris mentioned. Specifically, mine could reach as many as 500_000 cycles which meant the resolution within a cycle was reduced to about 34 bits (assuming Float64) , which I believed was good enough.

[EDIT] Even if you are simulating long enough for it to matter, you could relax on how often the callback has to fire to reset the phase - only fire it every 1000 cycles or 100_000 cycles, etc.