Approaches to Runge Kutta step size adjustments

[barker7]

I've written a Runge Kutta algorithm that uses a fixed step size. It works fine for many of my problems dealing with chemical reaction kinetics but I would like to make it faster and more robust by implementing automatic step size adjustments.

I have looked at the at the step-doubling algorithm in Numerical Recipes (1986 Ed.), and I have also looked at the Runge Kutta Fehlberg (RKF45) algorithm in Numerical Mathematics textbook by Cheney and Kincaid.

1) Does anyone have any comments on advantages of either method?

2) I have found many references to RKF45 on the web, but strangely NR doesn't mention it. RKF45 predates my 1986 copy of NR by many years. NR is usually very comprehensive, so I'm wondering why this omission.

Thanks for any info on choosing a step size algorithm,

Mitch

 

[AlephZero]

The basic purpose of NR is to present methods that work reliably enough to be widely useful (see the preface and introduction). It isn't, and couldn't possibly be, a comprehensive survey of numerical methods.

My personal view (based on quite a lot of experience) is this:

If you use a numerical method to solve an ODE of lower than the method, inveitably the numerical solution contains functions that are not solutions of the ODE. If the extra functions decay with time, that is not a problem, but often they grow with time, which is the cause of contitional stability.

Most "real world" ODEs are first and second order, and high order RK methods are almost always only conditionally stable. However because of the high accuracy, in practice you find that if the results are stable they are also accurate.

So, the real task of auto time step control is RK is not to control accuracy, but to maintain stability. The very simple algorithm "Start with a fixed large step size, and if the numerical solution is obviously nonsense (i.e. diverging to infinity) halve the step size and start again from the beginning" often works as well as anything else.

If you want a really efficient ODE solver, you need one selected for the equations you are trying to solve, not a general method like RK. In some cases that can give you a speed-up of several orders of magnitude, using step size control with an unconditionally stable method.

 

[barker7]

I agree about the purpose of NR. However when I saw numerous mentions of RKF45 on the web, and in another textbook, I assumed RKF45 was a "common" technique. And that is why I though it was strange that it wasn't even mentioned in NR.

My ODEs are all first order.

Thanks for your insight on stability. In fact, I did implement something as your said where i would integrate and then evaluate if the step size was too large or too small, and then repeat with an adjusted step size. This worked OK, but since I was running the integration for hundreds of initial conditions to build a surface plot, I was looking for some speed boost. Also "most" chemical reaction are fast at the beginning but slow way down towards the end (sorta exponential decay shape). Thus I was hoping to really extend the step size as the reaction progressed.

BTW, I'm a ChemE/programmer so that's why I gravitated towards the general purpose RK algo.

Mitch

 

[AlephZero]

Also "most" chemical reaction are fast at the beginning but slow way down towards the end (sorta exponential decay shape). Thus I was hoping to really extend the step size as the reaction progressed.

If you want to do that, you need an unconditionally stable method with step size control.

For a first order equation and unconditional stability, you are limited to first or second-order methods. But first order methods are too inaccurate, unless you have a VERY nonlinear problem and the stability requirement wins over accuracy.

I would try something like the modified Euler method (not the simple Euler forward difference method which is only first order).

1. Guess the solution at the end of the timestep, by straight line extrapolation from the last two values.
2. Find the derivative at the end of the step from the ODE.
3. Find a better solution from x(h) = x(0) + h(x'(0)+x'(h)/2. You already know x'(0) from the previous step, so you don't need to recalculate it.
4. Repeat steps 2 and 3 till it converges. If it doesn't converge in say 5 iterations, halve the timestep and try again.

When you have a solution for 2 equal sized steps, estimate the error from the "curvature" of the solution, i.e. x(0) - 2 x(h) + x(2h). If this is outside some tolerance band, either double the next timestep, or reject the current step and halve it.

That should work well for linear or mildly nonlinear problems. Reducing the time step will hardly ever happen, unless there is a sharp change in an external forcing term in the equations. In that case you might want to reduce h to a small value when the change occurs, rather than letting the algorithm figure out what to do on its own.

If you have a very nonlinear problem, then try a backward difference method - exactly the same algorithm, except
x(h) = x(0) + hx'(h).

Some people like using a weighted average of these two methods in an attempt to get the best of both worlds in accuracy and stability, i.e.
x(h) = x(0) + h((1-k)x'(0) + kx'(h)) where k is a magic constant between 0.5 and 1. Personally I'm not convinced, despite have heard some very heated arguments about what value of k gives the best magic.