QUAM- A novel algorithm for the numerical integration of stiff ordinary differential equations

Download

Microsoft Word Download (203Kb)

Description/Abstract

This work proposes a novel algorithm for the numerical computation of the solution of ordinary differential equations, particularly for the case of ‘stiff’ equations. Existing methods such as implicit Euler and implicit Trapezoidal algorithms, and backward difference formulas, are effective but do not integrate the fast modes correctly and at each timestep must undertake a Newtonian search to get the next solution value. Here we present the QUAM algorithm which advances to the next time step Xr+1 by making a first order Taylor expansion of gradient function F(x,t) about the current value Xr, and uses exact analytical expressions to derive Xr+1. Two analytic approaches are possible. One either derives an analytic matrix expression requiring the inverse of the gradient matrix A, or one performs an eigendecomposition of A to get the same result. An adaptive procedure with relatively low overheads, that adjusts timestep to keep one step error within bounds, is integrated into the algorithm. The algorithm was first tested on a stiff linear problem, and error analysis confirmed that the method is basically second order accurate. Next the adaptive QUAM algorithm was tested on the classic Robertson problem, where its performance compared very favourably with the MATLAB routine ode23s. The algorithm is particularly fast when the gradient matrix is either slowly varying or even constant.