Prev Next base2ad.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: base2ad Example and Test

See Also
taylor_ode.cpp , mul_level_ode.cpp

Purpose
This is a realistic example using base2ad to create an AD<Base> function from an Base function. The function represents an ordinary differential equation. It is differentiated with respect to its variables . These derivatives are used by the taylor_ode method. This solution is then differentiated with respect to the functions dynamic parameters .

ODE
For this example the function @(@ y : \B{R} \times \B{R}^n \rightarrow \B{R}^n @)@ is defined by @(@ y(0, x) = 0 @)@ and @(@ \partial_t y(t, x) = g(y, x) @)@ where @(@ g : \B{R}^n \times \B{R}^n \rightarrow \B{R}^n @)@ is defined by @[@ g(y, x) = \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) @]@


We define the function @(@ z(t, x) @)@ by the equation @[@ z ( t , x ) = g[ y ( t , x ), x ] @]@ see taylor_ode for the method used to compute the Taylor coefficients w.r.t @(@ t @)@ of @(@ y(t, x) @)@.

Source 


# include <cppad/cppad.hpp>

// =========================================================================
namespace { // BEGIN empty namespace

typedef CppAD::AD<double>                  a_double;

typedef CPPAD_TESTVECTOR(double)           d_vector;
typedef CPPAD_TESTVECTOR(a_double)         a_vector;

typedef CppAD::ADFun<double>               fun_double;
typedef CppAD::ADFun<a_double, double>     afun_double;

// -------------------------------------------------------------------------
// class definition for C++ function object that defines ODE
class Ode {
private:
    // copy of x that is set by constructor and used by g(y)
    a_vector x_;
public:
    // constructor
    Ode(const a_vector& x) : x_(x)
    { }
    // the function g(y) given the parameter vector x
    a_vector operator() (const a_vector& y) const
    {   size_t n = y.size();
        a_vector g(n);
        g[0] = x_[0];
        for(size_t i = 1; i < n; i++)
            g[i] = x_[i] * y[i-1];
        //
        return g;
    }
};

// -------------------------------------------------------------------------
// Routine that uses Taylor's method to solve ordinary differential equaitons
a_vector taylor_ode(
    afun_double&     fun_g   ,  // function that defines the ODE
    size_t           order   ,  // order of Taylor's method used
    size_t           nstep   ,  // number of steps to take
    const a_double&  dt      ,  // Delta t for each step
    const a_vector&  y_ini)     // y(t) at the initial time
{
    // number of variables in the ODE
    size_t n = y_ini.size();

    // initialize y
    a_vector y = y_ini;

    // loop with respect to each step of Taylors method
    for(size_t s = 0; s < nstep; s++)
    {
        // initialize
        a_vector y_k   = y;
        a_double dt_k  = a_double(1.0);
        a_vector next  = y;

        for(size_t k = 0; k < order; k++)
        {
            // evaluate k-th order Taylor coefficient z^{(k)} (t)
            a_vector z_k = fun_g.Forward(k, y_k);

            // dt^{k+1}
            dt_k *= dt;

            // y^{(k+1)}
            for(size_t i = 0; i < n; i++)
            {   // y^{(k+1)}
                y_k[i] = z_k[i] / a_double(k + 1);

                // add term for k+1 Taylor coefficient
                // to solution for next y
                next[i] += y_k[i] * dt_k;
            }
        }

        // take step
        y = next;
    }
    return y;
}
} // END empty namespace

// ==========================================================================
// Routine that tests alogirhtmic differentiation of solutions computed
// by the routine taylor_ode.
bool base2ad(void)
{   bool ok = true;
    double eps = 100. * std::numeric_limits<double>::epsilon();

    // number of components in differential equation
    size_t n = 4;

    // record function g(y, x)
    // with y as the independent variables and x as dynamic parameters
    a_vector  ay(n), ax(n);
    for(size_t i = 0; i < n; i++)
        ay[i] = ax[i] = double(i + 1);
    CppAD::Independent(ay, ax);

    // fun_g
    Ode G(ax);
    a_vector ag = G(ay);
    fun_double fun_g(ay, ag);


    // afun_g
    afun_double afun_g( fun_g.base2ad() ); // differential equation

    // other arguments to taylor_ode
    size_t   order = n;       // order of Taylor's method used
    size_t   nstep = 2;       // number of steps to take
    a_double adt   = 1.;      // Delta t for each step
    a_vector ay_ini(n);       // initial value of y
    for(size_t i = 0; i < n; i++)
        ay_ini[i] = 0.;

    // declare x as independent variables
    CppAD::Independent(ax);

    // the independent variables if this function are
    // the dynamic parameters in afun_g
    afun_g.new_dynamic(ax);

    // integrate the differential equation
    a_vector ay_final;
    ay_final = taylor_ode(afun_g, order, nstep, adt, ay_ini);

    // define differentiable fucntion object f(x) = y_final(x)
    // that computes its derivatives in double
    CppAD::ADFun<double> fun_f(ax, ay_final);

    // double version of ax
    d_vector x(n);
    for(size_t i = 0; i < n; i++)
        x[i] = Value( ax[i] );

    // check function values
    double check = 1.;
    double t     = double(nstep) * Value(adt);
    for(size_t i = 0; i < n; i++)
    {   check *= x[i] * t / double(i + 1);
        ok &= CppAD::NearEqual(Value(ay_final[i]), check, eps, eps);
    }

    // There appears to be a bug in g++ version 4.4.2 because it generates
    // a warning for the equivalent form
    // d_vector jac = fun_f.Jacobian(x);
    d_vector jac ( fun_f.Jacobian(x) );

    // check Jacobian
    for(size_t i = 0; i < n; i++)
    {   for(size_t j = 0; j < n; j++)
        {   double jac_ij = jac[i * n + j];
            if( i < j )
                check = 0.;
            else
                check = Value( ay_final[i] ) / x[j];
            ok &= CppAD::NearEqual(jac_ij, check, eps, eps);
        }
    }
    return ok;
}
Input File: example/general/base2ad.cpp