Add dt to system arguments

[ilja.j.honkonen@…]

A popular algorithm for calculating the change in velocity of a charged particle in a magnetic field (1) requires the length of the time step while integrating (dx/dt := f(x,t,dt)) and cannot be used with odeint as far as I could tell so I suggest modifying the definition of system from sys( x , dxdt , t ) to sys( x , dxdt , t , dt ). Some enable_if magic should allow all existing code using the current API to work.

(1) Equations 25 and 26 in (h)ttp://www.dtic.mil/dtic/tr/fulltext/u2/a023511.pdf

[karsten]

This sound interesting.

I had a look at the paper you mention. If one implements a solver as described in Eq 25 and Eq 26 one could easily pass the step size to the system function. The requirements of such solver on the system function is then different than for the usual RK steppers and it might require a system in the form sys( x , dxdt , t , dt ).

But I am not sure if it make sense to add dt to the Runge-Kutta steppers. For example, for the classical RK4, one time step is calculated from 4 individual Euler steps each with a different step size. Do you think it is possible to use your suggestion with Runge Kutta methods?

[ilja.j.honkonen@…]

Do you think it is possible to use your suggestion with Runge Kutta methods?

Should be possible as eqs 25 and 26 only present the accelerator. Instead of (h)ttps://github.com/nasailja/pamhd/blob/master/tests/particle/odeint.cpp#L70 the acceleration step would be (h)ttps://github.com/nasailja/pamhd/blob/master/source/particle/solve.hpp#L79 Currently the problem is that the former using rkf78 allows long steps but still conserves particle energy worse than the latter with any step size while the latter requires very short steps in order to get a reasonable particle trajectory.

I thought about this a bit more and basically at the moment odeint assumes that there is a linear relation between time step and the resulting change but I guess there are many cases where this doesn't hold. One that came to mind is corner transport when solving the advection equation in multiple dimensions. First hit on google has a nice picture (6.2): (h)ttp://ocw.mit.edu/courses/earth-atmospheric-and-planetary-sciences/12-950-atmospheric-and-oceanic-modeling-spring-2004/lecture-notes/lec10.pdf In that case it might not be too difficult to come up with a dt-less solver but with more dimensions and complex geometries that will quickly become not worth the trouble.