Prev Next Index-> contents reference index search external Up-> CppAD utility Rosen34 CppAD-> Install Introduction AD ADFun preprocessor multi_thread utility ipopt_solve Example speed Appendix utility-> ErrorHandler NearEqual speed_test SpeedTest time_test NumericType CheckNumericType SimpleVector CheckSimpleVector nan pow_int Poly LuDetAndSolve RombergOne RombergMul Runge45 Rosen34 OdeErrControl OdeGear OdeGearControl CppAD_vector thread_alloc index_sort to_string set_union Rosen34-> rosen_34.cpp Headings-> Syntax Description Include xf Fun ---..t ---..x ---..f ---..f_t ---..f_x ---..Nan ---..Warning ---..Optimization M ti tf xi e Scalar Vector Parallel Mode Example Source Code

A 3rd and 4th Order Rosenbrock ODE Solver

Syntax
# include <cppad/utility/rosen_34.hpp>  xf = Rosen34(F, M, ti, tf, xi)  xf = Rosen34(F, M, ti, tf, xi, e) 
Description
This is an embedded 3rd and 4th order Rosenbrock ODE solver (see Section 16.6 of Numerical Recipes for a description of Rosenbrock ODE solvers). In particular, we use the formulas taken from page 100 of Shampine, L.F. (except that the fraction 98/108 has been correction to be 97/108).  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 not stiff an explicit method may be better (perhaps Runge45 .)

Include
The file cppad/rosen_34.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^5 )$$ where $h = (tf - ti) / M$ is the step size. If xf contains not a number nan , see the discussion of f .

Fun
The class Fun and the object F satisfy the prototype       Fun &F  This must support the following set of calls       F.Ode(t, x, f)      F.Ode_ind(t, x, f_t)      F.Ode_dep(t, x, f_x) 
t
In all three cases, the argument t has prototype       const Scalar &t  (see description of Scalar below).

x
In all three cases, the argument x 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)$ (see F(t, x) in Description ).

f_t
The argument f_t to F.Ode_ind has prototype       Vector &f_t  On input and output, f_t is a vector of size n and the input values of the elements of f_t do not matter. On output, the i-th element of f_t is set equal to $\partial_t F_i (t, x)$ (see F(t, x) in Description ).

f_x
The argument f_x to F.Ode_dep has prototype       Vector &f_x  On input and output, f_x is a vector of size n*n and the input values of the elements of f_x do not matter. On output, the [ i*n+j ] element of f_x is set equal to $\partial_{x(j)} F_i (t, x)$ (see F(t, x) in Description ).

Nan
If any of the elements of f , f_t , or f_x have the value not a number nan, the routine Rosen34 returns with all the elements of xf and e equal to nan.

Warning
The arguments f , f_t , and f_x must have a call by reference in their prototypes; i.e., do not forget the & in the prototype for f , f_t and f_x .

Optimization
Every call of the form       F.Ode_ind(t, x, f_t)  is directly followed by a call of the form       F.Ode_dep(t, x, f_x)  where the arguments t and x have not changed between calls. In many cases it is faster to compute the values of f_t and f_x together and then pass them back one at a time.

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^4 )$$ where $h = (tf - ti) / M$ is the step size.

Scalar
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 operations must be defined for Scalar objects a and b :
 Operation Description a < b less than operator (returns a bool object)

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 Rosen34 must not be parallel execution mode.

Example
The file rosen_34.cpp contains an example and test a test of using this routine. It returns true if it succeeds and false otherwise.

Source Code
The source code for this routine is in the file cppad/rosen_34.hpp.
Input File: cppad/utility/rosen_34.hpp