# include <cppad/cppad.hpp>
// =========================================================================
// define types for each level
namespace { // BEGIN empty namespace
typedef CppAD::AD<double> a1type;
typedef CppAD::AD<a1type> a2type;
// -------------------------------------------------------------------------
// class definition for C++ function object that defines ODE
class Ode {
private:
// copy of a that is set by constructor and used by g(y)
CPPAD_TESTVECTOR(a1type) a1x_;
public:
// constructor
Ode(const CPPAD_TESTVECTOR(a1type)& a1x) : a1x_(a1x)
{ }
// the function g(y) is evaluated with two levels of taping
CPPAD_TESTVECTOR(a2type) operator()
( const CPPAD_TESTVECTOR(a2type)& a2y) const
{ size_t n = a2y.size();
CPPAD_TESTVECTOR(a2type) a2g(n);
size_t i;
a2g[0] = a1x_[0];
for(i = 1; i < n; i++)
a2g[i] = a1x_[i] * a2y[i-1];
return a2g;
}
};
// -------------------------------------------------------------------------
// Routine that uses Taylor's method to solve ordinary differential equaitons
// and allows for algorithmic differentiation of the solution.
CPPAD_TESTVECTOR(a1type) taylor_ode(
Ode G , // function that defines the ODE
size_t order , // order of Taylor's method used
size_t nstep , // number of steps to take
const a1type& a1dt , // Delta t for each step
const CPPAD_TESTVECTOR(a1type)& a1y_ini) // y(t) at the initial time
{
// some temporary indices
size_t i, k, ell;
// number of variables in the ODE
size_t n = a1y_ini.size();
// copies of x and g(y) with two levels of taping
CPPAD_TESTVECTOR(a2type) a2y(n), a2z(n);
// y, y^{(k)} , z^{(k)}, and y^{(k+1)}
CPPAD_TESTVECTOR(a1type) a1y(n), a1y_k(n), a1z_k(n), a1y_kp(n);
// initialize x
for(i = 0; i < n; i++)
a1y[i] = a1y_ini[i];
// loop with respect to each step of Taylors method
for(ell = 0; ell < nstep; ell++)
{ // prepare to compute derivatives using a1type
for(i = 0; i < n; i++)
a2y[i] = a1y[i];
CppAD::Independent(a2y);
// evaluate ODE in a2type
a2z = G(a2y);
// define differentiable version of a1g: y -> z
// that computes its derivatives using a1type objects
CppAD::ADFun<a1type> a1g(a2y, a2z);
// Use Taylor's method to take a step
a1y_k = a1y; // initialize y^{(k)}
a1type a1dt_kp = a1dt; // initialize dt^(k+1)
for(k = 0; k <= order; k++)
{ // evaluate k-th order Taylor coefficient of y
a1z_k = a1g.Forward(k, a1y_k);
for(i = 0; i < n; i++)
{ // convert to (k+1)-Taylor coefficient for x
a1y_kp[i] = a1z_k[i] / a1type(k + 1);
// add term for to this Taylor coefficient
// to solution for y(t, x)
a1y[i] += a1y_kp[i] * a1dt_kp;
}
// next power of t
a1dt_kp *= a1dt;
// next Taylor coefficient
a1y_k = a1y_kp;
}
}
return a1y;
}
} // END empty namespace
// ==========================================================================
// Routine that tests alogirhtmic differentiation of solutions computed
// by the routine taylor_ode.
bool mul_level_ode(void)
{ bool ok = true;
double eps = 100. * std::numeric_limits<double>::epsilon();
// number of components in differential equation
size_t n = 4;
// some temporary indices
size_t i, j;
// parameter vector in both double and a1type
CPPAD_TESTVECTOR(double) x(n);
CPPAD_TESTVECTOR(a1type) a1x(n);
for(i = 0; i < n; i++)
a1x[i] = x[i] = double(i + 1);
// declare the parameters as the independent variable
CppAD::Independent(a1x);
// arguments to taylor_ode
Ode G(a1x); // function that defines the ODE
size_t order = n; // order of Taylor's method used
size_t nstep = 2; // number of steps to take
a1type a1dt = double(1.); // Delta t for each step
// value of y(t, x) at the initial time
CPPAD_TESTVECTOR(a1type) a1y_ini(n);
for(i = 0; i < n; i++)
a1y_ini[i] = 0.;
// integrate the differential equation
CPPAD_TESTVECTOR(a1type) a1y_final(n);
a1y_final = taylor_ode(G, order, nstep, a1dt, a1y_ini);
// define differentiable fucntion object f : x -> y_final
// that computes its derivatives in double
CppAD::ADFun<double> f(a1x, a1y_final);
// check function values
double check = 1.;
double t = double(nstep) * Value(a1dt);
for(i = 0; i < n; i++)
{ check *= x[i] * t / double(i + 1);
ok &= CppAD::NearEqual(Value(a1y_final[i]), check, eps, eps);
}
// evaluate the Jacobian of h at a
CPPAD_TESTVECTOR(double) jac ( f.Jacobian(x) );
// There appears to be a bug in g++ version 4.4.2 because it generates
// a warning for the equivalent form
// CPPAD_TESTVECTOR(double) jac = f.Jacobian(x);
// check Jacobian
for(i = 0; i < n; i++)
{ for(j = 0; j < n; j++)
{ double jac_ij = jac[i * n + j];
if( i < j )
check = 0.;
else check = Value( a1y_final[i] ) / x[j];
ok &= CppAD::NearEqual(jac_ij, check, eps, eps);
}
}
return ok;
}