An adaptive theta method for the solution of stiff and nonstiff differential equations (Q1184718)

The paper deals with numerical methods for the solution of systems of ordinary differential equations (1) \(\dot y=f(x,y)\), \(a\leq x\leq b\), \(y(a)\) given. The theta method is widely used to solve such initial value problems. Suppose \(y_ n\) is an approximation to \(y(x_ n)\) and \(\dot y_ n\) to \(\dot y(x_ n)\), then the numerical solution at \(x_{n+1}\) where \(x_{n+1}=x_ n+h\) is denoted by \(y_{n+1}\) and is defined by \((2) y_{n+1}=y_ n+(1+\theta)h\dot y_ n+\theta hf(x_{n+1},y_{n+1})\). The theta method is \(S\)-stable for values of \(\theta > 0.5\) but is only conditionally stable for values of \(\theta < 0.5\). The formula is implicit; in order to calculate \(y_{n+1}\) is is necessary to choose values for \(h\) and \(\theta\), to select an iteration method and to have a predictor for the iteration method. In order to take reasonably sized steps when solving stiff problems the convergence criterion requires that a Newton-type iteration be used. The Jacobian evaluations required for the Newton iteration are expensive and often unnecessary for the solution of nonstiff problems where the stepsize restrictions are less severe and a cheaper iteration method such as functional iteration can be used instead. An obvious requirement is a single algorithm that can handle both stiff and nonstiff problems by automatically selecting the best value of theta and the most efficient iteration method with respect to step size and cost per step. There have been a number of attempts to find the best value of theta so that the integration is as accurate and efficient as possible. The authors derive an adaptive theta method integrator by directly choosing theta so as to try to minimize the number of time steps and combine this with the algorithm for switching between iteration methods. Numerical results are presented and these results illustrate the efficiency of the strategy both on standard test problems and on problems arising from the method of lines treatment of partial differential equations.