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â€¦

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?

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.