\(\newcommand{\W}[1]{ \; #1 \; }\) \(\newcommand{\R}[1]{ {\rm #1} }\) \(\newcommand{\B}[1]{ {\bf #1} }\) \(\newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} }\) \(\newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} }\) \(\newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} }\) \(\newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }\)

Runge45

View page source

An Embedded 4th and 5th Order Runge-Kutta ODE Solver

Syntax

# include <cppad/utility/runge_45.hpp>

xf = Runge45 ( F , M , ti , tf , xi )

xf = Runge45 ( F , M , ti , tf , xi , e )

Purpose

This is an implementation of the Cash-Karp embedded 4th and 5th order Runge-Kutta ODE solver described in Section 16.2 of Numerical Recipes . We use \(n\) for the size of the vector xi . Let \(\B{R}\) denote the real numbers and let \(F : \B{R} \times \B{R}^n \rightarrow \B{R}^n\) be a smooth function. The return value xf contains a 5th order approximation for the value \(X(tf)\) where \(X : [ti , tf] \rightarrow \B{R}^n\) is defined by the following initial value problem:

\begin{eqnarray} X(ti) & = & xi \\ X'(t) & = & F[t , X(t)] \end{eqnarray}

If your set of ordinary differential equations are stiff, an implicit method may be better (perhaps Rosen34 .)

Operation Sequence

The operation sequence for Runge does not depend on any of its Scalar input values provided that the operation sequence for

F . Ode ( t , x , f )

does not on any of its Scalar inputs (see below).

Include

The file cppad/utility/runge_45.hpp is included by cppad/cppad.hpp but it can also be included separately with out the rest of the CppAD routines.

xf

The return value xf has the prototype

Vector xf

and the size of xf is equal to n (see description of Vector below).

\[X(tf) = xf + O( h^6 )\]

where \(h = (tf - ti) / M\) is the step size. If xf contains not a number nan , see the discussion for f .

Fun

The class Fun and the object F satisfy the prototype

Fun & F

The object F (and the class Fun ) must have a member function named Ode that supports the syntax

F . Ode ( t , x , f )

t

The argument t to F . Ode has prototype

const Scalar & t

(see description of Scalar below).

x

The argument x to F . Ode has prototype

const Vector & x

and has size n (see description of Vector below).

f

The argument f to F . Ode has prototype

Vector & f

On input and output, f is a vector of size n and the input values of the elements of f do not matter. On output, f is set equal to \(F(t, x)\) in the differential equation. If any of the elements of f have the value not a number nan the routine Runge45 returns with all the elements of xf and e equal to nan .

Warning

The argument f to F . Ode must have a call by reference in its prototype; i.e., do not forget the & in the prototype for f .

M

The argument M has prototype

size_t M

It specifies the number of steps to use when solving the differential equation. This must be greater than or equal one. The step size is given by \(h = (tf - ti) / M\), thus the larger M , the more accurate the return value xf is as an approximation for \(X(tf)\).

ti

The argument ti has prototype

const Scalar & ti

(see description of Scalar below). It specifies the initial time for t in the differential equation; i.e., the time corresponding to the value xi .

tf

The argument tf has prototype

const Scalar & tf

It specifies the final time for t in the differential equation; i.e., the time corresponding to the value xf .

xi

The argument xi has the prototype

const Vector & xi

and the size of xi is equal to n . It specifies the value of \(X(ti)\)

e

The argument e is optional and has the prototype

Vector & e

If e is present, the size of e must be equal to n . The input value of the elements of e does not matter. On output it contains an element by element estimated bound for the absolute value of the error in xf

\[e = O( h^5 )\]

where \(h = (tf - ti) / M\) is the step size. If on output, e contains not a number nan , see the discussion for f .

Scalar

The type Scalar must satisfy the conditions for a NumericType . The routine CheckNumericType will generate an error message if this is not the case.

fabs

In addition, the following function must be defined for Scalar objects a and b

a = fabs ( b )

Note that this operation is only used for computing e ; hence the operation sequence for xf can still be independent of the arguments to Runge45 even if

fabs ( b ) = std::max ( - b , b )

.

Vector

The type Vector must be a SimpleVector class with elements of type Scalar . The routine CheckSimpleVector will generate an error message if this is not the case.

Parallel Mode

For each set of types Scalar , Vector , and Fun , the first call to Runge45 must not be parallel execution mode.

Example

The file runge45_1.cpp contains a simple example and test of Runge45 .

The file runge_45.cpp contains an example using Runge45 in the context of algorithmic differentiation. It also returns true if it succeeds and false otherwise.

Source Code

The source code for this routine is in the file cppad/runge_45.hpp .