Taylor's Ode Solver: An Example and Test

taylor_ode.cpp

@(@\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}} }@)@ This is cppad-20221105 documentation. Here is a link to its current documentation . Taylor's Ode Solver: An Example and Test
Purpose
This example uses the method described in taylor_ode to solve an ODE.

ODE
The ODE is defined by @(@ y(0) = 0 @)@ and @(@ y^1 (t) = g[ y(t) ] @)@ where the function @(@ g : \B{R}^n \rightarrow \B{R}^n @)@ is defined by @[@ g(y) = \left( \begin{array}{c} 1 \\ y_1 \\ \vdots \\ y_{n-1} \end{array} \right) @]@ and the initial condition is @(@ z(0) = 0 @)@.

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) = \left( \begin{array}{c} t \\ t^2 / 2 \\ \vdots \\ t^n / n ! \end{array} \right) @]@


# include <cppad/cppad.hpp>

// =========================================================================
// define types for each level
namespace { // BEGIN empty namespace

    typedef CppAD::AD<double>          a_double;
    typedef CPPAD_TESTVECTOR(double)   d_vector;
    typedef CPPAD_TESTVECTOR(a_double) a_vector;

    a_vector ode(const a_vector y)
    {   size_t n = y.size();
        a_vector g(n);
        g[0] = 1;
        for(size_t k = 1; k < n; k++)
            g[k] = y[k-1];
        return g;
    }

}

// -------------------------------------------------------------------------
// use Taylor's method to solve this ordinary differential equaiton
bool taylor_ode(void)
{   // initialize the return value as true
    bool ok = true;

    // 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
    a_vector ay(n), ag(n);
    CppAD::Independent( ay );
    ag = ode(ay);
    CppAD::ADFun<double> g(ay, ag);

    // initialize the solution vector at time zero
    d_vector y(n);
    for(size_t j = 0; j < n; j++)
        y[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
    d_vector yk(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
        yk           = y;     // initialize y^{(k)}  for k = 0
        double dt_kp = dt;    // initialize dt^(k+1) for k = 0
        for(size_t k = 0; k < order; k++)
        {   // evaluate k-th order Taylor coefficient of z(t) = g(y(t))
            zk = g.Forward(k, yk);

            for(size_t j = 0; j < n; j++)
            {   // convert to (k+1)-Taylor coefficient for y
                yk[j] = zk[j] / double(k + 1);

                // add term for to this Taylor coefficient
                // to solution for y(t, x)
                y[j] += yk[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     = double(n_step) * dt;
    for(size_t i = 0; i < n; i++)
    {   check *= t / double(i + 1);
        ok &= CppAD::NearEqual(y[i], check, eps, eps);
    }

    return ok;
}

Input File: example/general/taylor_ode.cpp