1 # ifndef CPPAD_UTILITY_ODE_GEAR_CONTROL_HPP
2 # define CPPAD_UTILITY_ODE_GEAR_CONTROL_HPP
373 template <
class Scalar,
class Vector,
class Fun>
390 CheckSimpleVector<Scalar, Vector>();
393 size_t n = size_t(xi.size());
397 "Error in OdeGearControl: M is less than one"
401 "Error in OdeGearControl: smin is greater than smax"
405 "Error in OdeGearControl: sini is greater than smax"
408 size_t(eabs.size()) == n,
409 "Error in OdeGearControl: size of eabs is not equal to n"
412 size_t(maxabs.size()) == n,
413 "Error in OdeGearControl: size of maxabs is not equal to n"
417 const Scalar zero(0);
419 const Scalar one_plus( Scalar(3) / Scalar(2) );
421 const Scalar ten(10);
427 Scalar step, sprevious, lambda, axi, a, root, r;
431 Vector X( (M + 1) * n );
443 for(i = 0; i < n; i++)
444 for(i = 0; i < n; i++)
450 else maxabs[i] = - xi[i];
469 else if( step <= smin )
473 if( tf <= T[m-1] + one_plus * step )
475 else T[m] = T[m-1] + step;
480 step = T[m] - T[m-1];
483 lambda = Scalar(10) * sprevious / step;
484 for(i = 0; i < n; i++)
485 { axi = X[m * n + i];
488 a = eabs[i] + erel * axi;
491 root = (a / e[i]) / ten;
493 { r = ( a / e[i] ) * step / (tf - ti);
494 root =
exp(
log(r) / Scalar(m-1) );
503 advance = one <= lambda || step <= one_plus * smin;
504 else advance = one <= lambda || step <= one_plus * sini;
512 for(k = 0; k < m; k++)
514 for(i = 0; i < n; i++)
515 X[k*n + i] = X[(k+1)*n + i];
519 for(i = 0; i < n; i++)
520 { ef[i] = ef[i] + e[i];
524 if( axi > maxabs[i] )
532 step = std::min(lambda , ten) * step / two;
534 for(i = 0; i < n; i++)
535 xf[i] = X[(m-1) * n + i];
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
void OdeGear(Fun &F, size_t m, size_t n, const Vector &T, Vector &X, Vector &e)
AD< Base > log(const AD< Base > &x)
AD< Base > exp(const AD< Base > &x)
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
Vector OdeGearControl(Fun &F, size_t M, const Scalar &ti, const Scalar &tf, const Vector &xi, const Scalar &smin, const Scalar &smax, Scalar &sini, const Vector &eabs, const Scalar &erel, Vector &ef, Vector &maxabs, size_t &nstep)