Example Code

Here is a simple example code that implements the Fourth-Order Runge-Kutta:

SUBROUTINE rk4(y, dydx, x, h, yout, derivs) USE nrtype; USE nrutil, ONLY : assert_eq IMPLICIT NONE REAL(sp), DIMENSION(:), INTENT(in) :: y, dydx REAL(sp), INTENT(in) :: x, h REAL(sp), DIMENSION(:), INTENT(out) :: yout INTERFACE SUBROUTINE derivs(x, y, dydx) USE nrtype IMPLICIT NONE REAL(sp), INTENT(in) :: x REAL(sp), DIMENSION(:), INTENT(in) :: y REAL(sp), DIMENSION(:), INTENT(out) :: dydx END SUBROUTINE derivs END INTERFACE INTEGER(i4b) :: ndum REAL(sp) :: h6, hh, xh REAL(sp), DIMENSION(SIZE(y)) :: dym, dyt, yt ndum = assert_eq(SIZE(y), SIZE(dydx), SIZE(yout), 'rk4') hh = h*0.5_sp h6 = h/6.0_sp xh = x + hh yt = y + hh * dydx CALL derivs(xh, yt, dyt) yt = y + hh * dyt CALL derivs(xh, yt, dym) yt = y + h * dym dym = dyt + dym CALL derivs(x + h, yt, dyt) yout = y + h6 * (dydx + dyt + 2.0_sp * dym) END SUBROUTINE rk4

The subroutine applies Runge-Kutta solver to a system
of equations:

The computation is carried on in parallel, assuming a parallelising Fortran compiler.

Note that here the lower index numbers the equations not
time steps. Furthermore the notation itself is
**y**(*x*) rather than
**x**(*t*).

The arguments that must be passed to the subroutine are

**y**- which is the vector of initial values of
**y**_{n}(in our notation this would be**x**_{n}) at ``time''*x*_{n}(in our notation this would be*t*_{n}); **dydx**- which is the vector of what in our notation
would be
**f**(**x**_{n},*t*_{n}), and in this program's notation, these are the values that the right hand side assumes at ``time''*x*_{n}; **x**- which is the value of time, in our notation
*t*_{n}; **h**- which is the length of the time step, in our notation ;
**yout**- which is a place for the new, updates values
of
**y**(*x*+*h*); in our notation this would be ; **derivs**- which is the name of the function that
corresponds to our function
*f*(**x**,*t*), and whose job is to evaluate dydx at various stages of the computation.

The first operation simply checks if arrays
`y`

, `dydx`

, and `yout`

are of matching sizes,
and, if not, an appropriate error message is written on
standard output. Otherwise the size is written on a dummy
variable `ndum`

that is not used any more:

ndum = assert_eq(SIZE(y), SIZE(dydx), SIZE(yout), 'rk4')The next three statements:

hh = h*0.5_sp h6 = h/6.0_sp xh = x + hhevaluate (

`hh`

),
(`h6`

), and
,
which goes into `xh`

.
And now we evaluate

and store it on

`yt`

.
This will go into yt = y + hh * dydx

The next step is to evaluate

where

`yt`

, and `derivs`

reverses the order of
arguments, compared with our notation, i.e., time (`xh`

)
goes first, then position (`yt`

), and the last argument
is used for the answer, once the subroutine returns:CALL derivs(xh, yt, dyt)What this function actually returns, is not

`dyt`

by anything yet.
This multiplication takes place in the next line, when we evaluate

yt = y + hh * dytBut remember that

`hh`

is really
,
so what we are evaluating here, in fact is:
so that the next call:

CALL derivs(xh, yt, dym)evaluates:

which really is , and returns it in

`dym`

.
Now we evaluate *x*_{n} + *k*_{3} thusly:

yt = y + h * dymRemember that whereas

`hh`

stands for
,
the single `h`

is just ,
so `h * dym`

is indeed plain `x + h`

stands
for, in our notation,
.
With these two
we can now evaluate CALL derivs(x + h, yt, dyt)What's returned in

`dyt`

is
.
Observe
that just before calling `derivs`

we have save
the previous value of
`dyt + dym`

on `dym`

. That previous
value was
.
So
this means that `h6 * 2.0_sp * dym`

stands
for

In turn

`h6 * (dydx + dyt)`

, stands for
Consequently, you can now see clearly that:

yout = y + h6 * (dydx + dyt + 2.0_sp * dym)evaluates

which is the Runge-Kutta formula.