Prev | Next |

`# include <cppad/runge_45.hpp>`

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

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

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{array}{rcl}

X(ti) & = & xi \\

X'(t) & = & F[t , X(t)]

\end{array}

\]

If your set of ordinary differential equations
are stiff, an implicit method may be better
(perhaps Rosen34
.)
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).
The file

`cppad/runge_45.hpp`

is included by `cppad/cppad.hpp`

but it can also be included separately with out the rest of
the `CppAD`

routines.
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
.
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*)

The argument

*t*

to
*F*.Ode

has prototype

const *Scalar* &*t*

(see description of Scalar
below).
The argument

*x*

to
*F*.Ode

has prototype

const *Vector* &*x*

and has size
*n*

(see description of Vector
below).
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`

.
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*

.
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)
```

.
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*

.
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*

.
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)
```

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
.
The type

*Scalar*

must satisfy the conditions
for a NumericType
type.
The routine CheckNumericType
will generate an error message
if this is not the case.
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*)

.
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.
For each set of types Scalar , Vector , and Fun , the first call to

`Runge45`

must not be parallel
execution mode.
The file runge45_1.cpp contains a simple example and test of

`Runge45`

.
It returns true if it succeeds and false otherwise.

The file
runge45_2.cpp
contains an example using `Runge45`

in the context of algorithmic differentiation.
It also returns true if it succeeds and false otherwise.
The source code for this routine is in the file

`cppad/runge_45.hpp`

.
Input File: cppad/runge_45.hpp