Procedure for adaptive step size for Runge Kutta 4?

Question

I am writing a Runga Kutta 4 algorithm in MATLAB. I would like to add adaptive step sizing to this algorithm.

From what I've read it seems you calculate the value of the function for two step sizes on each iteration and then from the size of the error terms you deduce which one to use as $y_{i+1}$. But I can't see any advantage in efficiency if you are calculating two potential $y_{i+1}$'s on each iteration.

Obviously I am missing/misinterpreted something as I have read that adaptive step sizing can lead to huge performance increases. So how does adaptive step sizing work for RK methods and can someone tell me the steps I need to implement?

Answer 1

An old question, but one that I'm working on right now. As the person in the answer suggested, I've made an RKF45 ODE integrator and am trying to implement an adaptive step-size stepper.

The idea behind the stepper is that you calculate the fourth order step and the fifth order step. You compare the two of them to see what the error looks like with the fifth order taken as "correct." You will then check to see how the error (in my case, the maximum error from a vector of protein concentrations) compares to the desired error. If it is greater than the desired error, you reduce the step size and calculate the 4th and 5th order steps again, repeating the process until the error is smaller than the desired value.

Once the error is smaller than the desired value, or if it is smaller than the desired value, you can take a time step and advance all your dependent variables. Earlier literature that I read used the 4th order calculation to take the step, but later literature suggested just using the 5th order since you've already gone through the trouble of calculating it.

If your error is much smaller than the desired value, an adaptive step size algorithm should increase the step size to try and speed up the overall calculation. In Numerical Recipes, they suggest a stepper like this :

h0 = h1 * (DesiredError/CalculatedError)^(.2)

where h1 is the step size you tried to take and h0 is the theoretically appropriate step size.

Answer 2

I think for adaptive size, typically what is called a embedded method is used. Basically one set of coefficients gives the next step, and another set gives the error estimate. Look for Fehlberg method http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method

The advantage comes from the fact that the overall number of steps may be less, even though there are more evaluations for each step. Obviously this is not always the case and a fixed step RK is the best bet for many types of problems.