Prev Next

@(@\newcommand{\B}[1]{{\bf #1}} \newcommand{\R}[1]{{\rm #1}}@)@
Fourth Order Runge Kutta

Syntax
Purpose
f
ti
yi
dt
yf
Motivation
Source Code
Example

Syntax
 yf = runge_kutta_4(f, ti, yi, dt)

Purpose
See Motivation below as well as the purpose described here. We are given a function @(@ f : \B{R}^n \rightarrow \B{R}^n @)@, and a point @(@ yi \in \B{R}^n @)@ such that an unknown function @(@ y : \B{R} \rightarrow \B{R}^n @)@ satisfies the equations @[@ \begin{array}{rcl} y( ti ) & = & yi \\ y'(t) & = & f[t, y(t) ] \\ \end{array} @]@ We use the Fourth order Runge-Kutta formula (see equation 16.1.2 of Numerical Recipes in Fortran, 2nd ed.) wish to approximate the value of @[@ yf = y( ti + \Delta t ) @]@

f
If t is a scalar and y is a vector with size @(@ n @)@,
     k = f(t, y)
 returns a vector of size @(@ n @)@ that is the value of @(@ f(t, y) @)@ at the specified values.

ti
is a scalar that specifies the value of @(@ ti @)@ in the problem above.

yi
is a vector of size @(@ n @)@ that specifies the value of @(@ yi @)@ in the problem above.

dt
is a scalar that specifies the value of @(@ \Delta t @)@ in the problem above.

yf
is a vector of size @(@ n @)@ that is the approximation for @(@ y( t + \Delta t ) @)@.

Motivation
This routine makes very few assumptions about the objects used to do these calculations. Thus, smart objects can be used for all sorts of purposes; for example, computing derivatives of the solution of an ODE. The table below lists the assumptions on the objects passed into runge_kutta_4. In this table, s and t are scalars, d is a decimal number (i.e., a float) and u and v are vectors with size @(@ n @)@.
operation result
d * s scalar
s + t scalar
s * u vector with size @(@ n @)@
d * u vector with size @(@ n @)@
s * u vector with size @(@ n @)@
u + v vector with size @(@ n @)@

Source Code 
 
def runge_kutta_4(f, ti, yi, dt) :
     k1 = dt * f(ti         , yi)
     k2 = dt * f(ti + .5*dt , yi + .5*k1) 
     k3 = dt * f(ti + .5*dt , yi + .5*k2) 
     k4 = dt * f(ti + dt    , yi + k3)
     yf = yi + (1./6.) * ( k1 + 2.*k2 + 2.*k3 + k4 )
     return yf

Example
  1. The file runge_kutta_4_correct.py contains an example and test of using runge_kutta_4 to solve an ODE.
  2. The file runge_kutta_4_ad.py contains an example and test of differentiating the numerical solution of an ODE.
  3. The file runge_kutta_4_cpp.py contains an example and test of using pycppad adfun object to evaluate python functions with C++ speed of execution.
Input File: pycppad/runge_kutta_4.py