$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
Taylor's Ode Solver: A Multi-Level Adolc Example and Test

Purpose
This is a realistic example using two levels of AD; see mul_level . The first level uses Adolc's adouble type to tape the solution of an ordinary differential equation. This solution is then differentiated with respect to a parameter vector. The second level uses CppAD's type AD<adouble> to take derivatives during the solution of the differential equation. These derivatives are used in the application of Taylor's method to the solution of the ODE. The example mul_level_ode.cpp computes the same values using AD<double> and AD< AD<double> >. The example ode_taylor.cpp is a simpler applications of Taylor's method for solving an ODE.

ODE
For this example the ODE's are defined by the function $h : \B{R}^n \times \B{R}^n \rightarrow \B{R}^n$ where $$h[ x, y(t, x) ] = \left( \begin{array}{c} x_0 \\ x_1 y_0 (t, x) \\ \vdots \\ x_{n-1} y_{n-2} (t, x) \end{array} \right) = \left( \begin{array}{c} \partial_t y_0 (t , x) \\ \partial_t y_1 (t , x) \\ \vdots \\ \partial_t y_{n-1} (t , x) \end{array} \right)$$ and the initial condition $y(0, x) = 0$. The value of $x$ is fixed during the solution of the ODE and the function $g : \B{R}^n \rightarrow \B{R}^n$ is used to define the ODE where $$g(y) = \left( \begin{array}{c} x_0 \\ x_1 y_0 \\ \vdots \\ x_{n-1} y_{n-2} \end{array} \right)$$

ODE Solution
The solution for this example can be calculated by starting with the first row and then using the solution for the first row to solve the second and so on. Doing this we obtain $$y(t, x ) = \left( \begin{array}{c} x_0 t \\ x_1 x_0 t^2 / 2 \\ \vdots \\ x_{n-1} x_{n-2} \ldots x_0 t^n / n ! \end{array} \right)$$

Derivative of ODE Solution
Differentiating the solution above, with respect to the parameter vector $x$, we notice that $$\partial_x y(t, x ) = \left( \begin{array}{cccc} y_0 (t,x) / x_0 & 0 & \cdots & 0 \\ y_1 (t,x) / x_0 & y_1 (t,x) / x_1 & 0 & \vdots \\ \vdots & \vdots & \ddots & 0 \\ y_{n-1} (t,x) / x_0 & y_{n-1} (t,x) / x_1 & \cdots & y_{n-1} (t,x) / x_{n-1} \end{array} \right)$$

An m-th order Taylor method for approximating the solution of an ordinary differential equations is $$y(t + \Delta t , x) \approx \sum_{k=0}^p \partial_t^k y(t , x ) \frac{ \Delta t^k }{ k ! } = y^{(0)} (t , x ) + y^{(1)} (t , x ) \Delta t + \cdots + y^{(p)} (t , x ) \Delta t^p$$ where the Taylor coefficients $y^{(k)} (t, x)$ are defined by $$y^{(k)} (t, x) = \partial_t^k y(t , x ) / k !$$ We define the function $z(t, x)$ by the equation $$z ( t , x ) = g[ y ( t , x ) ] = h [ x , y( t , x ) ]$$ It follows that $$\begin{array}{rcl} \partial_t y(t, x) & = & z (t , x) \\ \partial_t^{k+1} y(t , x) & = & \partial_t^k z (t , x) \\ y^{(k+1)} ( t , x) & = & z^{(k)} (t, x) / (k+1) \end{array}$$ where $z^{(k)} (t, x)$ is the k-th order Taylor coefficient for $z(t, x)$. In the example below, the Taylor coefficients $$y^{(0)} (t , x) , \ldots , y^{(k)} ( t , x )$$ are used to calculate the Taylor coefficient $z^{(k)} ( t , x )$ which in turn gives the value for $y^{(k+1)} y ( t , x)$.

The file base_adolc.hpp is implements the Base type requirements where Base is adolc.
Source  // suppress conversion warnings before other includes # include <cppad/wno_conversion.hpp> // # include <adolc/adouble.h> # include <adolc/taping.h> # include <adolc/drivers/drivers.h> // definitions not in Adolc distribution and required to use CppAD::AD<adouble> # include <cppad/example/base_adolc.hpp> # include <cppad/cppad.hpp> // ========================================================================== namespace { // BEGIN empty namespace // define types for each level typedef adouble 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_adolc( 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), Z(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 using a2type Z = G(a2y); // define differentiable version of g: X -> Y // that computes its derivatives using a1type CppAD::ADFun<a1type> a1g(a2y, Z); // Use Taylor's method to take a step a1y_k = a1y; // initialize y^{(k)} a1type dt_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] * dt_kp; } // next power of t dt_kp *= a1dt; // next Taylor coefficient a1y_k = a1y_kp; } } return a1y; } } // END empty namespace // ========================================================================== // Routine that tests algorithmic differentiation of solutions computed // by the routine taylor_ode. bool mul_level_adolc_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; // set up for thread_alloc memory allocator using CppAD::thread_alloc; // the allocator size_t capacity; // capacity of an allocation // the vector x with length n (or greater) in double double* x = thread_alloc::create_array<double>(n, capacity); // the vector x with length n in a1type 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 int tag = 0; // Adolc setup int keep = 1; trace_on(tag, keep); for(i = 0; i < n; i++) a1x[i] <<= double(i + 1); // a1x is independent for adouble type // arguments to taylor_ode_adolc 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 = 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_adolc(G, order, nstep, a1dt, a1y_ini); // declare the differentiable fucntion f : x -> y_final // (corresponding to the tape of adouble operations) double* y_final = thread_alloc::create_array<double>(n, capacity); for(i = 0; i < n; i++) a1y_final[i] >>= y_final[i]; trace_off(); // check function values double check = 1.; double t = nstep * a1dt.value(); for(i = 0; i < n; i++) { check *= x[i] * t / double(i + 1); ok &= CppAD::NearEqual(y_final[i], check, eps, eps); } // memory where Jacobian will be returned double* jac_ = thread_alloc::create_array<double>(n * n, capacity); double** jac = thread_alloc::create_array<double*>(n, capacity); for(i = 0; i < n; i++) jac[i] = jac_ + i * n; // evaluate Jacobian of h at a size_t m = n; // # dependent variables jacobian(tag, int(m), int(n), x, jac); // check Jacobian for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { if( i < j ) check = 0.; else check = y_final[i] / x[j]; ok &= CppAD::NearEqual(jac[i][j], check, eps, eps); } } // make memroy avaiable for other use by this thread thread_alloc::delete_array(x); thread_alloc::delete_array(y_final); thread_alloc::delete_array(jac_); thread_alloc::delete_array(jac); return ok; }