Ignore:
Timestamp:
Aug 8, 2007 2:07:48 AM (13 years ago)
Author:
bradbell
Message:

trunk: Add nan to return specifications for Rosen34.

whats_new_07.omh: user's view of the changes.
example/ode_err_control.cpp: check x instead of f.
test_more/ode_err_control.cpp: test nan returns in Rosen34.
cppad/rosen_34.hpp: add nan return case.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/test_more/ode_err_control.cpp

    r982 r983  
    213213        return ok;
    214214}
    215 
     215/*
     216Define
     217$latex X : \R \rightarrow \R^2$$ by
     218$latex \[
     219\begin{array}{rcl}
     220        X_0 (0)       & = & 1  \\
     221        X_1 (0)       & = & 0  \\
     222        X_0^{(1)} (t) & = & 2 \alpha t X_0 (t)  \\
     223        X_1^{(1)} (t) & = &  1 / X_0 (t)
     224\end{array}
     225\] $$
     226It follows that
     227$latex \[
     228\begin{array}{rcl}
     229X_0 (t) & = &  \exp ( \alpha t^2 )  \\
     230X_1 (t) & = &  \int_0^t \exp( - \alpha s^2 ) {\bf d} s \\
     231& = & 
     232\frac{1}{ \sqrt{\alpha} \int_0^{\sqrt{\alpha} t} \exp( - r^2 ) {\bf d} r
     233\\
     234& = & \frac{\sqrt{\pi}}{ 2 \sqrt{\alpha} {\rm erf} ( \sqrt{\alpha} t )
     235\end{array}
     236\] $$
     237If $latex X_0 (t) < 0$$,
     238we return $code nan$$ in order to inform
     239$code OdeErrControl$$ that its is taking to large a step.
     240
     241*/
     242
     243# include <cppad/rosen_34.hpp>          // CppAD::Rosen34
     244# include <cppad/cppad.h>
     245
     246namespace {
     247        // --------------------------------------------------------------
     248        class Fun_three {
     249        private:
     250                const double alpha_;
     251                bool was_negative_;
     252        public:
     253                // constructor
     254                Fun_three(double alpha) : alpha_(alpha), was_negative_(false)
     255                { }
     256
     257                // set f = x'(t)
     258                void Ode(
     259                        const double                &t,
     260                        const CppAD::vector<double> &x,
     261                        CppAD::vector<double>       &f)
     262                {       f[0] = 2. * alpha_ * t * x[0];
     263                        f[1] = 1. / x[0];       
     264                        // case where ODE does not make sense
     265                        if( x[0] < 0. || x[1] < 0. )
     266                        {       was_negative_ = true;
     267                                f[0] = CppAD::nan(0.);
     268                        }
     269                }
     270                // set f_t = df / dt
     271                void Ode_ind(
     272                        const double                &t,
     273                        const CppAD::vector<double> &x,
     274                        CppAD::vector<double>       &f_t)
     275                {
     276                        f_t[0] =  2. * alpha_ * x[0];
     277                        f_t[1] = 0.;   
     278                        if( x[0] < 0. || x[1] < 0. )
     279                        {       was_negative_ = true;
     280                                f_t[0] = CppAD::nan(0.);
     281                        }
     282                }
     283                // set f_x = df / dx
     284                void Ode_dep(
     285                        const double                &t,
     286                        const CppAD::vector<double> &x,
     287                        CppAD::vector<double>       &f_x)
     288                {       double x0_sq = x[0] * x[0];
     289                        f_x[0 * 2 + 0] = 2. * alpha_ * t;   // f0 w.r.t. x0
     290                        f_x[0 * 2 + 1] = 0.;                // f0 w.r.t. x1
     291                        f_x[1 * 2 + 0] = -1./x0_sq;         // f1 w.r.t. x0
     292                        f_x[1 * 2 + 1] = 0.;                // f1 w.r.t. x1
     293                        if( x[0] < 0. || x[1] < 0. )
     294                        {       was_negative_ = true;
     295                                f_x[0] = CppAD::nan(0.);
     296                        }
     297                }
     298                bool was_negative(void)
     299                {       return was_negative_; }
     300
     301        };
     302
     303        // --------------------------------------------------------------
     304        class Method_three {
     305        public:
     306                Fun_three F;
     307
     308                // constructor
     309                Method_three(double alpha) : F(alpha)
     310                { }
     311                void step(
     312                        double ta,
     313                        double tb,
     314                        CppAD::vector<double> &xa ,
     315                        CppAD::vector<double> &xb ,
     316                        CppAD::vector<double> &eb )
     317                {       xb = CppAD::Rosen34(F, 1, ta, tb, xa, eb);
     318                }
     319                size_t order(void)
     320                {       return 3; }
     321        };
     322}
     323
     324bool OdeErrControl_three(void)
     325{       bool ok = true;     // initial return value
     326
     327        double alpha = 10.;
     328        Method_three method(alpha);
     329
     330        CppAD::vector<double> xi(2);
     331        xi[0] = 1.;
     332        xi[1] = 0.;
     333
     334        CppAD::vector<double> eabs(2);
     335        eabs[0] = 1e-4;
     336        eabs[1] = 1e-4;
     337
     338        // inputs
     339        double ti   = 0.;
     340        double tf   = 1.;
     341        double smin = 1e-4;
     342        double smax = 1.;
     343        double scur = 1.;
     344        double erel = 0.;
     345
     346        // outputs
     347        CppAD::vector<double> ef(2);
     348        CppAD::vector<double> xf(2);
     349        CppAD::vector<double> maxabs(2);
     350        size_t nstep;
     351
     352       
     353        xf = OdeErrControl(method,
     354                ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);
     355
     356
     357        double x0       = exp( alpha * tf * tf );
     358        ok &= CppAD::NearEqual(x0, xf[0], 1e-4, 1e-4);
     359        ok &= CppAD::NearEqual(0., ef[0], 1e-4, 1e-4);
     360
     361        double root_pi    = sqrt( 4. * atan(1.));
     362        double root_alpha = sqrt( alpha );
     363        double x1 = CppAD::erf(alpha * tf) * root_pi / (2 * root_alpha);
     364        ok &= CppAD::NearEqual(x1, xf[1], 1e-4, 1e-4);
     365        ok &= CppAD::NearEqual(0., ef[1], 1e-4, 1e-4);
     366
     367        ok &= method.F.was_negative();
     368
     369        return ok;
     370}
     371
     372// ==========================================================================
    216373bool OdeErrControl(void)
    217374{       bool ok = true;
    218375        ok     &= OdeErrControl_one();
    219376        ok     &= OdeErrControl_two();
     377        ok     &= OdeErrControl_three();
    220378        return ok;
    221379}
Note: See TracChangeset for help on using the changeset viewer.