|
|
Prev | Next | ode_taylor.cpp |
\[
Z(t + \Delta t)
\approx
Z^{(0)} (t) +
\frac{ Z^{(1)} (t) }{ 1 !} \Delta t + \cdots +
\frac{ Z^{(p)} (t) }{ p !} ( \Delta t )^p )
\]
h : \B{R}^n \rightarrow \B{R}^n
,
which for this example is given by
\[
Z^{(1)} (t) = H[ Z(t) ] =
\left( \begin{array}{c}
1 \\
Z_1 (t) \\
\vdots \\
Z_{n-1} (t)
\end{array} \right)
\]
and the initial condition is
z(0) = 0
.
\[
Z(t) =
\left( \begin{array}{c}
t \\
t^2 / 2 \\
\vdots \\
t^n / n !
\end{array} \right)
\]
k = 0 , \ldots , K
\[
z^{(k)} = \frac{ Z^{(k)} }{ k !} (t)
\]
we note that
\[
\begin{array}{rcl}
Z^{(1)} (t)
& = &
H( z^{(0)} + z^{(1)} t + \cdots + z^{(K)} t^K ) + O( t^{K+1} )
\\
& = &
h^{(0)} + h^{(1)} t + \cdots + h^{(K)} t^K + O( t^{K+1} )
\end{array}
\]
where
h^{(k)}
is the k-th order Taylor coefficient for
H( Z(t) )
.
Taking K-th order derivatives of both sides we obtain
\[
\begin{array}{rcl}
Z^{(K+1)} (t) & = & K ! h^{(K)} \\
z^{(K+1)} & = & h^{(K)} / K
\end{array}
\]
The code below uses this relationship to implement Taylor's
method for approximating the solution of an ODE.
# include <cppad/cppad.hpp>
// =========================================================================
// define types for each level
namespace { // BEGIN empty namespace
using CppAD::AD;
CPPAD_TESTVECTOR( AD<double> ) ode(
const CPPAD_TESTVECTOR( AD<double> )& Z )
{ size_t n = Z.size();
CPPAD_TESTVECTOR( AD<double> ) y(n);
y[0] = 1;
for(size_t k = 1; k < n; k++)
y[k] = Z[k-1];
return y;
}
}
// -------------------------------------------------------------------------
// Example that uses Taylor's method to solve ordinary differential equaitons
bool ode_taylor(void)
{ // initialize the return value as true
bool ok = true;
// some temporary indices
size_t i, j, k;
// The ODE does not depend on the arugment values
// so only tape once, also note that ode does not depend on t
size_t n = 5; // number of independent and dependent variables
CPPAD_TESTVECTOR( AD<double> ) a_x(n), a_y(n);
CppAD::Independent( a_x );
a_y = ode(a_x);
CppAD::ADFun<double> H(a_x, a_y);
// initialize the solution vector at time zero
CPPAD_TESTVECTOR( double ) z(n);
for(j = 0; j < n; j++)
z[j] = 0.0;
size_t order = n; // order of the Taylor method
size_t n_step = 4; // number of time steps
double dt = 0.5; // step size in time
// Taylor coefficients of order k
CPPAD_TESTVECTOR( double ) hk(n), zk(n);
// loop with respect to each step of Taylor's method
for(size_t i_step = 0; i_step < n_step; i_step++)
{ // Use Taylor's method to take a step
zk = z; // initialize z^{(k)} for k = 0
double dt_kp = dt; // initialize dt^(k+1) for k = 0
for(k = 0; k < order; k++)
{ // evaluate k-th order Taylor coefficient of H
hk = H.Forward(k, zk);
for(j = 0; j < n; j++)
{ // convert to (k+1)-Taylor coefficient for z
zk[j] = hk[j] / double(k + 1);
// add term for to this Taylor coefficient
// to solution for y(t, x)
z[j] += zk[j] * dt_kp;
}
// next power of t
dt_kp *= dt;
}
}
// check solution of the ODE,
// Taylor's method should have no truncation error for this case
double eps = 100. * std::numeric_limits<double>::epsilon();
double check = 1.;
double t = n_step * dt;
for(i = 0; i < n; i++)
{ check *= t / double(i + 1);
ok &= CppAD::NearEqual(z[i], check, eps, eps);
}
return ok;
}