Ignore:
Timestamp:
Sep 6, 2007 9:36:08 AM (13 years ago)
Author:
bradbell
Message:

trunk: Fix infinite loop in OdeErrControl?.

whats_new_07.omh: user's view of the changes.
ode_err_control.cpp: add test case for the infinite loop.
config.h: update version for subversion distribution option.
ode_err_control.hpp: fix infinite loop here.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/test_more/ode_err_control.cpp

    r985 r990  
    370370}
    371371
     372namespace {
     373        // --------------------------------------------------------------
     374        class Fun_four {
     375        private:
     376                 size_t n;   // dimension of the state space
     377        public:
     378                // constructor
     379                Fun_four(size_t n_) : n(n_)
     380                { }
     381
     382                // given x(0) = 0
     383                // solution is x_i (t) = t^(i+1)
     384                void Ode(
     385                        const double                &t,
     386                        const CppAD::vector<double> &x,
     387                        CppAD::vector<double>       &f)
     388                {       size_t i;
     389                        f[0] = CppAD::nan(0.);
     390                        for(i = 1; i < n; i++)
     391                                f[i] = (i+1) * x[i-1];
     392                }
     393        };
     394
     395        // --------------------------------------------------------------
     396        class Method_four {
     397        private:
     398                Fun_four F;
     399        public:
     400                // constructor
     401                Method_four(size_t n_) : F(n_)
     402                { }
     403                void step(
     404                        double ta,
     405                        double tb,
     406                        CppAD::vector<double> &xa ,
     407                        CppAD::vector<double> &xb ,
     408                        CppAD::vector<double> &eb )
     409                {       xb = CppAD::Runge45(F, 1, ta, tb, xa, eb);
     410                }
     411                size_t order(void)
     412                {       return 4; }
     413        };
     414}
     415
     416bool OdeErrControl_four(void)
     417{       bool   ok = true;     // initial return value
     418       
     419        // construct method for n component solution
     420        size_t  n = 6;       
     421        Method_four method(n);
     422
     423        // inputs to OdeErrControl
     424
     425        // special case where scur is converted to ti - tf
     426        // (so it is not equal to smin)
     427        double ti   = 0.;
     428        double tf   = .9;
     429        double smin = .8;
     430        double smax = 1.;
     431        double scur = smin;
     432        double erel = 1e-7;
     433
     434        CppAD::vector<double> xi(n);
     435        CppAD::vector<double> eabs(n);
     436        size_t i;
     437        for(i = 0; i < n; i++)
     438        {       xi[i]   = 0.;
     439                eabs[i] = 0.;
     440        }
     441
     442        // outputs from OdeErrControl
     443        CppAD::vector<double> ef(n);
     444        CppAD::vector<double> xf(n);
     445       
     446        xf = OdeErrControl(method,
     447                ti, tf, xi, smin, smax, scur, eabs, erel, ef);
     448
     449        // check that Fun_four always returning nan results in nan
     450        for(i = 0; i < n; i++)
     451        {       ok &= CppAD::isnan(xf[i]);
     452                ok &= CppAD::isnan(ef[i]);
     453        }
     454
     455        return ok;
     456}
     457
    372458// ==========================================================================
    373459bool OdeErrControl(void)
     
    376462        ok     &= OdeErrControl_two();
    377463        ok     &= OdeErrControl_three();
     464        ok     &= OdeErrControl_four();
    378465        return ok;
    379466}
Note: See TracChangeset for help on using the changeset viewer.