Prev Next runge_45.cpp Headings

@(@\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 .
Runge45: Example and Test
Define @(@ X : \B{R} \times \B{R} \rightarrow \B{R}^n @)@ by @[@ X_j (b, t) = b \left( \sum_{k=0}^j t^k / k ! \right) @]@ for @(@ j = 0 , \ldots , n-1 @)@. It follows that @[@ \begin{array}{rcl} X_j (b, 0) & = & b \\ \partial_t X_j (b, t) & = & b \left( \sum_{k=0}^{j-1} t^k / k ! \right) \\ \partial_t X_j (b, t) & = & \left\{ \begin{array}{ll} 0 & {\rm if} \; j = 0 \\ X_{j-1} (b, t) & {\rm otherwise} \end{array} \right. \end{array} @]@ For a fixed @(@ t_f @)@, we can use Runge45 to define @(@ f : \B{R} \rightarrow \B{R}^n @)@ as an approximation for @(@ f(b) = X(b, t_f ) @)@. We can then compute @(@ f^{(1)} (b) @)@ which is an approximation for @[@ \partial_b X(b, t_f ) = \sum_{k=0}^j t_f^k / k ! @]@ 

# include <cstddef>              // for size_t
# include <limits>               // for machine epsilon
# include <cppad/cppad.hpp>      // for all of CppAD

namespace {

    template <class Scalar>
    class Fun {
    public:
        // constructor
        Fun(void)
        { }

        // set return value to X'(t)
        void Ode(
            const Scalar                    &t,
            const CPPAD_TESTVECTOR(Scalar) &x,
            CPPAD_TESTVECTOR(Scalar)       &f)
        {   size_t n  = x.size();
            f[0]      = 0.;
            for(size_t k = 1; k < n; k++)
                f[k] = x[k-1];
        }
    };
}

bool runge_45(void)
{   typedef CppAD::AD<double> Scalar;
    using CppAD::NearEqual;

    bool ok = true;     // initial return value
    size_t j;           // temporary indices

    size_t     n = 5;   // number components in X(t) and order of method
    size_t     M = 2;   // number of Runge45 steps in [ti, tf]
    Scalar ad_ti = 0.;  // initial time
    Scalar ad_tf = 2.;  // final time

    // value of independent variable at which to record operations
    CPPAD_TESTVECTOR(Scalar) ad_b(1);
    ad_b[0] = 1.;

    // declare b to be the independent variable
    Independent(ad_b);

    // object to evaluate ODE
    Fun<Scalar> ad_F;

    // xi = X(0)
    CPPAD_TESTVECTOR(Scalar) ad_xi(n);
    for(j = 0; j < n; j++)
        ad_xi[j] = ad_b[0];

    // compute Runge45 approximation for X(tf)
    CPPAD_TESTVECTOR(Scalar) ad_xf(n), ad_e(n);
    ad_xf = CppAD::Runge45(ad_F, M, ad_ti, ad_tf, ad_xi, ad_e);

    // stop recording and use it to create f : b -> xf
    CppAD::ADFun<double> f(ad_b, ad_xf);

    // evaluate f(b)
    CPPAD_TESTVECTOR(double)  b(1);
    CPPAD_TESTVECTOR(double) xf(n);
    b[0] = 1.;
    xf   = f.Forward(0, b);

    // check that f(b) = X(b, tf)
    double tf    = Value(ad_tf);
    double term  = 1;
    double sum   = 0;
    double eps   = 10. * CppAD::numeric_limits<double>::epsilon();
    for(j = 0; j < n; j++)
    {   sum += term;
        ok &= NearEqual(xf[j], b[0] * sum, eps, eps);
        term *= tf;
        term /= double(j+1);
    }

    // evalute f'(b)
    CPPAD_TESTVECTOR(double) d_xf(n);
    d_xf = f.Jacobian(b);

    // check that f'(b) = partial of X(b, tf) w.r.t b
    term  = 1;
    sum   = 0;
    for(j = 0; j < n; j++)
    {   sum += term;
        ok &= NearEqual(d_xf[j], sum, eps, eps);
        term *= tf;
        term /= double(j+1);
    }

    return ok;
}
Input File: example/utility/runge_45.cpp