Changeset 982 for trunk/test_more


Ignore:
Timestamp:
Aug 7, 2007 9:18:27 AM (13 years ago)
Author:
bradbell
Message:

trunk: Add nan to return specifications for ode eqautions.

nan.hpp: routines for returning and detecting nan.
example/ode_err_control.cpp: a case that uses nan in ODE return.
text_more/ode_err_control.cpp: move old ode_err_control example to here.
makefile.am: add nan.hpp to distribution.
library.omh: add nan routines to library.
whats_new_07.omh: user's view of the changes.
fun_construct.hpp: fix call to Dependent in documentation.
track_new_del.hpp: improve TrackDelVec? error message.
ode_err_control.hpp: add use of nan for backing up and using smaller steps.
runge_45.hpp: use nan to communicate when step is too large.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/test_more/ode_err_control.cpp

    r683 r982  
    1 // BEGIN SHORT COPYRIGHT
    21/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-06 Bradley M. Bell
     2CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
    43
    54CppAD is distributed under multiple licenses. This distribution is under
     
    109Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
    1110-------------------------------------------------------------------------- */
    12 // END SHORT COPYRIGHT
    13 
    14 /*
     11
     12# include <cstddef>                    // for size_t
     13# include <cmath>                      // for exp
     14# include <cppad/ode_err_control.hpp>  // CppAD::OdeErrControl
     15# include <cppad/near_equal.hpp>       // CppAD::NearEqual
     16# include <cppad/vector.hpp>           // CppAD::vector
     17# include <cppad/runge_45.hpp>         // CppAD::Runge45
     18
     19/* -------------------------------------------------------------------------
    1520Test relative error with zero initial conditions.
    1621(Uses minimum step size to integrate).
    1722*/
    18 // BEGIN PROGRAM
    19 
    20 # include <cstddef>                 // for size_t
    21 # include <cppad/ode_err_control.hpp>   // CppAD::OdeErrControl
    22 # include <cppad/near_equal.hpp>       // CppAD::NearEqual
    23 # include <cppad/vector.hpp>    // CppAD::vector
    24 # include <cppad/runge_45.hpp>         // CppAD::Runge45
    2523
    2624namespace {
    2725        // --------------------------------------------------------------
    28         class Fun {
     26        class Fun_one {
    2927        private:
    3028                 size_t n;   // dimension of the state space
    3129        public:
    3230                // constructor
    33                 Fun(size_t n_) : n(n_)
     31                Fun_one(size_t n_) : n(n_)
    3432                { }
    3533
     
    4846
    4947        // --------------------------------------------------------------
    50         class Method {
    51         private:
    52                 Fun F;
    53         public:
    54                 // constructor
    55                 Method(size_t n_) : F(n_)
     48        class Method_one {
     49        private:
     50                Fun_one F;
     51        public:
     52                // constructor
     53                Method_one(size_t n_) : F(n_)
    5654                { }
    5755                void step(
     
    6866}
    6967
    70 bool OdeErrControl(void)
     68bool OdeErrControl_one(void)
    7169{       bool   ok = true;     // initial return value
    7270
     
    7573       
    7674        // construct method for n component solution
    77         Method method(n);
     75        Method_one method(n);
    7876
    7977        // inputs to OdeErrControl
     
    111109}
    112110
    113 // END PROGRAM
     111/*
     112Old example new just a test
     113Define
     114$latex X : \R \rightarrow \R^2$$ by
     115$latex \[
     116\begin{array}{rcl}
     117X_0 (t) & = &  - \exp ( - w_0 t )  \\
     118X_1 (t) & = & \frac{w_0}{w_1 - w_0} [ \exp ( - w_0 t ) - \exp( - w_1 t )]
     119\end{array}
     120\] $$
     121It follows that $latex X_0 (0) = 1$$, $latex X_1 (0) = 0$$ and
     122$latex \[
     123\begin{array}{rcl}
     124        X_0^{(1)} (t) & = & - w_0 X_0 (t)  \\
     125        X_1^{(1)} (t) & = & + w_0 X_0 (t) - w_1 X_1 (t)
     126\end{array}
     127\] $$
     128*/
     129
     130namespace {
     131        // --------------------------------------------------------------
     132        class Fun_two {
     133        private:
     134                 CppAD::vector<double> w;
     135        public:
     136                // constructor
     137                Fun_two(const CppAD::vector<double> &w_) : w(w_)
     138                { }
     139
     140                // set f = x'(t)
     141                void Ode(
     142                        const double                &t,
     143                        const CppAD::vector<double> &x,
     144                        CppAD::vector<double>       &f)
     145                {       f[0] = - w[0] * x[0];
     146                        f[1] = + w[0] * x[0] - w[1] * x[1];     
     147                }
     148        };
     149
     150        // --------------------------------------------------------------
     151        class Method_two {
     152        private:
     153                Fun_two F;
     154        public:
     155                // constructor
     156                Method_two(const CppAD::vector<double> &w_) : F(w_)
     157                { }
     158                void step(
     159                        double ta,
     160                        double tb,
     161                        CppAD::vector<double> &xa ,
     162                        CppAD::vector<double> &xb ,
     163                        CppAD::vector<double> &eb )
     164                {       xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
     165                }
     166                size_t order(void)
     167                {       return 4; }
     168        };
     169}
     170
     171bool OdeErrControl_two(void)
     172{       bool ok = true;     // initial return value
     173
     174        CppAD::vector<double> w(2);
     175        w[0] = 10.;
     176        w[1] = 1.;
     177        Method_two method(w);
     178
     179        CppAD::vector<double> xi(2);
     180        xi[0] = 1.;
     181        xi[1] = 0.;
     182
     183        CppAD::vector<double> eabs(2);
     184        eabs[0] = 1e-4;
     185        eabs[1] = 1e-4;
     186
     187        // inputs
     188        double ti   = 0.;
     189        double tf   = 1.;
     190        double smin = 1e-4;
     191        double smax = 1.;
     192        double scur = .5;
     193        double erel = 0.;
     194
     195        // outputs
     196        CppAD::vector<double> ef(2);
     197        CppAD::vector<double> xf(2);
     198        CppAD::vector<double> maxabs(2);
     199        size_t nstep;
     200
     201       
     202        xf = OdeErrControl(method,
     203                ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);
     204
     205        double x0 = exp(-w[0]*tf);
     206        ok &= CppAD::NearEqual(x0, xf[0], 1e-4, 1e-4);
     207        ok &= CppAD::NearEqual(0., ef[0], 1e-4, 1e-4);
     208
     209        double x1 = w[0] * (exp(-w[0]*tf) - exp(-w[1]*tf))/(w[1] - w[0]);
     210        ok &= CppAD::NearEqual(x1, xf[1], 1e-4, 1e-4);
     211        ok &= CppAD::NearEqual(0., ef[1], 1e-4, 1e-4);
     212
     213        return ok;
     214}
     215
     216bool OdeErrControl(void)
     217{       bool ok = true;
     218        ok     &= OdeErrControl_one();
     219        ok     &= OdeErrControl_two();
     220        return ok;
     221}
     222
Note: See TracChangeset for help on using the changeset viewer.