------------------------------------------------------- lines 9-285 of file: include/cppad/utility/rosen_34.hpp ------------------------------------------------------- {xrst_begin Rosen34} {xrst_spell dep rosenbrock test test tf xf } A 3rd and 4th Order Rosenbrock ODE Solver ######################################### Syntax ****** | # ``include `` | *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 :ref:`Bib@Numerical Recipes` for a description of Rosenbrock ODE solvers). In particular, we use the formulas taken from page 100 of :ref:`Bib@Shampine, L.F.` (except that the fraction 98/108 has been correction to be 97/108). We use :math:`n` for the size of the vector *xi* . Let :math:`\B{R}` denote the real numbers and let :math:`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 :math:`X(tf)` where :math:`X : [ti , tf] \rightarrow \B{R}^n` is defined by the following initial value problem: .. math:: :nowrap: \begin{eqnarray} X(ti) & = & xi \\ X'(t) & = & F[t , X(t)] \end{eqnarray} If your set of ordinary differential equations are not stiff an explicit method may be better (perhaps :ref:`Runge45-name` .) Include ******* The file ``cppad/utility/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 :ref:`Rosen34@Vector` below). .. math:: X(tf) = xf + O( h^5 ) where :math:`h = (tf - ti) / M` is the step size. If *xf* contains not a number :ref:`nan-name` , see the discussion of :ref:`f` . Fun *** The class *Fun* and the object *F* satisfy the prototype *Fun* & *F* This must support the following set of calls | |tab| *F* . ``Ode`` ( *t* , *x* , *f* ) | |tab| *F* . ``Ode_ind`` ( *t* , *x* , *f_t* ) | |tab| *F* . ``Ode_dep`` ( *t* , *x* , *f_x* ) t = In all three cases, the argument *t* has prototype ``const`` *Scalar* & *t* (see description of :ref:`Rosen34@Scalar` below). x = In all three cases, the argument *x* has prototype ``const`` *Vector* & *x* and has size *n* (see description of :ref:`Rosen34@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 :math:`F(t, x)` (see *F* ( *t* , *x* ) in :ref:`Rosen34@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 :math:`\partial_t F_i (t, x)` (see *F* ( *t* , *x* ) in :ref:`Rosen34@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 :math:`\partial_{x(j)} F_i (t, x)` (see *F* ( *t* , *x* ) in :ref:`Rosen34@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 :math:`h = (tf - ti) / M`, thus the larger *M* , the more accurate the return value *xf* is as an approximation for :math:`X(tf)`. ti ** The argument *ti* has prototype ``const`` *Scalar* & *ti* (see description of :ref:`Rosen34@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 :math:`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* .. math:: e = O( h^4 ) where :math:`h = (tf - ti) / M` is the step size. Scalar ****** The type *Scalar* must satisfy the conditions for a :ref:`NumericType-name` . The routine :ref:`CheckNumericType-name` 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* : .. list-table:: :widths: auto * - **Operation** - **Description** * - *a* < *b* - less than operator (returns a ``bool`` object) Vector ****** The type *Vector* must be a :ref:`SimpleVector-name` class with :ref:`elements of type Scalar` . The routine :ref:`CheckSimpleVector-name` will generate an error message if this is not the case. Parallel Mode ************* For each set of types :ref:`Rosen34@Scalar` , :ref:`Rosen34@Vector` , and :ref:`Rosen34@Fun` , the first call to ``Rosen34`` must not be :ref:`parallel` execution mode. Example ******* {xrst_toc_hidden example/utility/rosen_34.cpp } The file :ref:`rosen_34.cpp-name` contains an example and test a test of using this routine. Source Code *********** The source code for this routine is in the file ``cppad/rosen_34.hpp`` . {xrst_end Rosen34}