4th Order RungeKutta Methods

1
4th- Order Runge-Kutta Methods
The most commonly used set is

Four Evaluations of the right hand side, f(x,y),
are needed
2
Subroutine RK4 (x, y, h, yout) Real, intent(in)
x, y, h Real, intent(out) yout
Real hh, xh Real
dydx, yt Real
k1, k2, k3, k4 hh h / 2. xh
x hh Call derivatives ( x, y, dydx) k1 h
dydx yt y k1 / 2. Call derivatives ( xh, yt,
dydx) k2 h dydx yt y k2 / 2. Call
derivatives ( xh, yt, dydx) k3 h dydx yt y
k3 Call derivatives ( xh, yt, dydx) k4 h
dydx yout y (k1 2.k2 2.k3 k4) /
6. Return End subroutine RK4
Subroutine for one 4th order Runge Kutta step
(one 1st order ODE)

3
Subroutine RK4 (x, y, h, yout) Real, intent(in)
x, y, h Real, intent(out) yout
Real hh, xh Real
dydx, yt Real
k1, k2, k3, k4 hh h / 2. xh
x hh Call derivatives ( x, y, dydx) k1 h
dydx yt y k1 / 2. Call derivatives ( xh, yt,
dydx) k2 h dydx yt y k2 / 2. Call
derivatives ( xh, yt, dydx) k3 h dydx yt y
k3 Call derivatives ( xh, yt, dydx) k4 h
dydx yout y (k1 2.k2 2.k3 k4) /
6. Return End subroutine RK4
Subroutine for one 4th order Runge Kutta step
(one 1st order ODE)

4
Program to solve one 1st order differential
equation
Program Example_for_Runge_Kutta Declarations
x x_initial y y_initiaL h (x_end -
x_start)/ nstep Do i 1, nstep call rk4 (x,
y, h, yout) x x h y yout write(,)
x, y End do End Program Example_for_Runge_Kutta

5
Subroutine rk4 (x, y, n, h, yout) Integer,
intent(in) n Real, intent(in)
y(n), h, x Real, intent(out) yout(n)
Real hh, xh Real
dydx(n), yt(n) Real
k1(n), k2(n), k3(n),
k4(n) Integer i hh h /
2. xh x hh Call derivatives( x, y, n,
dydx) Do i 1, n k1(i) h dydx(i) yt(i)
y(i) k1(i) / 2. End do Call derivatives( xh,
yt, n, dydx) Do i 1, n k2(i) h dydx(i)
yt(i) y(i) k2(i) / 2. End do
Subroutine for one 4th order Runge Kutta step
(set of n 1st ODEs)
In this case y, yout, dydx, are all vectors of
length n.

6
Program to solve a set of n 1st order
differential equation

7
Adaptive Stepsize Control
Consider the following set of ODEs
with the Boundary Conditions
The analytical solution is

8
Adaptive Stepsize Control
The characteristic scale of the solution
decreases rapidly
Question What should we chose for our stepsize
?t?

Answer ?t should adjust to the changing
characteristic scales
---gtgt Adaptive Stepsize Control
9
Adaptive Stepsize Control
Goal of adaptive stepsize control
Achieve some predetermined accuracy in the
solution with minimum computational effort!
Numerical Recipes Many small steps should
tiptoe through treacherous terrain, while a few
great strides should speed through smooth
uninteresting countryside.
Adaptive stepsize control --gt requires
constantly monitor the solution and return
information about its performance
Most important is an estimation of the truncation
error

Truncation error is used to control the stepsize!
10
Step Doubling
A simple implementation of an adaptive stepsize
control is step doubling
Take each step twice
1. As a full step
2. As two half steps
Difference between the two answers is used to
estimate the local truncation error.

11
Adaptive Stepsize Control
We know that the 4th order RK has an error O(h5)
Lets denote the exact solution for an advance
from x to x2h by y(x2h)

The function ? remains to order h5 constant over
the step

The difference of the two numerical estimates
ybig and ysmall is a convenient indicator of the
truncation error
12
Adaptive Stepsize Control


13
Adaptive Stepsize Control
How can we use this information to control the
stepsize?
Suppose that at the current stepsize hc we found
that the error to be
This is our estimate of the truncation error.

We want this error to be less or equal to our
specified ideal error, call it ?i

We can estimate the required stepsize to be
14
Some Practical Considerations
Since hest is only an estimate we want to set the
actual new stepsize a little smaller than the
estimated value
A typical value for S1 is S1 0.9
A second safety factor, S2 gt 1, is often used to
ensure that the program does not raise or lower
the stepsize too enthusiastically
Prevents a too large increase
Prevents a too large decrease
Reasonable Variations
A typical value for S2 is S2 4.0
15
Some Practical Considerations
Our notation of ? is a little misleading.
For a set of n 1st order ODEs it is actually a
vector of the desired accuracy for each equation
in the set of ODEs.
In general all equations need to be within their
respective accuracies.
We need to scale our stepsize according to the
worst-offender equation
16
Some Practical Considerations
How to chose ?I ?
Often we say something like The solution has to
be good to within one part in 106.
In this case you might want to use fractional
errors
where eps is for example a number like 10-6
If your solution, however passes through zero
than this scheme can become problematic
In this case a useful trick is
17
Global Error Constraints
How can we chose a cautious stepsize adjustor
that will control the global error?
The local error due to one step is related to the
stepsize by
If we take m steps of size hest across the entire
interval x1, x2
We can estimate the required stepsize to achieve
a global error control
Slightly more stringent than local error control
18
Some Practical Considerations
(evaluates one single step)
Subroutine Adaptive_Runge_Kutta Inputs x,
y(n), h_initial, ?i, n Output xout, yout(n),
h_new Set initial variables (e.g., max number
of attempts, S1, S2, ) Loop over maximum
number of attempts to satisfy error bounds -
Take the two small steps (call RK4) - Take the
single big step (call RK4) - Compute the
estimated truncation error - Estimate the new
hnew - If error is acceptable, return computed
values Issue error message if error bound is
never satisfied