The Algorithm

The method used for integration is an embedded RK2(3) scheme due to Fehlberg, and reported on page 170 of Hairer et al. (1987). Lagged variables (and gradients) are stored in a ringbuffer at each step of the integrator. Interpolation is required to estimate values of the lagged variables between storage times. For numerical probity it is essential that the interpolation of lagged variables is of a higher order of approximation than the integrator, otherwise the assumptions underlying the error estimate from the RK pair will not be met. The algorithm used in Solv95 uses cubic hermite interpolation (e.g. Burden and Faires 1987) to achieve this (which is the reason that gradients need to be stored along with lagged values). The consequences of not using consistent interpolation and integration schemes are vividly illustrated in Highman (1993). Paul (1992) was also influential in the design of the method used here, and the step size selection is straight out of Press et al. (1992) (method, not code!). The RK2(3) pair used is not actually optimal - it should be possible to derive an improved scheme - see Butcher (1987) for an explanation of how to go about it.

A deficiency of the algorithm is that switches are not tracked automatically - if a switch impacts on a lagged variable the resulting derivative discontinuities are only dealt with by adaptive timestepping through them (Note, however, that lagged variables and gradients are stored immediately before and immediately after the application of a switch - i.e. the interpolator is not applied blindly through switches!).

  1. Burden, R.L. and J.D. Faires (1985) Numerical Analysis. Pridle Weber and Schmidt, Boston.
  2. Butcher, J.C. (1987) The Numerical Analysis of Ordinary Differential Equations. John Wiley & Sons, Chichester.
  3. Hairer, E., S.P.Norsett & G.Wanner (1987) Solving Ordinary differential Equations I. Springer-Verlag Berlin. p170 RKF2(3)B
  4. Highman, D.J. (1993) Appl. Numer. Math. 12:403-414
  5. Paul, C.A.H (1992) Appl. Numer. Math. 9:403-414
  6. Press et al. (1992) Numerical Recipes in C. CUP

Back to index