The problem of stiffness arises in the solution of ordinary differential equations (ODEs) when the solution varies very rapidly with time in some time intervals and much more slowly in others. Many ODE problems in science and engineering applications exhibit some degree of stiffness. For a stiff ODE problem, explicit solution methods, even with adaptive step size, may fail, or be unacceptably slow due to the extremely small step size that may be required in the intervals with rapid change. To address this issue generally requires use of implicit solution methods, often involving a backwards difference approach, using appropriate backwards difference formulas (BDFs) or numerical differentiation formulas (NDFs) [1]. Adaptive step sizes can also be used in this context, usually by comparing the implicit solution estimate with some form of higher-order truncated Taylor approximation to obtain an error estimate, but there still may be performance and/or accuracy issues since there is no guarantee that either the implicit solution estimate or the truncated approximation is near the true solution.
In this paper we introduce an adaptive step size approach based on a dynamic error estimate obtained from a pair of âsqueezeâ functionals that bound the true solution and its Taylor approximation (of desired order). By monitoring the size and growth rate of the bound gap, relative to specified tolerances, it can be determined whether the Taylor approximations remain sufficiently accurate, and then the length of the time step can be adapted accordingly. Results of preliminary tests of this approach on a variety of stiff ODE problems will be presented, indicating that it is more efficient and reliable than standard Matlab solvers (e.g., ode15s, ode23s, ode23t) for stiff problems. Test problems include a flame propagation problem [2], the Robertsonâs chemical reaction model [3], the van der Pol oscillator, a double pendulum problem and a three-body problem.
[1] Shampine, L. F. and M. W. Reichelt, âThe MATLAB ODE Suite,â SIAM Journal on Scientific Computing, 18, 1-22 (1998).
[2] Jannelli A., Fazio R. âAdaptive stiff solvers at low accuracy and complexityâ, J. of Computational and Applied Mathematics, 191, 246-258 (2006).
[3] Robertson, H.H. The solution of a set of reaction rate equations. In Numerical Analysis: An Introduction, J. Walsh, Ed., Academic Press: London, UK, 1966.