You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
428 lines
11 KiB
428 lines
11 KiB
# ifndef CPPAD_UTILITY_RUNGE_45_HPP
|
|
# define CPPAD_UTILITY_RUNGE_45_HPP
|
|
|
|
/* --------------------------------------------------------------------------
|
|
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
|
|
|
|
CppAD is distributed under multiple licenses. This distribution is under
|
|
the terms of the
|
|
Eclipse Public License Version 1.0.
|
|
|
|
A copy of this license is included in the COPYING file of this distribution.
|
|
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
|
|
-------------------------------------------------------------------------- */
|
|
|
|
/*
|
|
$begin Runge45$$
|
|
$spell
|
|
std
|
|
fabs
|
|
cppad.hpp
|
|
bool
|
|
xf
|
|
templated
|
|
const
|
|
Runge-Kutta
|
|
CppAD
|
|
xi
|
|
ti
|
|
tf
|
|
Karp
|
|
$$
|
|
|
|
|
|
$section An Embedded 4th and 5th Order Runge-Kutta ODE Solver$$
|
|
$mindex Runge45 Runge Kutta solve differential equation$$
|
|
|
|
$head Syntax$$
|
|
$codei%# include <cppad/utility/runge_45.hpp>
|
|
%$$
|
|
$icode%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%)
|
|
%$$
|
|
$icode%xf% = Runge45(%F%, %M%, %ti%, %tf%, %xi%, %e%)
|
|
%$$
|
|
|
|
|
|
$head Purpose$$
|
|
This is an implementation of the
|
|
Cash-Karp embedded 4th and 5th order Runge-Kutta ODE solver
|
|
described in Section 16.2 of $cref/Numerical Recipes/Bib/Numerical Recipes/$$.
|
|
We use $latex n$$ for the size of the vector $icode xi$$.
|
|
Let $latex \B{R}$$ denote the real numbers
|
|
and let $latex F : \B{R} \times \B{R}^n \rightarrow \B{R}^n$$
|
|
be a smooth function.
|
|
The return value $icode xf$$ contains a 5th order
|
|
approximation for the value $latex X(tf)$$ where
|
|
$latex X : [ti , tf] \rightarrow \B{R}^n$$ is defined by
|
|
the following initial value problem:
|
|
$latex \[
|
|
\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 $cref Rosen34$$.)
|
|
|
|
$head Operation Sequence$$
|
|
The $cref/operation sequence/glossary/Operation/Sequence/$$ for $icode Runge$$
|
|
does not depend on any of its $icode Scalar$$ input values provided that
|
|
the operation sequence for
|
|
$codei%
|
|
%F%.Ode(%t%, %x%, %f%)
|
|
%$$
|
|
does not on any of its $icode Scalar$$ inputs (see below).
|
|
|
|
$head Include$$
|
|
The file $code cppad/runge_45.hpp$$ is included by $code cppad/cppad.hpp$$
|
|
but it can also be included separately with out the rest of
|
|
the $code CppAD$$ routines.
|
|
|
|
$head xf$$
|
|
The return value $icode xf$$ has the prototype
|
|
$codei%
|
|
%Vector% %xf%
|
|
%$$
|
|
and the size of $icode xf$$ is equal to $icode n$$
|
|
(see description of $cref/Vector/Runge45/Vector/$$ below).
|
|
$latex \[
|
|
X(tf) = xf + O( h^6 )
|
|
\] $$
|
|
where $latex h = (tf - ti) / M$$ is the step size.
|
|
If $icode xf$$ contains not a number $cref nan$$,
|
|
see the discussion for $cref/f/Runge45/Fun/f/$$.
|
|
|
|
$head Fun$$
|
|
The class $icode Fun$$
|
|
and the object $icode F$$ satisfy the prototype
|
|
$codei%
|
|
%Fun% &%F%
|
|
%$$
|
|
The object $icode F$$ (and the class $icode Fun$$)
|
|
must have a member function named $code Ode$$
|
|
that supports the syntax
|
|
$codei%
|
|
%F%.Ode(%t%, %x%, %f%)
|
|
%$$
|
|
|
|
$subhead t$$
|
|
The argument $icode t$$ to $icode%F%.Ode%$$ has prototype
|
|
$codei%
|
|
const %Scalar% &%t%
|
|
%$$
|
|
(see description of $cref/Scalar/Runge45/Scalar/$$ below).
|
|
|
|
$subhead x$$
|
|
The argument $icode x$$ to $icode%F%.Ode%$$ has prototype
|
|
$codei%
|
|
const %Vector% &%x%
|
|
%$$
|
|
and has size $icode n$$
|
|
(see description of $cref/Vector/Runge45/Vector/$$ below).
|
|
|
|
$subhead f$$
|
|
The argument $icode f$$ to $icode%F%.Ode%$$ has prototype
|
|
$codei%
|
|
%Vector% &%f%
|
|
%$$
|
|
On input and output, $icode f$$ is a vector of size $icode n$$
|
|
and the input values of the elements of $icode f$$ do not matter.
|
|
On output,
|
|
$icode f$$ is set equal to $latex F(t, x)$$ in the differential equation.
|
|
If any of the elements of $icode f$$ have the value not a number $code nan$$
|
|
the routine $code Runge45$$ returns with all the
|
|
elements of $icode xf$$ and $icode e$$ equal to $code nan$$.
|
|
|
|
$subhead Warning$$
|
|
The argument $icode f$$ to $icode%F%.Ode%$$
|
|
must have a call by reference in its prototype; i.e.,
|
|
do not forget the $code &$$ in the prototype for $icode f$$.
|
|
|
|
$head M$$
|
|
The argument $icode M$$ has prototype
|
|
$codei%
|
|
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 $latex h = (tf - ti) / M$$, thus
|
|
the larger $icode M$$, the more accurate the
|
|
return value $icode xf$$ is as an approximation
|
|
for $latex X(tf)$$.
|
|
|
|
$head ti$$
|
|
The argument $icode ti$$ has prototype
|
|
$codei%
|
|
const %Scalar% &%ti%
|
|
%$$
|
|
(see description of $cref/Scalar/Runge45/Scalar/$$ below).
|
|
It specifies the initial time for $icode t$$ in the
|
|
differential equation; i.e.,
|
|
the time corresponding to the value $icode xi$$.
|
|
|
|
$head tf$$
|
|
The argument $icode tf$$ has prototype
|
|
$codei%
|
|
const %Scalar% &%tf%
|
|
%$$
|
|
It specifies the final time for $icode t$$ in the
|
|
differential equation; i.e.,
|
|
the time corresponding to the value $icode xf$$.
|
|
|
|
$head xi$$
|
|
The argument $icode xi$$ has the prototype
|
|
$codei%
|
|
const %Vector% &%xi%
|
|
%$$
|
|
and the size of $icode xi$$ is equal to $icode n$$.
|
|
It specifies the value of $latex X(ti)$$
|
|
|
|
$head e$$
|
|
The argument $icode e$$ is optional and has the prototype
|
|
$codei%
|
|
%Vector% &%e%
|
|
%$$
|
|
If $icode e$$ is present,
|
|
the size of $icode e$$ must be equal to $icode n$$.
|
|
The input value of the elements of $icode e$$ does not matter.
|
|
On output
|
|
it contains an element by element
|
|
estimated bound for the absolute value of the error in $icode xf$$
|
|
$latex \[
|
|
e = O( h^5 )
|
|
\] $$
|
|
where $latex h = (tf - ti) / M$$ is the step size.
|
|
If on output, $icode e$$ contains not a number $code nan$$,
|
|
see the discussion for $cref/f/Runge45/Fun/f/$$.
|
|
|
|
$head Scalar$$
|
|
The type $icode Scalar$$ must satisfy the conditions
|
|
for a $cref NumericType$$ type.
|
|
The routine $cref CheckNumericType$$ will generate an error message
|
|
if this is not the case.
|
|
|
|
$subhead fabs$$
|
|
In addition, the following function must be defined for
|
|
$icode Scalar$$ objects $icode a$$ and $icode b$$
|
|
$codei%
|
|
%a% = fabs(%b%)
|
|
%$$
|
|
Note that this operation is only used for computing $icode e$$; hence
|
|
the operation sequence for $icode xf$$ can still be independent of
|
|
the arguments to $code Runge45$$ even if
|
|
$codei%
|
|
fabs(%b%) = std::max(-%b%, %b%)
|
|
%$$.
|
|
|
|
$head Vector$$
|
|
The type $icode Vector$$ must be a $cref SimpleVector$$ class with
|
|
$cref/elements of type Scalar/SimpleVector/Elements of Specified Type/$$.
|
|
The routine $cref CheckSimpleVector$$ will generate an error message
|
|
if this is not the case.
|
|
|
|
$head Parallel Mode$$
|
|
For each set of types
|
|
$cref/Scalar/Runge45/Scalar/$$,
|
|
$cref/Vector/Runge45/Vector/$$, and
|
|
$cref/Fun/Runge45/Fun/$$,
|
|
the first call to $code Runge45$$
|
|
must not be $cref/parallel/ta_in_parallel/$$ execution mode.
|
|
|
|
|
|
$head Example$$
|
|
$children%
|
|
example/utility/runge45_1.cpp%
|
|
example/general/runge45_2.cpp
|
|
%$$
|
|
The file
|
|
$cref runge45_1.cpp$$
|
|
contains a simple example and test of $code Runge45$$.
|
|
It returns true if it succeeds and false otherwise.
|
|
$pre
|
|
|
|
$$
|
|
The file
|
|
$cref runge45_2.cpp$$ contains an example using $code Runge45$$
|
|
in the context of algorithmic differentiation.
|
|
It also returns true if it succeeds and false otherwise.
|
|
|
|
$head Source Code$$
|
|
The source code for this routine is in the file
|
|
$code cppad/runge_45.hpp$$.
|
|
|
|
$end
|
|
--------------------------------------------------------------------------
|
|
*/
|
|
# include <cstddef>
|
|
# include <cppad/core/cppad_assert.hpp>
|
|
# include <cppad/utility/check_simple_vector.hpp>
|
|
# include <cppad/utility/check_numeric_type.hpp>
|
|
# include <cppad/utility/nan.hpp>
|
|
|
|
// needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
|
|
# include <cppad/utility/thread_alloc.hpp>
|
|
|
|
namespace CppAD { // BEGIN CppAD namespace
|
|
|
|
template <typename Scalar, typename Vector, typename Fun>
|
|
Vector Runge45(
|
|
Fun &F ,
|
|
size_t M ,
|
|
const Scalar &ti ,
|
|
const Scalar &tf ,
|
|
const Vector &xi )
|
|
{ Vector e( xi.size() );
|
|
return Runge45(F, M, ti, tf, xi, e);
|
|
}
|
|
|
|
template <typename Scalar, typename Vector, typename Fun>
|
|
Vector Runge45(
|
|
Fun &F ,
|
|
size_t M ,
|
|
const Scalar &ti ,
|
|
const Scalar &tf ,
|
|
const Vector &xi ,
|
|
Vector &e )
|
|
{
|
|
CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
|
|
|
|
// check numeric type specifications
|
|
CheckNumericType<Scalar>();
|
|
|
|
// check simple vector class specifications
|
|
CheckSimpleVector<Scalar, Vector>();
|
|
|
|
// Cash-Karp parameters for embedded Runge-Kutta method
|
|
// are static to avoid recalculation on each call and
|
|
// do not use Vector to avoid possible memory leak
|
|
static Scalar a[6] = {
|
|
Scalar(0),
|
|
Scalar(1) / Scalar(5),
|
|
Scalar(3) / Scalar(10),
|
|
Scalar(3) / Scalar(5),
|
|
Scalar(1),
|
|
Scalar(7) / Scalar(8)
|
|
};
|
|
static Scalar b[5 * 5] = {
|
|
Scalar(1) / Scalar(5),
|
|
Scalar(0),
|
|
Scalar(0),
|
|
Scalar(0),
|
|
Scalar(0),
|
|
|
|
Scalar(3) / Scalar(40),
|
|
Scalar(9) / Scalar(40),
|
|
Scalar(0),
|
|
Scalar(0),
|
|
Scalar(0),
|
|
|
|
Scalar(3) / Scalar(10),
|
|
-Scalar(9) / Scalar(10),
|
|
Scalar(6) / Scalar(5),
|
|
Scalar(0),
|
|
Scalar(0),
|
|
|
|
-Scalar(11) / Scalar(54),
|
|
Scalar(5) / Scalar(2),
|
|
-Scalar(70) / Scalar(27),
|
|
Scalar(35) / Scalar(27),
|
|
Scalar(0),
|
|
|
|
Scalar(1631) / Scalar(55296),
|
|
Scalar(175) / Scalar(512),
|
|
Scalar(575) / Scalar(13824),
|
|
Scalar(44275) / Scalar(110592),
|
|
Scalar(253) / Scalar(4096)
|
|
};
|
|
static Scalar c4[6] = {
|
|
Scalar(2825) / Scalar(27648),
|
|
Scalar(0),
|
|
Scalar(18575) / Scalar(48384),
|
|
Scalar(13525) / Scalar(55296),
|
|
Scalar(277) / Scalar(14336),
|
|
Scalar(1) / Scalar(4),
|
|
};
|
|
static Scalar c5[6] = {
|
|
Scalar(37) / Scalar(378),
|
|
Scalar(0),
|
|
Scalar(250) / Scalar(621),
|
|
Scalar(125) / Scalar(594),
|
|
Scalar(0),
|
|
Scalar(512) / Scalar(1771)
|
|
};
|
|
|
|
CPPAD_ASSERT_KNOWN(
|
|
M >= 1,
|
|
"Error in Runge45: the number of steps is less than one"
|
|
);
|
|
CPPAD_ASSERT_KNOWN(
|
|
e.size() == xi.size(),
|
|
"Error in Runge45: size of e not equal to size of xi"
|
|
);
|
|
size_t i, j, k, m; // indices
|
|
|
|
size_t n = xi.size(); // number of components in X(t)
|
|
Scalar ns = Scalar(int(M)); // number of steps as Scalar object
|
|
Scalar h = (tf - ti) / ns; // step size
|
|
Scalar zero_or_nan = Scalar(0); // zero (nan if Ode returns has a nan)
|
|
for(i = 0; i < n; i++) // initialize e = 0
|
|
e[i] = zero_or_nan;
|
|
|
|
// vectors used to store values returned by F
|
|
Vector fh(6 * n), xtmp(n), ftmp(n), x4(n), x5(n), xf(n);
|
|
|
|
xf = xi; // initialize solution
|
|
for(m = 0; m < M; m++)
|
|
{ // time at beginning of this interval
|
|
// (convert to int to avoid MS compiler warning)
|
|
Scalar t = ti * (Scalar(int(M - m)) / ns)
|
|
+ tf * (Scalar(int(m)) / ns);
|
|
|
|
// loop over integration steps
|
|
x4 = x5 = xf; // start x4 and x5 at same point for each step
|
|
for(j = 0; j < 6; j++)
|
|
{ // loop over function evaluations for this step
|
|
xtmp = xf; // location for next function evaluation
|
|
for(k = 0; k < j; k++)
|
|
{ // loop over previous function evaluations
|
|
Scalar bjk = b[ (j-1) * 5 + k ];
|
|
for(i = 0; i < n; i++)
|
|
{ // loop over elements of x
|
|
xtmp[i] += bjk * fh[i * 6 + k];
|
|
}
|
|
}
|
|
// ftmp = F(t + a[j] * h, xtmp)
|
|
F.Ode(t + a[j] * h, xtmp, ftmp);
|
|
|
|
// if ftmp has a nan, set zero_or_nan to nan
|
|
for(i = 0; i < n; i++)
|
|
zero_or_nan *= ftmp[i];
|
|
|
|
for(i = 0; i < n; i++)
|
|
{ // loop over elements of x
|
|
Scalar fhi = ftmp[i] * h;
|
|
fh[i * 6 + j] = fhi;
|
|
x4[i] += c4[j] * fhi;
|
|
x5[i] += c5[j] * fhi;
|
|
x5[i] += zero_or_nan;
|
|
}
|
|
}
|
|
// accumulate error bound
|
|
for(i = 0; i < n; i++)
|
|
{ // cant use abs because cppad.hpp may not be included
|
|
Scalar diff = x5[i] - x4[i];
|
|
e[i] += fabs(diff);
|
|
e[i] += zero_or_nan;
|
|
}
|
|
|
|
// advance xf for this step using x5
|
|
xf = x5;
|
|
}
|
|
return xf;
|
|
}
|
|
|
|
} // End CppAD namespace
|
|
|
|
# endif
|
|
|