A system of differential equations can be solved using a numerical integrator with adaptive stepsize $ h$, e.g. Runge Kutta Fehlberg

Now suppose I want to integrate a set of differential equations $\{\dot{x}, \dot{y}, \dot{z} \}$ numerically (e.g. the 3D trajectory of some particle). Each differential equation uses the same adaptive stepsize $h$.

I want to evaluate $x,y$ at some specific target value of $z = z_t$

Now, due to the adaptive stepsize, it is unlikely that $z$ will ever be exactly $z_t$. That is to say, I might get for integration step $n$ and step $n+1$, $z_n=z_t + \epsilon$, $z_{n+1}=z_t - \delta$ for small values $\epsilon, \delta$

Obviously, by setting a more strict tolerance on the integrator, I can get close to $z_t$, but this would be at the expense of speed.

Given that I know a priori that I want to evaluate $x,y$ at $z=z_t$, is there a method which permits both an adaptive stepsize, and an exact $z=z_t$ result?


  • $\begingroup$ See for example stackoverflow.com/q/48892663/3088138 where zero-crossings of an oscillatory process are computed. The last attempt uses Newton's method for finding the best t. Note that odeint uses lsoda which uses adaptive step sizes. $\endgroup$ Feb 21 '18 at 16:59
  • $\begingroup$ Why do you want to control the position of the last step indirectly via the tolerances and not directly by cutting it to the desired length? $\endgroup$ Feb 21 '18 at 17:01
  • $\begingroup$ @LutzL Thanks for your reply. Can you explain what you mean by 'cutting it to the desired length'? $\endgroup$ Feb 21 '18 at 17:02
  • 1
    $\begingroup$ something like if t+1.01*dt > tf then dt = tf-t. The factor 1.01 is to avoid/preempt to cut to a ridiculously small step. $\endgroup$ Feb 21 '18 at 17:06

Once you choose a method, it has an order which approximates how the error depends on stepsize. The popular fourth order Runge-Kutta has an error that scales as $h^4$. You can get an approximation to the accuracy by doing your step once as a full step and once as two half steps. You would expect the error to be reduced by a factor $16$ so the difference between them is about $15$ times the error of the two step approach. You can use this to update your stepsize. There is a discussion in chapter 17.2 of Numerical Recipes and probably in many other numerical analysis texts.

  • $\begingroup$ The question was IMO not about implementing the adaptive step size, but how to implement what is sometimes called "events" on top of an adaptive scheme. $\endgroup$ Feb 22 '18 at 9:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.