A Stiff Ode: Example and Test

ode_stiff.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 . A Stiff Ode: Example and Test Define @(@ x : \B{R} \rightarrow \B{R}^2 @)@ by @[@ \begin{array}{rcl} x_0 (0) & = & 1 \\ x_1 (0) & = & 0 \\ x_0^\prime (t) & = & - a_0 x_0 (t) \\ x_1^\prime (t) & = & + a_0 x_0 (t) - a_1 x_1 (t) \end{array} @]@ If @(@ a_0 \gg a_1 > 0 @)@, this is a stiff Ode and the analytic solution is @[@ \begin{array}{rcl} x_0 (t) & = & \exp( - a_0 t ) \\ x_1 (t) & = & a_0 [ \exp( - a_1 t ) - \exp( - a_0 t ) ] / ( a_0 - a_1 ) \end{array} @]@ The example tests Rosen34 using the relations above:


# include <cppad/cppad.hpp>

// To print the comparison, change the 0 to 1 on the next line.
# define CPPAD_ODE_STIFF_PRINT 0

namespace {
    // --------------------------------------------------------------
    class Fun {
    private:
        CPPAD_TESTVECTOR(double) a;
    public:
        // constructor
        Fun(const CPPAD_TESTVECTOR(double)& a_) : a(a_)
        { }
        // compute f(t, x)
        void Ode(
            const double                    &t,
            const CPPAD_TESTVECTOR(double) &x,
            CPPAD_TESTVECTOR(double)       &f)
        {   f[0]  = - a[0] * x[0];
            f[1]  = + a[0] * x[0] - a[1] * x[1];
        }
        // compute partial of f(t, x) w.r.t. t
        void Ode_ind(
            const double                    &t,
            const CPPAD_TESTVECTOR(double) &x,
            CPPAD_TESTVECTOR(double)       &f_t)
        {   f_t[0] = 0.;
            f_t[1] = 0.;
        }
        // compute partial of f(t, x) w.r.t. x
        void Ode_dep(
            const double                    &t,
            const CPPAD_TESTVECTOR(double) &x,
            CPPAD_TESTVECTOR(double)       &f_x)
        {   f_x[0] = -a[0];
            f_x[1] = 0.;
            f_x[2] = +a[0];
            f_x[3] = -a[1];
        }
    };
    // --------------------------------------------------------------
    class RungeMethod {
    private:
        Fun F;
    public:
        // constructor
        RungeMethod(const CPPAD_TESTVECTOR(double) &a_) : F(a_)
        { }
        void step(
            double                     ta ,
            double                     tb ,
            CPPAD_TESTVECTOR(double) &xa ,
            CPPAD_TESTVECTOR(double) &xb ,
            CPPAD_TESTVECTOR(double) &eb )
        {   xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
        }
        size_t order(void)
        {   return 5; }
    };
    class RosenMethod {
    private:
        Fun F;
    public:
        // constructor
        RosenMethod(const CPPAD_TESTVECTOR(double) &a_) : F(a_)
        { }
        void step(
            double                     ta ,
            double                     tb ,
            CPPAD_TESTVECTOR(double) &xa ,
            CPPAD_TESTVECTOR(double) &xb ,
            CPPAD_TESTVECTOR(double) &eb )
        {   xb = CppAD::Rosen34(F, 1, ta, tb, xa, eb);
        }
        size_t order(void)
        {   return 4; }
    };
}

bool OdeStiff(void)
{   bool ok = true;     // initial return value

    CPPAD_TESTVECTOR(double) a(2);
    a[0] = 1e3;
    a[1] = 1.;
    RosenMethod rosen(a);
    RungeMethod runge(a);
    Fun          gear(a);

    CPPAD_TESTVECTOR(double) xi(2);
    xi[0] = 1.;
    xi[1] = 0.;

    CPPAD_TESTVECTOR(double) eabs(2);
    eabs[0] = 1e-6;
    eabs[1] = 1e-6;

    CPPAD_TESTVECTOR(double) ef(2);
    CPPAD_TESTVECTOR(double) xf(2);
    CPPAD_TESTVECTOR(double) maxabs(2);
    size_t                nstep;

    size_t k;
    for(k = 0; k < 3; k++)
    {
        size_t M    = 5;
        double ti   = 0.;
        double tf   = 1.;
        double smin = 1e-7;
        double sini = 1e-7;
        double smax = 1.;
        double scur = .5;
        double erel = 0.;

        if( k == 0 )
        {   xf = CppAD::OdeErrControl(rosen, ti, tf,
            xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);
        }
        else if( k == 1 )
        {   xf = CppAD::OdeErrControl(runge, ti, tf,
            xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);
        }
        else if( k == 2 )
        {   xf = CppAD::OdeGearControl(gear, M, ti, tf,
            xi, smin, smax, sini, eabs, erel, ef, maxabs, nstep);
        }
        double x0 = exp(-a[0]*tf);
        ok &= CppAD::NearEqual(x0, xf[0], 0., eabs[0]);
        ok &= CppAD::NearEqual(0., ef[0], 0., eabs[0]);

        double x1 = a[0] *
            (exp(-a[1]*tf) - exp(-a[0]*tf))/(a[0] - a[1]);
        ok &= CppAD::NearEqual(x1, xf[1], 0., eabs[1]);
        ok &= CppAD::NearEqual(0., ef[1], 0., eabs[0]);
# if CPPAD_ODE_STIFF_PRINT
        const char* method[]={ "Rosen34", "Runge45", "Gear5" };
        std::cout << std::endl;
        std::cout << "method     = " << method[k] << std::endl;
        std::cout << "nstep      = " << nstep  << std::endl;
        std::cout << "x0         = " << x0 << std::endl;
        std::cout << "xf[0]      = " << xf[0] << std::endl;
        std::cout << "x0 - xf[0] = " << x0 - xf[0] << std::endl;
        std::cout << "ef[0]      = " << ef[0] << std::endl;
        std::cout << "x1         = " << x1 << std::endl;
        std::cout << "xf[1]      = " << xf[1] << std::endl;
        std::cout << "x1 - xf[1] = " << x1 - xf[1] << std::endl;
        std::cout << "ef[1]      = " << ef[1] << std::endl;
# endif
    }

    return ok;
}

Input File: example/general/ode_stiff.cpp