1 # ifndef CPPAD_UTILITY_ODE_ERR_CONTROL_HPP
2 # define CPPAD_UTILITY_ODE_ERR_CONTROL_HPP
422 template <
typename Scalar,
typename Vector,
typename Method>
438 CheckSimpleVector<Scalar, Vector>();
440 size_t n = size_t(xi.size());
444 "Error in OdeErrControl: smin > smax"
447 size_t(eabs.size()) == n,
448 "Error in OdeErrControl: size of eabs is not equal to n"
451 size_t(maxabs.size()) == n,
452 "Error in OdeErrControl: size of maxabs is not equal to n"
454 size_t m = method.order();
457 "Error in OdeErrControl: m is less than or equal one"
463 Vector xa(n), xb(n), eb(n), nan_vec(n);
470 Scalar m1(
double(m-1));
472 for(i = 0; i < n; i++)
473 { nan_vec[i] =
nan(zero);
478 else maxabs[i] = - xi[i];
483 Scalar tb, step, lambda, axbi, a, r, root;
484 while( ! (ta == tf) )
493 minimum_step = step <= smin;
498 if( tf <= ta + step * three / two )
504 method.step(ta, tb, xa, xb, eb);
509 if( (! ok) && minimum_step )
515 lambda = Scalar(10) * scur / step;
516 for(i = 0; i < n; i++)
517 {
if( zero <= xb[i] )
520 a = eabs[i] + erel * axbi;
521 if( ! (eb[i] == zero) )
522 { r = ( a / eb[i] ) * step / (tf - ti);
523 root =
exp(
log(r) / m1 );
528 if( ok && ( one <= lambda || step <= smin * three / two) )
532 for(i = 0; i < n; i++)
534 ef[i] = ef[i] + eb[i];
538 if( axbi > maxabs[i] )
546 else if( ! (ta == tf) )
549 scur = lambda * step / two;
555 template <
typename Scalar,
typename Vector,
typename Method>
567 { Vector maxabs(xi.size());
570 method, ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep
574 template <
typename Scalar,
typename Vector,
typename Method>
589 method, ti, tf, xi, smin, smax, scur, eabs, erel, ef, maxabs, nstep
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
AD< Base > log(const AD< Base > &x)
Scalar nan(const Scalar &zero)
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
AD< Base > exp(const AD< Base > &x)
bool hasnan(const Vector &v)
Vector OdeErrControl(Method &method, const Scalar &ti, const Scalar &tf, const Vector &xi, const Scalar &smin, const Scalar &smax, Scalar &scur, const Vector &eabs, const Scalar &erel, Vector &ef, Vector &maxabs, size_t &nstep)