An Accurate Block Solver for Stiff Initial Value Problems

New implicit block formulae that compute solution of stiff initial value problems at two points simultaneously are derived and implemented in a variable step size mode. The strategy for changing the step size for optimum performance involves halving, increasing by a multiple of 1.7, or maintaining the current step size. The stability analysis of the methods indicates their suitability for solving stiff problems. Numerical results are given and compared with some existing backward differentiation formula algorithms. The results indicate an improvement in terms of accuracy. 1. Introduction Stiff problems occur in many fields of engineering science and they represent coupled physical systems having components varying with very different time scales [1]. Some physical examples include chemical kinetics where reactions go at very different speeds and in circuit simulation where components respond at widely different time scales. A considerable research has been done on numerical methods for solving stiff initial value problems (IVPs) of the following form: The most successful ones are implicit in nature. Examples include single step methods like the diagonally implicit Runge Kutta (DIRK) methods, single diagonally implicit Runge-Kutta (SDIRK), and singly implicit Runge-Kutta methods [2–4]. Also well known are methods based on Richardson extrapolation and deferred correction schemes [5–8]. Others are multistep methods like the backward differentiation formula (BDF) (see, e.g., [4]). Initially, the BDF methods are known to be -stable only up to order 2 [9]. However, Cash in [10] proved that the barrier in [9] can be circumvented by adding extra future point to produce higher order -stable methods. Methods of this type are known as extended backward differentiation formula [10–12]. Some block methods for stiff problems are also known to be A-stable and have order greater than two (see, e.g., [13–18]). One of the motivation of this paper is to develop a new block method based on the block backward differentiation formula [13, 14, 17]. The method is similar in fashion with that in [15] but differs greatly by the choice of the constant . Furthermore, our choice of gives rise to different formulae with better stability regions than in [15], even with the choice of the same step ratio . In the following sections of the paper, we will present the derivation of the method, its stability analysis, and implementation. Standard test problems will also be presented and a comparison of their numerical results with some known stiff solvers will be given. 2. Formulation of