# Changeset 983 for trunk/test_more/ode_err_control.cpp

Ignore:
Timestamp:
Aug 8, 2007 2:07:48 AM (13 years ago)
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.
 r982 return ok; } /* Define $latex X : \R \rightarrow \R^2$$by latex $\begin{array}{rcl} X_0 (0) & = & 1 \\ X_1 (0) & = & 0 \\ X_0^{(1)} (t) & = & 2 \alpha t X_0 (t) \\ X_1^{(1)} (t) & = & 1 / X_0 (t) \end{array}$$$ It follows that$latex $\begin{array}{rcl} X_0 (t) & = & \exp ( \alpha t^2 ) \\ X_1 (t) & = & \int_0^t \exp( - \alpha s^2 ) {\bf d} s \\ & = & \frac{1}{ \sqrt{\alpha} \int_0^{\sqrt{\alpha} t} \exp( - r^2 ) {\bf d} r \\ & = & \frac{\sqrt{\pi}}{ 2 \sqrt{\alpha} {\rm erf} ( \sqrt{\alpha} t ) \end{array}$ $$If latex X_0 (t) < 0$$, we return \$code nan$$in order to inform code OdeErrControl$$ that its is taking to large a step. */ # include           // CppAD::Rosen34 # include namespace { // -------------------------------------------------------------- class Fun_three { private: const double alpha_; bool was_negative_; public: // constructor Fun_three(double alpha) : alpha_(alpha), was_negative_(false) { } // set f = x'(t) void Ode( const double                &t, const CppAD::vector &x, CppAD::vector       &f) {       f[0] = 2. * alpha_ * t * x[0]; f[1] = 1. / x[0]; // case where ODE does not make sense if( x[0] < 0. || x[1] < 0. ) {       was_negative_ = true; f[0] = CppAD::nan(0.); } } // set f_t = df / dt void Ode_ind( const double                &t, const CppAD::vector &x, CppAD::vector       &f_t) { f_t[0] =  2. * alpha_ * x[0]; f_t[1] = 0.; if( x[0] < 0. || x[1] < 0. ) {       was_negative_ = true; f_t[0] = CppAD::nan(0.); } } // set f_x = df / dx void Ode_dep( const double                &t, const CppAD::vector &x, CppAD::vector       &f_x) {       double x0_sq = x[0] * x[0]; f_x[0 * 2 + 0] = 2. * alpha_ * t;   // f0 w.r.t. x0 f_x[0 * 2 + 1] = 0.;                // f0 w.r.t. x1 f_x[1 * 2 + 0] = -1./x0_sq;         // f1 w.r.t. x0 f_x[1 * 2 + 1] = 0.;                // f1 w.r.t. x1 if( x[0] < 0. || x[1] < 0. ) {       was_negative_ = true; f_x[0] = CppAD::nan(0.); } } bool was_negative(void) {       return was_negative_; } }; // -------------------------------------------------------------- class Method_three { public: Fun_three F; // constructor Method_three(double alpha) : F(alpha) { } void step( double ta, double tb, CppAD::vector &xa , CppAD::vector &xb , CppAD::vector &eb ) {       xb = CppAD::Rosen34(F, 1, ta, tb, xa, eb); } size_t order(void) {       return 3; } }; } bool OdeErrControl_three(void) {       bool ok = true;     // initial return value double alpha = 10.; Method_three method(alpha); CppAD::vector xi(2); xi[0] = 1.; xi[1] = 0.; CppAD::vector eabs(2); eabs[0] = 1e-4; eabs[1] = 1e-4; // inputs double ti   = 0.; double tf   = 1.; double smin = 1e-4; double smax = 1.; double scur = 1.; double erel = 0.; // outputs CppAD::vector ef(2); CppAD::vector xf(2); CppAD::vector maxabs(2); size_t nstep; xf = OdeErrControl(method, ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep); double x0       = exp( alpha * tf * tf ); ok &= CppAD::NearEqual(x0, xf[0], 1e-4, 1e-4); ok &= CppAD::NearEqual(0., ef[0], 1e-4, 1e-4); double root_pi    = sqrt( 4. * atan(1.)); double root_alpha = sqrt( alpha ); double x1 = CppAD::erf(alpha * tf) * root_pi / (2 * root_alpha); ok &= CppAD::NearEqual(x1, xf[1], 1e-4, 1e-4); ok &= CppAD::NearEqual(0., ef[1], 1e-4, 1e-4); ok &= method.F.was_negative(); return ok; } // ========================================================================== bool OdeErrControl(void) {       bool ok = true; ok     &= OdeErrControl_one(); ok     &= OdeErrControl_two(); ok     &= OdeErrControl_three(); return ok; }