ACCELERATION OF RUNGE-KUTTA INTEGRATION SCHEMES

A simple accelerated third-order Runge-Kutta-type, fixed time 
step, integration scheme that uses just two function evaluations 
per step is developed. Because of the lower number of function 
evaluations, the scheme proposed herein has a lower computational 
cost than the standard third-order Runge-Kutta scheme while 
maintaining the same order of local accuracy. Numerical examples 
illustrating the computational efficiency and accuracy are 
presented and the actual speedup when the accelerated algorithm 
is implemented is also provided.


Introduction
The behavior of many systems in nature and society is described by sets of differential equations that are nonlinear in character.Obtaining closed-form solutions for these equations is therefore, generally speaking, impossible, especially for systems where the state vector has a high dimension.Thus numerical methods for solving such equations become very important.
One of the most common methods for numerically integrating sets of nonlinear differential equations is the Runge-Kutta (RK) family of single-step methods.The RK integrator of order n has a local error over a time step h of O(h n+1 ).For the third-and fourth-order (RK3 and RK4, respectively) methods, the number of function evaluations required over a single time step are three and four, respectively.In this paper, we propose a simple scheme that uses an RK-type third-order procedure [1,2] but which utilizes just two function evaluations at each step.The scheme has the same order of local accuracy as the standard RK3 method and hence appears to provide improvement in computational efficiency without any significant effect on computational accuracy.The drawback of the scheme, however, is that it is, in reality, a two-step scheme and therefore needs to be initially "started" using an accurate one-step method, like a higher-order RK integrator.

Numerical scheme
Consider the ordinary differential equation where y ∈ R n .We will assume that the function f is Lipschitz and the solution of the equation starting with the initial conditions y(0) = y 0 is unique.In this sequel, we will consider the following RK-type scheme to integrate this equation.We denote (2. 2) The approximate solution of (2.1) at time t n+1 is then given by where h is the fixed time step.We will determine the four parameters a 1 , a −1 , b, and β in relation (2.3) by comparing relation (2.3) with the Taylor series of y(t n + h).Using (2.1), we have (2.4) The Taylor's series expansion of y(t n + h) is simply Using relations (2.1), (2.2), and (2.4), we get (2.6) Using (2.6), we obtain (2.7) Comparing (2.7) with (2.5), we require (2.8) so that our numerical scheme is O(h 4 ).Using b as a free parameter, we now obtain some sample formulae using the scheme proposed in (2.2) and (2.3): ) ) ) Each set of values given in (2.9) represents an O(h 4 ) RK-type scheme.We note that our integration scheme (2.2) and (2.3) requires just two function evaluations at each time step.To compute y n+1 , we need to calculate at time t = t n only the two quantities k 1 and k 2 ; the quantities k −1 and k −2 are already known at time t n since they were computed at the previous time t n−1 .A two-step computational algorithm is thus obtained.To start the integration process, one would need a one-step integrator like RK4.However, after the start-up, only two function evaluations are required cutting the computation cost when compared to the standard RK3 algorithm, which is also accurate to O(h 4 ), by about a third.
It should be noted that each of the specific accelerated RK3 schemes obtained by using the numerical values illustrated in (2.9) will have, in general, different truncation errors (which will also depend on the specific differential equations being integrated) and different numerical stability characteristics.The examples described here seem to indicate that the parameters given in (2.9a) and (2.9b) may be an adequate set to use.Further investigation of these aspects will be taken up in a later communication.

Numerical examples
We consider two numerical examples to illustrate the accuracy of the method proposed herein.Comparisons of the error with the standard RK3 (see [1]) method are presented.The two equations considered are (3.1) whose solutions are y = e 1−t and y = √ 2/ (1 + t 2 ), respectively.The differential equation (3.2) is nonautonomous.We represent it, in the usual way, as a coupled set of ordinary nonlinear differential equations as where the dots refer to differentiation with respect to the parameter τ, and the initial conditions are taken to be y(τ = 1) = z(τ = 1) = 1.In fact, in the results presented below, (3.1) and (3.2) are both simultaneously solved as a system of three first-order differential equations by augmenting (3.1) with (3.3).The results described below refer to the implementation of the standard RK3 algorithm [1] and the accelerated RK3 algorithm proposed in (2.2) and (2.3).All computations are done in the Matlab environment.
In each case, the accelerated RK3 scheme was started by using an RK 45 variable step size scheme (ODE45 in Matlab) with a relative error tolerance of 10 −13 to compute the solution at the end of the first time step of duration h.
Figures 3.1 and 3.2 show the relative errors in the solution y(t n ) over a 10-second interval of time using a time step of 0.01.Our accelerated RK3 scheme is compared with the standard RK3 scheme provided in [1].Both schemes give the same order of magnitude of the relative local error.Figure 3.1 corresponds to the solution of the first example   (3.1).The pluses indicate results for the accelerated RK3 scheme and the circles indicate those for the standard RK3.The linear least-squares fits to the computed data points are shown by the dashed lines.The slopes of these dashed lines are 3.004 (for the pluses) and 3.005 (for the circles) in close conformity with the expected accuracy of O(h 4 ) per

Figure 3 . 1 .
Figure 3.1.The solid line shows the accelerated RK3 scheme given in (2.9b) for (3.1).The dashed line is for the standard RK3 algorithm that uses three function evaluations.The time step for the integration is h = 0.01.

Figure 3 . 2 .
Figure 3.2.The solid line shows the accelerated RK3 scheme given in (2.9b) for (3.2).The dashed line is for the standard RK3 algorithm that uses three function evaluations per step.
Figure 3.3.Absolute value of error at t = 10 in solution of (3.1) versus step size.The pluses show results for the accelerated RK3 scheme using the parameters given in (2.9b) and the circles show the standard RK3.