Application of second derivative RungeKutta collocation methods to stiff systems of initial value problems
Samaila MARKUS^{1} and Dauda Gulibur YAKUBU^{2*}
^{1}Department of Mathematics and Statistics, University of Maiduguri, PMB 1069, Maiduguri, Borno State, Nigeria
^{2}Department of Mathematical Sciences, Abubakar Tafawa Balewa University, PMB 0248, Bauchi, Bauchi State , Nigeria
Emails: daudagyakubu@gmail.com &smvalleh@gmail.com
*Corresponding author: Phone: 2348134168950 & 2348052168433
Abstract
In this paper we consider the construction of implicit second derivative RungeKutta collocation methods designed for the continuous numerical solution of stiff systems of first order initial value problems in ordinary differential equations. These methods are obtained based on the multistep collocation technique, which are shown to be convergent, with improved regions of absolute stability. Although the implementation of the second derivative RungeKutta collocation methods remains iterative due to the implicit nature of the methods, the advantage gained makes them suitable for solving stiff systems with eigenvalues of large modulus lying close to the imaginary axis. Some absolute stability characteristics and order of accuracy of the methods are studied. Finally, we show two possible ways of implementing the methods and compare them on some numerical examples found in the literature to demonstrate the high order of accuracy and reliability of the methods.
Keywords
Block hybrid scheme; Continuous scheme; Multistep collocation formula; Secondderivative RungeKutta method; Stiff system
Introduction
In the past few decades many authors proposed different numerical integration methods to obtain more accurate approximate solution to ordinary differential equations (ODEs), especially stiff systems of initial value problems in ordinary differential equations. However, extensive numerical experiments have shown that implicit integrators have the numerical superiority over the explicit methods when applied to solve stiff systems of ODEs [14]. In particular, implicit RungeKutta methods play an important role in the numerical solution of stiff systems of ODEs, because they have good stability properties and high order of convergence. In this paper we derive implicit secondderivative RungeKutta (SDRK) collocation methods designed for the continuous numerical solution of stiff systems of initial value problems in ODEs of the form

(1.1) 
Here the unknown function y is a mapping y:[x_{o},T]®R^{d} and the righthand side function f is f:[x_{0},T]_{x}R^{d}®R^{d} which is assumed to be sufficiently smooth, y_{0}®R^{d} is the given initial value and dy/dx is the differential operator. Let h>0 be a constant step size and we define the grid by x_{n}=x_{0}+nh, n=0,1,2,..,N where Nh=Tx_{0}and a set of uniformly spaced points on the integration interval is defined by x_{0}<x_{1}<x_{2}<…..<x_{n+1}=T The main reason for considering the second derivative terms is to derive a set of methods which are suitable for the continuous numerical solution of stiff differential systems with Jacobians having large eigenvalues lying close to the imaginary axis. However, for some important classes of problems it is necessary, for the sake of efficiency, to allow secondderivative integration methods to be used (see [58]) and in this case existing numerical methods tend to be much less satisfactory. Further, we will examine in detail the problem of implementing the secondderivative RungeKutta integration methods with fixed time steps. We show that, even though enormous gains in efficiency can be achieved if the methods are implemented in an appropriate way, there are still some important practical problems to be overcome, for example, the calculation of the second derivative terms in the methods which costs little higher than the first derivative terms [9].
There are several eminent authors who derived methods that can handle nonstiff, stiff, periodic and oscillatory problems handy [24,9], but little seems to have been done in deriving implicit secondderivative RungeKutta collocation methods with minimal function evaluations, uniformly accurate orders and maximal gain in efficiency. Hence, these motivate us to consider the construction of the implicit second derivative RungeKutta collocation methods in this paper.
Definition 1.1: A numerical method is said to be A(α)stable αÎ(0,π/2) if its region of absolute stability contains the infinite wedge W_{α} where

see [10]. 
Theorem 1.1: If f satisfies Lipschitz condition with constant L then the initial value problem y’(x)=f(x,y(x)), y(x_{0})=y_{0} possesses a unique solution on the interval [x_{0}, T] see [11].
Definition 1.2: A solution y(x) of (1.1) is said to be stable if given any Î>0 there is δ>0such that any other solution y(x) of (1.1) which satisfies

(1.2a) 
also satisfies

(1.2b) 
for all x>a.
The solution y(x) is asymptotically stable if in addition to (1.2b) as x®∞.
Material and method
The implicit second derivative RungeKutta collocation methods
In this section our objective is to describe the construction of the implicit secondderivative RungeKutta collocation methods based on the multistep collocation technique. In this regard we seek an approximate solution to the exact solution of (1.1) by the interpolant of the form

(1.3) 
which is twice continuously differentiable. We set the sum r+s+t to be equal to p so as to be able to determine {α_{i}} in (1.3). In this formulation r denotes the number of interpolation points used and s>0, t>0are distinct collocation points. Interpolating (1.3) at the points {x_{n+j}}, and collocating y’(x) and y”(x) at the points {c_{n+j}} we obtain the following system of equations

(1.4) 

(1.5) 

(1.6) 
In fact equations (1.4) to (1.6) can be expressed in the matrixvector form as

(1.7) 
where the psquare matrix D, p vectors α and y are defined as follows:

(1.8) 


where D’=(p1) and D’’=(p1)(p2) in (1.8) represent the first and second derivatives respectively and correspond to the differentiation with respect to x. Similar to the Vandermonde matrix, the matrix D in (1.7) is nonsingular. Consequently, equation (1.7) has the unique solution given by
α=C_{y} where C=D1 (1.9)
Rearranging equations (1.7) to (1.9) we obtain the multistep collocation formula of the type in [12] which was a generalization of [13] and here we extend to second derivative as follows.

(1.10) 
where: 







Here α_{j}(x), β_{j}(x), ω_{j}(x) are continuous coefficients of the method which are to be determined. They are assumed polynomials of the form

(1.11) 



The numerical constant coefficients α_{j,i+1}, β_{j,i+1}, ω_{j,i+1} in (1.11) are to be determined.
In fact, the above coefficients can be obtained from the components of the matrix D^{1}.
That is, if the identity (1.12) holds. Actual evaluations of the matrices C and D are carried out with a computer algebra system, for example, Maple, to determine the constant coefficients α_{j,i+1}, β_{j,i+1}, and ω_{j,i+1} in (1.11).

(1.12) 
To obtain the continuous scheme, we insert (1.11) into (1.10) to have

(1.13) 




where: 



From (1.13) we have

(1.14) 



Recall that p=r+s+t, such that (1.14) reduces to

(1.15) 
Thus, expanding (1.15) fully we obtain the proposed continuous scheme as follows:
_{} (1.16)
where T denotes transpose of the matrix C in (1.12) and the vector (1,x,x^{2},…,x^{r+s+t1}).
Remark 1.
We call D the multistep collocation and interpolation matrix which has a very simple structure. From (1.8), the columns of D which give the continuous coefficients α_{j}(x), β_{j}(x), ω_{j}(x) can be obtain from the corresponding columns of C. As can be seen the entries of C are the constant coefficients of the polynomial given in (1.11) which are to be determined. The matrix C is the solution vector (output) and D is termed the data (input), which is assumed to be nonsingular for the existence of the inverse matrix C.
In the secondderivative methods, we see that not only the function f(x,y) is evaluated at some intermediate points, but in addition the functions Df, D^{2}f where Dis the differential operator (see [14]). Hence, in addition to the computation of the fvalues at the internal stages in the standard RungeKutta methods, the modified methods involve computing gvalues, where f and g are as defined in (1.10).
According to [15] these methods can be practical if the costs of evaluating g are comparable to those in evaluating f and can even be more efficient than the standard RungeKutta methods if the number of function evaluations is fewer. It is convenient to rewrite the coefficients of the defining formula (1.10) evaluated at some certain points in the block matrix form as,

(1.17a) 

(1.17b) 
where A=[a_{ij}]_{sxs}, Â=[â_{ij}]_{sxs} indicate the dependence of the stages on the derivatives found on the other stages and b=[b_{i}]_{sx1}, _{}are vectors of quadrature weights showing how the final result depends on the derivatives computed at the various stages, I is the identity matrix of size equal to the differential equation system to be solved and N is the dimension of the system. Also Ä is the Kronecker product of two matrices. For simplicity, we write the method in (1.17) as follows:

(1.18a) 

(1.18b) 
and the block vectors in R^{sN} are defined by



(1.19) 
The coefficients of the implicit secondderivative RungeKutta collocation methods can be conveniently represented more compactly in an extended partitioned Butcher Tableau, of the form

(1.20) 
where c=[1]_{sx1} is the abscissae vectors which indicates the position within the step of the stage values.
Specification of the implicit second derivative RungeKutta collocation methods
A fourthorder implicit second derivative RungeKutta collocation method
In this section we develop the general form of conditions for the coefficients of the implicit secondderivative RungeKutta collocation methods. To obtain the coefficients of the first method, the matrix D in (1.8) takes the form,
_{} 
(1.21) 
Inverting the matrix D in (1.21) once, we obtain the continuous scheme of the form in (1.16) as follows
_{} 
(1.22) 
where
_{} 
_{} 
_{} 
_{} 
_{} 
To get the value of u in (1.21), we find the zero of the N^{th} degree Jacobi polynomial (see, [16]) defined by

(1.23a) 
where

(1.23b) 
with γ_{0}=1.
Next by substituting u = ⅓ into the continuous scheme (1.22) and evaluate at x=x_{n+1} and x_{n+u} we obtain the following block hybrid discrete scheme, which can be applied simultaneously as block method for dense output, if desired,

(1.24a) 

(1.24b) 
If we convert the block hybrid discrete scheme (1.24(ab)) to implicit secondderivative RungeKutta collocation method and write the method in the form of (1.18) we have the following

(1.25) 
The internal stage values at the n^{th} step are calculated as,



with the stage derivatives as follows: 



A sixthorder implicit second derivative RungeKutta collocation method
Here we derive implicit secondderivative RungeKutta collocation method of higher order for the numerical solution of stiff systems. The matrix D in (1.8) takes the following form
_{} 
(1.26) 
where u and v are the zeros of p2(x)=0. Jacobi polynomial in (1.23) of degree which are valid in the interval [x_{n},x_{n+1}]. Inverting the matrix in (1.26) once we obtain the continuous scheme as follows:
_{} 
(1.27) 
where
_{}, 
_{} 
_{} 
_{} 
_{} 
_{} 
_{} 
Evaluating the continuous scheme (1.27) at the points {x_{n+1},x_{n+u},x_{n+v}} we obtain the block hybrid discrete scheme which can also be applied simultaneously as block method for dense output, if desired,
_{} 

_{} 
(1.28) 
_{} 

Writing the method in the formalism of (1.18) we have
_{} 
(1.29) 
The internal stage values at the n^{th} step are calculated as,
_{} 
_{} 
_{} 
_{} 
where the stage derivatives are as follows:
_{} 
_{} 
_{} 
_{} 
Analysis of the properties of the implicit SDRK collocation methods
Order, Consistency, ZeroStability and Convergence of the Implicit SDRKCMs
With the multistep collocation formula (1.10) we associate the linear difference operator _{} defined by:
_{} 
(1.30) 
where y(x) is an arbitrary function, continuously differentiable on [x_{0}, T]. Following [17], we can write the terms in (1.30) as a Taylor series expansion about the point x to obtain the expression,
_{} 
(1.31) 
where the constant coefficients C_{p}, p = 0,1,2,… are given as follows:
_{} 
_{} 
_{} 
_{} 
According to [17], the multistep collocation formula (1.10) has order p if
_{} 
(1.32) 
Therefore, C_{p+1 }is the error constant and C_{p+1}h^{p+1}y^{(p+1)}(x_{n}) is the principal local truncation error at the point x_{n} (see [18]). Hence from our calculation the order and error constants for the constructed methods are presented in Table 1. It is clear from the Table that the implicit secondderivative RungeKutta collocation methods are of highorder. They have smaller error constants and hence more accurate than the conventional RadauRungeKutta methods of the same order of convergence.
Table 1. Order and error constants of the implicit SDRK collocation methods
Method 
Order 
Error constant 
Method (1.25) 
P=4 
C_{5}=9.8765_{}10^{3} 
Method (1.29) 
P=6 
C_{7}=2.5714_{}10^{4} 
Definition 1.3: Consistency
A linear multistep method is said to be consistent if the order of the method is greater or equal to one, that is if p>1 (see, [17]).(i) ρ(1)=0_{ }and(ii)ρ’(1)=σ(1), where ρ(z) and σ(z) are respectively the 1^{st} and 2^{nd} characteristic polynomials.
From Table 1 and definition 1.3, we can attest that the implicit secondderivative RungeKutta collocation methods are consistent.
Definition 1.4: Zerostability
A linear multistep method is said to be zerostable if the roots:
_{} 


see [17].
Based on definition 1.4, the newly constructed implicit secondderivative RungeKutta collocation methods are zerostable.
Definition 1.5: Convergence
The necessary and sufficient conditions for a linear multistep method to be convergent are that it must be consistent and zerostable (see [10]). Hence, from definitions 1.3 and 1.4 the implicit secondderivative RungeKutta collocation methods are convergent.
Regions of absolute stability of the implicit SDRK collocation methods
To study the stability properties of the implicit secondderivative RungeKutta collocation methods we reformulate them as general linear methods (see [19]). Hence, we use the notations introduced by [20] which a general linear method is represented by the partitioned (s+r)×(s+r) matrix, (containing A, U, B, V),
_{} n = 1, 2, …,N, 
(1.33a) 
where
_{} 
_{} 
_{} 
_{} 
_{} 
_{} 
_{} 
_{} 
and e=[1,…,1]^{t}ÎR. Hence (1.33a) takes the form
_{} 
(1.33b) 
where r denotes quantities as output from each step and input to the next step and s denotes stage values used in the computation of the step y_{1},…,y_{s}. The coefficients of these matrices indicate the relationship between the various numerical quantities that arise in the computation of stability regions. The elements of the matrices A, U, B and V are substituted into the stability matrix. In the sense of [21] we apply (1.33) to the linear test equation y’=λy,x>0 and λÎC but for the secondderivative highorder methods we use y’’=λ2y which leads to the recurrence relation y^{[n+1]}=M(z)y^{[n]}, n=1,2,…,N1, z=λh, where the stability matrix M(z) is defined by
_{} 
(1.34) 
We also define the stability polynomial p(η,z) by the relation
_{} 
(1.35) 
and the absolute stability region _{}of the method is given by
_{} 
To compute the regions of absolute stability we substitute the elements of the matrices A, U, B and V into the stability function (1.34) and finally into the stability polynomial (1.35) of the methods, which is plotted to produce the required graphs of the absolute stability regions of the methods as shown in Figure 1.

Method (1.25) is A(_{})stable 
Method(1.29) is A(_{})stable 
Figure 1. Regions of absolute stability of the implicit SDRK collocation methods
The methods constructed in this paper are A(α)stable, since their regions consist of the complex plane outside the enclosed Figures, with only very small portion inside the left handhalf plane.
Numerical experiments
Our aim in this section is to compare the performance of the new methods with other methods of the conventional type (Radau methods) which are widely known to be among the most efficient methods for the numerical solution of stiff ordinary differential equations. Therefore, we report the numerical results obtained from the applications of the new implicit secondderivative RungeKutta collocation methods. We present the computed results in the formalism of [22]. In the presentation we use nfe to denote the number of function evaluations and Ext to indicate the exact solutions.
Example 1. Consider the stiff system of first order differential equation,
_{}
The exact solution is
_{}
This system has eigenvalues of large modulus lying close to the imaginary axis 1±15i. We solve the system in the range of [0, 50] with h = 0.1and the results obtained are shown in Table 2. In this example we compare the results from the two new methods side by side in the Table. The solution curves obtained are compared with the exact (Ext) solutions in Figure 2. Thus, for fair comparison see the results from the third derivative methods of [23] Table 4 page 447 and results of [24] Table 2 page 166.
Table 2. Absolute errors in the numerical integration of example 1
x 
Y_{i} 
Method (1.25) 
Method (1.29) 
5 
Y_{1} 
3.46368813864673 ×10^{3} 
2.16834805259769 ×10^{3} 
Y_{2} 
3.94364780750056 ×10^{4} 
6.23671456063679 ×10^{4} 

50 
Y_{1} 
1.49495218771792 ×10^{4} 
7.38442455972734 ×10^{5} 
Y_{2} 
4.86742095474412 ×10^{9} 
2.19934956374418 ×10^{9} 

150 
Y_{1} 
1.85517639473028 ×10^{4} 
1.20177414793520 ×10^{4} 
Y_{2} 
1.45786198883123 ×10^{8} 
7.02697714351372 ×10^{9} 

250 
Y_{1} 
1.25223273382179 ×10^{13} 
1.99479371074799 ×10^{14} 
Y_{2} 
8.66743589934986 ×10^{13} 
3.55424556071686 ×10^{13} 

500 
Y_{1} 
1.64891543846432 ×10^{23} 
2.32684011997371 ×10^{24} 
Y_{2} 
4.31480861603385 ×10^{24} 
3.70791674680754 ×10^{24} 
a) b)
Figure 2. Graphical plots of example 1 using Second Derivative RungeKutta Collocation methods: a) Solution curve of example 1 using method (1.25), with nfe =500; b) Solution curve of example 1 using method (1.29), with nfe =500
Example 2: The second example is a stiff linear system of equation, see [25],
_{}
The exact solution is given by
_{}
We solve this system in the range [0, 10]. In this example we also compare the results from the two new methods in Table 3 and the efficiency curves obtained are compared with the exact solutions in Figure 3.
Table 3. Absolute errors in the numerical integration of example 2
x 
Y_{i} 
Method (1.25) 
Method (1.29) 
5 
Y_{1} 
7.00794977603891 × 10^{12} 
2.22044604925031 × 10^{16} 
Y_{2} 
7.04547531427124 × 10^{12} 
0 

50 
Y_{1} 
2.27042828981894 × 10^{11} 
1.11022302462516 × 10^{15} 
Y_{2} 
2.27533547558778 × 10^{11} 
1.11022302462516 × 10^{15} 

150 
Y_{1} 
7.35667637918880 × 10^{11} 
4.44089209850063 × 10^{16} 
Y_{2} 
7.35377314597940 × 10^{11} 
3.33066907387547 × 10^{16} 

250 
Y_{1} 
4.08534317486442 × 10^{11} 
1.11022302462516 × 10^{16} 
Y_{2} 
4.08287292863463 × 10^{11} 
1.11022302462516 × 10^{16} 

500 
Y_{1} 
6.84553524976650 × 10^{11} 
5.55111512312578 × 10^{16} 
Y_{2} 
6.84068357514889 × 10^{11} 
5.55111512312578 × 10^{16} 
a) b)
Figure 3. Graphical plots of example 2 using second derivative RungeKutta collocation methods: a) Solution curve of example 2 using method (1.25), with nfe =500;b) Solution curve of example 2 using method (1.29), with nfe =500
Example 3: The third test example is a stiff problem taken from [25].
The exact solution is,
_{}
The eigenvalues of the system are λ_{1}=50 and λ_{2,3} = 0.1±8i. The computed solutions of the example within the interval of [0, 1] are presented in Table 4 and the graphical plots are shown in Figure 4.
Example 4: We consider another linear problem which is particularly referred to by [27] as a troublesome problem for some existing methods because the eigenvalues of this problem are lying close to the imaginary axis where some stiff integrators were known to be inefficient. The eigenvalues of the Jacobian are λ_{1,2} = 10±100i, λ_{3} = 4, λ_{4} = 1, λ_{5 } = 0.5 and λ_{6} = 0.1.
_{}
Table 4. Absolute errors in the numerical integration of example 3
x 
Y_{i} 
Radau IIA method [20, 26] 
Method (1.29) 
5 
Y_{1} 
3.70353636647280 × 10^{10} 
3.12638803734444 × 10^{13} 
Y_{2} 
3.70363018031838 × 10^{10} 
3.12583292583213 × 10^{13} 

Y_{3} 
3.70344199751571 × 10^{10} 
3.12638803734444 × 10^{13} 

50 
Y_{1} 
4.98715513330694 × 10^{11} 
4.27435864480685 × 10^{14} 
Y_{2} 
5.00358643407139 × 10^{11} 
4.27435864480685 × 10^{14} 

Y_{3} 
4.97826224687969 × 10^{11} 
4.21884749357559 × 10^{14} 

150 
Y_{1} 
2.59459120854899 × 10^{13} 
8.88178419700125 × 10^{16} 
Y_{2} 
2.37143638059933 × 10^{13} 
5.99520433297585 × 10^{15} 

Y_{3} 
1.55916946020795 × 10^{14} 
6.77236045021345 × 10^{15} 

250 
Y_{1} 
4.30766533554561 × 10^{13} 
7.10542735760101 × 10^{15} 
Y_{2} 
4.39759340054025 × 10^{13} 
2.88657986402541 × 10^{15} 

Y3 
8.70636895911048 × 10^{13} 
4.44089209850063 × 10^{15} 

500 
Y_{1} 
1.29318777908338 × 10^{12} 
5.99520433297585 × 10^{15} 
Y_{2} 
9.13435993510348 × 10^{14} 
8.93729534823251 × 10^{15} 

Y3 
1.20203846876166 × 10^{12} 
3.55271367880050 × 10^{15} 
a) b)
Figure 4. Graphical plots of example 3 using Radau IIA method [20,26] and SDRKC method: a) Solution curve of example 3 using Radau IIA method [20,26], with nfe =500, b) Solution curve of example 3 using method (1.29), with nfe =500
We solve this problem in the range 0 ≤ x ≤ 1. The computed solutions are shown in Table 5, while the graphical plots are displayed in Figure 5. Though, only the first four components {y_{1},y_{2},y_{3},y_{4}} are shown in the Table of values.
Example 5: HIRES
The High irradiance responses problem consists of a stiff system of 8 nonlinear ordinary differential equation which originates from plant physiology and describes how light is involved in morphogenesis. It explains the high irradiance responses of photo morphogenesis on the basis of phytochrome by means of a chemical reaction involving eight reactants in the form of initial value problem see [28,29].
Table 5. Absolute errors in the numerical integration of example 5
x 
Y_{i} 
Radau IIA method [20, 26] 
Method (1.29) 
5 
Y_{1} 
3.83760043742853 × 10^{8} 
5.84532422465145 × 10^{11} 
Y_{2} 
2.90659031079721 × 10^{8} 
5.86538248525947 × 10^{11} 

Y_{3} 
2.22044604925031 × 10^{16} 
0 

Y_{4} 
0 
1.11022302462516 × 10^{16} 

50 
Y_{1} 
2.32388491339108 × 10^{7} 
1.51528067870998 × 10^{10} 
Y_{2} 
5.72332533588238 × 10^{8} 
3.82774117957396 × 10^{10} 

Y_{3} 
1.55431223447522 × 10^{15} 
3.33066907387547 × 10^{16} 

Y_{4} 
1.11022302462516 × 10^{16} 
4.44089209850063 × 10^{16} 

250 
Y_{1} 
1.91869619220464 × 10^{8} 
1.42258618079927 × 10^{11} 
Y_{2} 
1.09619847643572 × 10^{8} 
3.52476554887904 × 10^{11} 

Y_{3} 
1.38777878078145 × 10^{15} 
1.66533453693773 × 10^{16} 

Y_{4} 
7.77156117237610 × 10^{16} 
2.22044604925031 × 10^{16} 

500 
Y_{1} 
2.28890865933633 × 10^{10} 
2.65111307337423 × 10^{13} 
Y_{2} 
1.86748006761241 × 10^{10} 
4.33490943878637 × 10^{13} 

Y3 
3.57353036051222 × 10^{16} 
6.93889390390723 × 10^{17} 

Y_{4} 
4.44089209850063 × 10^{16} 

a) b)
Figure 5. Graphical plots of example 4 using Radau IIA method [20,26] and SDRKC method: a) Solution curve of example 4 using Radau IIA method [20,26], with nfe =500, b) Solution curve of example 4 using method (1.29), with nfe =500
_{}
We solve the problem in the interval [0, 1] and compare the graphical plots with the ODE code of Matlab in Figure 6.
a) b)
Figure 6. Graphical plots of example 5 using SDRKCM and the ODE code of Matlab: a) Solution curve of example 5 using method (1.25), with nfe = 500, b) Solution curve of example 5 using method (1.29), with nfe = 500
Concluding remarks
In conclusion, the methods so derived in this paper so far have been to introduce a class of implicit secondderivative RungeKutta collocation methods suitable for the approximate numerical integration of systems, particularly stiff systems of ordinary differential equations. The derived methods proved an efficient way to find numerical solution to systems of initial value problems when the second derivative terms are cheap to evaluate. We present two new methods of orders four and six that are intended for accurate integration of stiff systems of equations. We have also compared the numerical solution obtained with the conventional Radau RungeKutta methods which are widely recognized to be among the most efficient methods for stiff systems [3, 20, 26]. The solution curves are also compared with the exact solutions graphically in Figures.
Acknowledgement
This work was supported by Tertiary Education Trust Fund (TETF) Ref. No. TETF /DAST &D.D/6.13/NOMCA/BAS&BNAS. Therefore the authors gratefully acknowledged the financial support of the TETF.
The authors are also grateful to the anonymous referee for his/her useful comments and suggestions
References
1. Butcher J. C., Coefficients for study of RungeKutta integration processes, J. Austral. Math. Society, 1963, 3(2), p. 185 201.
2. Enright W. H., Secondderivative multistep methods for stiff ordinary differential equations, SIAM J. Numer. Anal., 1974, 11(2), p. 321341.
3. Hairer E., Nørsett S. P., Wanner G., Solving Ordinary Differential Equations II: Stiff and Differential Algebraic Problems, Second Revised Edition, Berlin Germany, Springer, 1996.
4. VigoAguiar J., Ramos H., A family of Astable RungeKutta collocation methods of higher order for initial value problems, IMA J. Numer Anal., 2007, 27(4), p.798817.
5. Cash J. R., High order methods for the numerical integration of ordinary differential equations, Numer. Math., 1978, 30, p. 385409.
6. Mitsui T., RungeKutta type integration formulas including the evaluation of the second derivative, Part I, Publ. RIMS, Kyoto University, Japan, 1982, 18, p. 325364.
7. Mitsui T., Yakubu D. G., TwoStep Family of LookAhead linear multistep methods for ODEs, The science and Engineering Review Doshisha University, Japan, 2011, 52(3), p. 181188.
8. Urabe M., An implicit onestep method of highorder accuracy for the numerical integration of ordinary differential equations, J. Numer. Math., 1970, 15, p. 151164.
9. Yakubu D. G., Kwami A. M., Implicit twoderivative RungeKutta collocation methods for systems of initial value problems, J. Nigerian Math. Society, 2015, 34, p. 128142.
10. Dahlquist G., A special stability problem for linear multistep methods, BIT Numer. Math., 1963, 3, p. 2743.
11. Butcher J. C., Order and stability of RungeKutta methods, Nanjing University Conf., 2014.
12. Onumanyi P., Awoyemi D. O., Jator S. N., Sirisena U. W., New linear multistep methods with continuous coefficients for first order initial value problems in ODEs, J. Nigerian Math. Society, 1994, 13, p. 3751.
13. Lie I., Nørsett S. P., Superconvergence for multistep collocation, Math. Comp., 1989, 52, p. 6580.
14. Kastlunger K. H., Wanner G., RungeKutta processes with multiple nodes, J. Computing, 1972, 9(9), p. 924.
15. Chan R. P. K., Tsai A. Y. J., On explicit two derivative RungeKutta methods, Numer. Algor., 2010, 53, p. 171194.
16. Villadsen J., Michelsen M. L., Solution of Differential Equation Models by Polynomial Approximations, PrenticeHall, INC, Eaglewood Cliffs, New Jersey, 1978, p. 112.
17. Lambert J. D., Computational Methods in Ordinary Differential Equations, John Willy, and Sons, New York, 1973.
18. Mitsui T., A modified version of Urabe’s implicit singlestep method, J. Comp. Appl. Math., 1987, 20, p. 325332.
19. Burrage K., Butcher J. C., Stability criteria for implicit RungeKutta methods, SIAM J. Numer. Anal., 1979, 16, p. 4657.
20. Butcher J. C., Numerical Methods for Ordinary Differential Equations, 2^{nd} ed. John Wiley & Sons Ltd, 2008, p. 225.
21. Chollom J. P., Jackiewicz Z., Construction of two step RungeKutta (TSRK) methods with large regions of absolute stability, J. Comp. and Appl. Maths., 2003, 157, p. 125137.
22. Butcher J. C., Hojjati G., Second derivative methods with RungeKutta stability, Numer. Algor., 2005, 40, p. 415429.
23. Ali K. E., Hojjati G., Third derivative multistep methods for stiff systems, Intern. J. of Nonlinear Science, 2012, 14(4), p. 443450.
24. Akinfenwa O. A., Akinnukawe B., Mudasiru S. B., A family of continuous third derivative block methods for solving stiff systems of first order ordinary differential equations, J. Nigerian Math. Society, 2015, 34, p. 160168.
25. Lambert J. D., Numerical Methods for Ordinary Differential System, John Willy, & Sons, New York, 1991.
26. Butcher J. C., Integration processes based on Radau quadrature formulas, Math. Compt., 1964, 18, p. 233244.
27. Fatunla S. O., Numerical integrators for stiff and highly oscillatory differential equations, J. Math. Compt., 1980, 34(150), p. 373390.
28. GonzalezPinto S., GonzalezConcepcion C., Montijino J. I., Iterative schemes for Gauss methods, J. Comput. Math. Applic., 1994, 27(7), p. 67–81.
29. Schäfer E., Freiburi Br., A new approach to explain the “High Irradiance Responses” of photomorphogenesis on the basis of phytochrome, J. Math. Bio., 1975, 2, p. 4156.