Proficiency of Second Derivative Schemes for the Numerical Solution of Stiff Systems

James Adewale^1*, Adesanya Olaide², Onsachi Oziohu², Sunday Joshua³, Moses Omuya¹

¹Mathematics Department, American University of Nigeria, Yola, Nigeria.
²Department of Mathematics, Modibbo Adama University of Technology, Yola, Nigeria.
³Department of Mathematics, Adamawa State University, Mubi, Nigeria.

DOI: 10.4236/ajcm.2018.81008                

Abstract

This paper presents a study on the development and implementation of a second derivative method for the solution of stiff first order initial value problems of ordinary differential equations using method of interpolation and collocation of polynomial approximate solution. The results of this paper bring some useful information. The constructed methods are A-stable up to order 8. As it is shown in the numerical examples, the new methods are superior for stiff systems.

Keywords

Second Derivative, Interpolation, Collocation, Continuous Scheme, Block Method, Stiff Problems, Initial Value, Linear Multistep Method

1. Introduction

We considered development of second derivative method for the solution of

$y^{\prime }=f\left(x,y\right),y\left(x_n\right)=y_0,x_n\le x\le x_N$ (1)

where $x_n$ is the initial points, $y\mathrm{<mo>\; :\; </mo>}\left[x_n\mathrm{,}{x}_N\right]\to R^m$ , $f\mathrm{<mo>\; :\; </mo>}\left[x_n\mathrm{,}{x}_N\right]\times R\to R^m$ is continuous and at least twice differentiable. We seek the solution on equidistant set of points defined on the integration interval $x_n=x_0<x_1<\cdots <x_N=b$ ,

$x_n=x_0+nh$ , $n=\mathrm{0,1,2,}\cdots \mathrm{,}N-1$ , $h=\frac{b-1}{N-1}$ , N is a positive integer.

A potentially good numerical method for the solution of stiff systems must have good accuracy and reasonably wide region of absolute stability. A-Stability requirement is the minimum criteria on the choice of suitable methods. The search for higher order A-stable linear multistep method is carried out in two ways; firstly, the use of higher derivatives of the approximate solution and secondly, the inclusion of additional stages of off grid points Ezzeddine and Hojjati [1] .

Several authors such as Enright [2] , Enright and Pryce [3] , Brown [4] , Cash [5] , Okunuga [6] , Abhilimen and Okunuga [7] , Ngwane and Jator [8] , and Yakubu and Markus [9] have developed second derivative methods for the solution of (1) whose solution has exponential functions.

The aim of this paper is to develop a class of second derivative linear multistep method with varying step-lengths which are A-stable with large region of absolute stability (see Figures 1-3). The three methods recovered are tested on some numerical examples and their results compared with each other in order to determine how to fix the varying step-lengths to obtain the best results as shown Tables 1-4.

2. Development of the Method

We considered the approximate solution of the form

$y\left(x\right)={\displaystyle \underset{n=0}{\overset{j}{\sum }}}{a}_n{x}^n$ (2)

where $a_n$ ’s are constants to be determined. The ith derivative of (2) gives

$y^{\left(i\right)}\left(x\right)={\displaystyle \underset{n-i}{\overset{j}{\sum }}}n\left(n-1\right)\left(n-2\right)\cdots \left(n-i\right)a_n{x}^{n-i}$ (3)

Imposing the following conditions on (2)

$y\left(x_n\right)=y_n,y^{\prime }\left(x_{n+j}\right)=f_{n+j},j=0,1,2,\cdots ,s,y^{″}\left(x_{n+j}\right)=g_{n+j},j=0,1,2,\cdots ,\tau$

evaluating (2) at $x_n$ , first derivative of (3) at $x_{n+j},j=0,1,\cdots ,s$ and second derivative of (3) at $x_{n+j},j=0,1,\cdots ,\tau$ give a system of equations in the form

$AX=U$ (4)

Figure 1. Showing RAS for Case 1.

Figure 2. Showing RAS for Case 11.

Figure 3. Showing RAS for Case III.

where

$X=\left[\begin{array}{cccccc}1& x_n& x_n^2& x_n^3& \cdots & x_n^j\\ 0& 1& 2x_n& 3x_n^2& \cdots & j{x}_n^{j-1}\\ 0& 1& 2x_{n+1}& 3x_{n+1}^2& \cdots & j{x}_{n+1}^{j-1}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 1& 2x_{n+s}& 3x_{n+s}^2& \cdots & j{x}_{n+s}^{j-1}\\ 0& 0& 2& 6x_n& \cdots & j\left(j-1\right)x_n^{j-2}\\ 0& 0& 2& 6x_{n+1}& \cdots & j\left(j-1\right)x_{n+1}^{j-2}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 2& 6x_{n+\tau }& \cdots & j\left(j-1\right)x_{n+\tau }^{j-2}\end{array}\right]$

$A={\left[a_0{a}_1\cdots a_{s+\tau +2}\right]}^{\text{T}},U={\left[y_n{f}_n\cdots f_{n+s}{g}_n\cdots g_{n+\tau }\right]}^{\text{T}}$

Table 1. Showing Results of Example 1.

Table 2. Showing Results for Example 2.

Solving (4) using cramer’s method for the unknown constants, substitute it into (2), we obtain continuous linear multistep method in the form

$y_{n+t}=y_n+h{\displaystyle \underset{j=0}{\overset{s}{\sum }}}{\beta }_j\left(t\right)f_{n+j}+h^2{\displaystyle \underset{j=0}{\overset{\tau }{\sum }}}{\eta }_j\left(t\right)g_{n+j}$ (5)

subject to $h{\displaystyle {\sum }_{j=0}^s}{\beta }_j\left(t\right)=th$ , ${\beta }_j\left(t\right)$ and ${\eta }_j\left(t\right)$ are polynomial of degree $s+\tau +1$ . Evaluating (5) at selected grid points give a block method in the form

${\zeta }^{\left(1\right)}{Y}_{m+1}={\zeta }^{\left(0\right)}{Y}_m+h\left({\eta }^{\left(0\right)}{F}_m+{\eta }^{\left(1\right)}{F}_{m+1}\right)+h^2\left({\gamma }^{\left(0\right)}{G}_m+{\gamma }^{\left(1\right)}{G}_{m+1}\right)$ (6)

where

Table 3. Showing Results for Example 3.

Table 4. Showing Results for Example 4.

$Y_{m+1}={\left[\begin{array}{cccc}{y}_{n+1}& y_{n+2}& \cdots & y_{n+s}\end{array}\right]}^{\text{T}},F_m={\left[\begin{array}{cccc}{f}_{n-1}& f_{n-2}& \cdots & f_n\end{array}\right]}^{\text{T}}$

$F_{m+1}={\left[\begin{array}{cccc}{f}_{n+1}& f_{n+2}& \cdots & f_{n+s}\end{array}\right]}^{\text{T}},Y_m={\left[\begin{array}{cccc}{y}_{n-1}& y_{n-2}& \cdots & y_n\end{array}\right]}^{\text{T}}$

$G_m={\left[\begin{array}{cccc}{g}_{n-1}& g_{n-2}& \cdots & g_n\end{array}\right]}^{\text{T}},G_{m+1}={\left[\begin{array}{cccc}{g}_{n+1}& g_{n+2}& \cdots & g_{n+\tau }\end{array}\right]}^{\text{T}}$

Writing (6) in the form

$F\left(Y_{m+1}\right)=Y_{m+1}-{\zeta }^{\left(0\right)}{Y}_m-h\left({\eta }^{\left(0\right)}{F}_m+{\eta }^{\left(1\right)}{F}_{m+1}\right)-h^2\left({\gamma }^{\left(0\right)}{G}_m+{\gamma }^{\left(1\right)}{G}_{m+1}\right)=0$ (7)

Solving the system of nonlinear equations using Newton Raphson’s method gives

$Y_{m+1}^{\kappa +1}=Y_{m+1}^{\kappa }-{\left(J_{\kappa }\right)}^{-1}F\left(Y_{m+1}^{\kappa }\right)$ (8)

where $J_{\kappa }$ is the Jacobian matrix. The necessary and sufficient condition for convergence of (8) is that the spectral radius of the inverse of the Jacobian matrix $\left|\rho \left(J_{\kappa }^{-1}\right)\right|<1$

$J=\left|\begin{array}{cccccccc}\frac{\partial F\left(y_{n+1}\right)}{\partial y_{n+1}}& \frac{\partial F\left(y_{n+2}\right)}{\partial y_{n+1}}& \cdots & \frac{\partial F\left(y_{n+s}\right)}{\partial y_{n+1}}& \frac{\partial G\left(y_{n+1}\right)}{\partial y_{n+1}}& \frac{\partial G\left(y_{n+1}\right)}{\partial y_{n+1}}& \cdots & \frac{\partial G\left(y_{n+\tau }\right)}{\partial y_{n+1}}\\ \frac{\partial F\left(y_{n+1}\right)}{\partial y_{n+2}}& \frac{\partial F\left(y_{n+2}\right)}{\partial y_{n+2}}& \cdots & \frac{\partial F\left(y_{n+s}\right)}{\partial y_{n+2}}& \frac{\partial G\left(y_{n+1}\right)}{\partial y_{n+2}}& \frac{\partial G\left(y_{n+2}\right)}{\partial y_{n+2}}& \cdots & \frac{\partial G\left(y_{n+\tau }\right)}{\partial y_{n+2}}\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ \frac{\partial F\left(y_{n+s}\right)}{\partial y_{n+s}}& \frac{\partial F\left(y_{n+s}\right)}{\partial y_{n+s}}& \cdots & \frac{\partial F\left(y_{n+s}\right)}{\partial y_{n+s}}& \frac{\partial G\left(y2\right)}{\partial y_{n+s}}& \frac{\partial G\left(y_{n+2}\right)}{\partial y_{n+s}}& \cdots & \frac{\partial G\left(y_{n+\tau }\right)}{\partial y_{n+s}}\end{array}\right|$

Specification of the Method

In this paper, we consider grid points $x_{n+j},j=0<u<v<w$ , hence the (6) reduces to

${\zeta }^{\left(1\right)}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],{\zeta }^{\left(0\right)}=\left[\begin{array}{ccc}0& 0& 1\\ 0& 0& 1\\ 0& 0& 1\end{array}\right],{\eta }^{\left(0\right)}=\left[\begin{array}{ccc}0& 0& {\theta }_{11}\\ 0& 0& {\theta }_{21}\\ 0& 0& {\theta }_{31}\end{array}\right]$

${\eta }^{\left(1\right)}=\left[\begin{array}{ccc}{\theta }_{12}& {\theta }_{13}& {\theta }_{14}\\ {\theta }_{22}& {\theta }_{23}& {\theta }_{24}\\ {\theta }_{32}& {\theta }_{33}& {\theta }_{34}\end{array}\right],{\gamma }^{\left(0\right)}=\left[\begin{array}{ccc}0& 0& {\theta }_{15}\\ 0& 0& {\theta }_{25}\\ 0& 0& {\theta }_{35}\end{array}\right],{\gamma }^{\left(1\right)}=\left[\begin{array}{ccc}{\theta }_{16}& {\theta }_{17}& {\theta }_{18}\\ {\theta }_{26}& {\theta }_{27}& {\theta }_{28}\\ {\theta }_{36}& {\theta }_{37}& {\theta }_{30}\end{array}\right]$

$Y_{m+1}={\left[\begin{array}{ccc}{y}_{n+u}& y_{n+v}& y_{n+w}\end{array}\right]}^{\text{T}},Y_m={\left[\begin{array}{ccc}{y}_{n-1}& y_{n-2}& y_n\end{array}\right]}^{\text{T}},$

$F_m={\left[\begin{array}{ccc}{f}_{n-1}& f_{n-2}& f_n\end{array}\right]}^{\text{T}},F_{m+1}={\left[\begin{array}{ccc}{f}_{n+u}& f_{n+v}& f_{n+w}\end{array}\right]}^{\text{T}},$

$G_m={\left[\begin{array}{ccc}{g}_{n-1}& g_{n-2}& g_n\end{array}\right]}^{\text{T}},G_{m+1}={\left[\begin{array}{ccc}{g}_{n+u}& g_{n+v}& g_{n+w}\end{array}\right]}^{\text{T}}$

${\theta }_{11}=u\frac{\begin{array}{l}[5u^{5}v+5u^{5}w+14u^3{v}^3-16u^4{v}^2+14u^3{w}^3-16u^4{w}^2+210v^3{w}^3-23u^{4}vw\\ -56u{v}^2{w}^3-56u{v}^3{w}^2-28u^{2}v{w}^3-28u^2{v}^{3}w+40u^{3}v{w}^2+40u^3{v}^{2}w]\end{array}}{420v^3{w}^3}$

${\theta }_{12}=u\frac{\begin{array}{l}[-350u^{5}v-350u^{5}w-140u^3{v}^3+388u^4{v}^2-140u^3{w}^3+388u^4{w}^2\\ +210v^3{w}^3+105u^6+1554u^2{v}^2{w}^2+1187u^{4}vw-574u{v}^2{w}^3\\ -574u{v}^3{w}^2+490u^{2}v{w}^3+490u^2{v}^{3}w-1342u^{3}v{w}^2-1342u^3{v}^{2}w]\end{array}}{420{\left(u-w\right)}^3{\left(u-v\right)}^3}$

${\theta }_{13}=-u^5\frac{\begin{array}{l}[10u^{4}v-5u^{4}w+28u^2{v}^3-35u^3{v}^2-14u^2{w}^3+16u^3{w}^2-70v^2{w}^3+98v^3{w}^2\\ +70uv{w}^3-98u{v}^{3}w-10u^{3}vw-42u{v}^2{w}^2-46u^{2}v{w}^2+98u^2{v}^{2}w]\end{array}}{420v^3{\left(v-w\right)}^3{\left(u-v\right)}^3}$

${\theta }_{14}=-u^5\frac{\begin{array}{l}[5u^{4}v-10u^{4}w+14u^2{v}^3-16u^3{v}^2-28u^2{w}^3+35u^3{w}^2{}^3-98v^{2}w+70v^3{w}^2\\ +98uv{w}^3-70u{v}^{3}w+10u^{3}vw+42u{v}^2{w}^2-98u^{2}v{w}^2+46u^2{v}^{2}w]\end{array}}{420w^3{\left(v-w\right)}^3{\left(u-w\right)}^3}$

${\theta }_{15}=u^2\frac{\begin{array}{l}[-16u^{3}v-16u^{3}w+14u^2{v}^2+14u^2{w}^2+70v^2{w}^2\\ +5u^4-56uv{w}^2-56u{v}^{2}w+56u^{2}vw]\end{array}}{840v^2{w}^2}$

${\theta }_{16}=-u^2\frac{\begin{array}{l}[-40u^{3}v-40u^{3}w+28u^2{v}^2+28u^2{w}^2+70v^2{w}^2\\ +15u^4-84uv{w}^2-84u{v}^{2}w+112u^{2}vw]\end{array}}{840{\left(u-w\right)}^2{\left(u-v\right)}^2}$

${\theta }_{17}=u^5\frac{-8u^{2}v+14u{w}^2-16u^{2}w-28v{w}^2+5u^3+28uvw}{840v^2{\left(v-w\right)}^2{\left(u-v\right)}^2}$

${\theta }_{18}=u^5\frac{14u{v}^2-16u^{2}v-8u^{2}w-28v^{2}w+5u^3+28uvw}{840w^2{\left(v-w\right)}^2{\left(u-w\right)}^2}$

${\theta }_{21}=v\frac{\begin{array}{l}[5u{v}^5+5v^{5}w-16u^2{v}^4+14u^3{v}^3+210u^3{w}^3+14v^3{w}^3\\ -16v^4{w}^2-23u{v}^{4}w-28u{v}^2{w}^3+40u{v}^3{w}^2-56u^{2}v{w}^3\\ +40u^2{v}^{3}w-56u^{3}v{w}^2-28u^3{v}^{2}w]\end{array}}{420u^3{w}^3}$

${\theta }_{22}=v^5\frac{\begin{array}{l}[10u{v}^4-5v^{4}w-35u^2{v}^3+28u^3{v}^2-70u^2{w}^3+98u^3{w}^2-14v^2{w}^3+16v^3{w}^2\\ +70uv{w}^3-10u{v}^{3}w-98u^{3}vw-46u{v}^2{w}^2-42u^{2}v{w}^2+98u^2{v}^{2}w]\end{array}}{420u^3{\left(u-w\right)}^3{\left(u-v\right)}^3}$

${\theta }_{23}=v\frac{\begin{array}{l}[350u{v}^5+350v^{5}w-388u^2{v}^4+140u^3{v}^3-210u^3{w}^3+140v^3{w}^3\\ -388v^4{w}^2-105v^6-1554u^2{v}^2{w}^2-1187u{v}^{4}w-490u{v}^2{w}^3\\ +1342u{v}^3{w}^2+574u^{2}v{w}^3+1342u^2{v}^{3}w+574u^{3}v{w}^2-490u^3{v}^{2}w]\end{array}}{420{\left(v-w\right)}^3{\left(u-v\right)}^3}$

${\theta }_{24}=-v^5\frac{\begin{array}{l}[5u{v}^4-10v^{4}w-16u^2{v}^3+14u^3{v}^2-98u^2{w}^3+70u^3{w}^2-28v^2{w}^3+35v^3{w}^2\\ +98uv{w}^3+10u{v}^{3}w-70u^{3}vw-98u{v}^2{w}^2+42u^{2}v{w}^2+46u^2{v}^{2}w]\end{array}}{420w^3{\left(v-w\right)}^3{\left(u-w\right)}^3}$

${\theta }_{25}=v^2\frac{\begin{array}{l}[-16u{v}^3-16v^{3}w+14u^2{v}^2+70u^2{w}^2+14v^2{w}^2\\ +5v^4-56uv{w}^2+56u{v}^{2}w-56u^{2}vw]\end{array}}{840u^2{w}^2}$

${\theta }_{26}=-v^5\frac{8u{v}^2+28u{w}^2-14v{w}^2+16v^{2}w-5v^3-28uvw}{840u^2{\left(u-w\right)}^2{\left(u-v\right)}^2}$

${\theta }_{27}=-v^2\frac{\begin{array}{l}[-40u{v}^3-40v^{3}w+28u^2{v}^2+70u^2{w}^2+28v^2{w}^2+15v^4\\ -84uv{w}^2+112u{v}^{2}w-84u^{2}vw]\end{array}}{840{\left(v-w\right)}^2{\left(u-v\right)}^2}$

${\theta }_{28}=v^5\frac{-16u{v}^2+14u^{2}v-28u^{2}w-8v^{2}w+5v^3+28uvw}{840w^2{\left(v-w\right)}^2{\left(u-w\right)}^2}$

${\theta }_{31}=w\frac{\begin{array}{l}[5u{w}^5+5v{w}^5+210u^3{v}^3-16u^2{w}^4+14u^3{w}^3-16v^2{w}^4+14v^3{w}^3-23uv{w}^4\\ +40u{v}^2{w}^3-28u{v}^3{w}^2+40u^{2}v{w}^3-56u^2{v}^{3}w-28u^{3}v{w}^2-56u^3{v}^{2}w]\end{array}}{420u^3{v}^3}$

${\theta }_{32}=w^5\frac{\begin{array}{l}[10u{w}^4-5v{w}^4-70u^2{v}^3+98u^3{v}^2-35u^2{w}^3+28u^3{w}^2+16v^2{w}^3-14v^3{w}^2\\ -10uv{w}^3+70u{v}^{3}w-98u^{3}vw-46u{v}^2{w}^2+98u^{2}v{w}^2-42u^2{v}^{2}w]\end{array}}{420u^3{\left(u-w\right)}^3{\left(u-v\right)}^3}$

${\theta }_{33}=w^5\frac{\begin{array}{l}[5u{w}^4-10v{w}^4-98u^2{v}^3+70u^3{v}^2-16u^2{w}^3+14u^3{w}^2+35v^2{w}^3-28v^3{w}^2\\ +10uv{w}^3+98u{v}^{3}w-70u^{3}vw-98u{v}^2{w}^2+46u^{2}v{w}^2+42u^2{v}^{2}w]\end{array}}{420v^3{\left(v-w\right)}^3{\left(u-v\right)}^3}$

${\theta }_{34}=w\frac{\begin{array}{l}[-350u{w}^5-350v{w}^5+210u^3{v}^3+388u^2{w}^4-140u^3{w}^3+388v^2{w}^4\\ -140v^3{w}^3+105w^6+1554u^2{v}^2{w}^2+1187uv{w}^4-1342u{v}^2{w}^3\\ +490u{v}^3{w}^2-1342u^{2}v{w}^3-574u^2{v}^{3}w+490u^{3}v{w}^2-574u^3{v}^{2}w]\end{array}}{420{\left(v-w\right)}^3{\left(u-w\right)}^3}$

${\theta }_{35}=w^2\frac{\begin{array}{l}[-16u{w}^3-16v{w}^3+70u^2{v}^2+14u^2{w}^2+14v^2{w}^2\\ +5w^4+56uv{w}^2-56u{v}^{2}w-56u^{2}vw]\end{array}}{840u^2{v}^2}$

${\theta }_{36}=-w^5\frac{28u{v}^2+8u{w}^2+16v{w}^2-14v^{2}w-5w^3-28uvw}{840u^2{\left(u-w\right)}^2{\left(u-v\right)}^2}$

${\theta }_{37}=-w^5\frac{28u^{2}v+16u{w}^2-14u^{2}w+8v{w}^2-5w^3-28uvw}{840v^2{\left(v-w\right)}^2{\left(u-v\right)}^2}$

${\theta }_{38}=-w^2\frac{\begin{array}{l}[-40u{w}^3-40v{w}^3+70u^2{v}^2+28u^2{w}^2+28v^2{w}^2\\ +15w^4+112uv{w}^2-84u{v}^{2}w-84u^{2}vw]\end{array}}{840{\left(v-w\right)}^2{\left(u-w\right)}^2}$

For our methods, in case 1, we considered one step method with two hybrid

points with equal interval where $u=\frac{1}{3},v=\frac{2}{3},w=1$ . For case II, we considered

two step method with one hybrid points with equal interval, where $u=1,v=3/2,w=2$ and for the lase case III, we considered three step method with equal interval where $u=1,v=2,w=3$ .

3. Stability Properties

In this section, we investigate the basic properties of the developed method vis-a-vis order, local truncation error, consistency, zero-stability, convergence, and region of absolute stability of the methods.

3.1. Order of Convergence

The operation l is associated with the linear method defined by

$\mathcal{l}\left[y\left(x\right)\mathrm{<mo>\; :\; </mo>}h\right]={\zeta }^{\left(1\right)}{Y}_{m+1}-{\zeta }^{\left(0\right)}{Y}_m-h\left({\eta }^{\left(0\right)}{F}_m+{\eta }^{\left(1\right)}{F}_{m+1}\right)-h^2\left({\gamma }^{\left(0\right)}{G}_m+{\gamma }^{\left(1\right)}{G}_{m+1}\right)$ (9)

where $y\left(x\right)$ is an arbitrary function, continuously differentiable on an interval $\left[x_n\mathrm{,}{x}_N\right]$ . Ehigie, J. O. and Okunuga [11] can be written in Taylor expansion as

$\mathcal{l}\left[y\left(x\right)\mathrm{<mo>\; :\; </mo>}h\right]=c_{0}y\left(x_n\right)+c_{1}hy\text{'}\left(x_n\right)+c_2{h}^2{y}^{″}\left(x_n\right)+\cdots +c_q{h}^q{y}^{\left(q\right)}\left(x_n\right)+\cdots$

where

$c_p=\frac{1}{p}\left[{\displaystyle \underset{j=1}{\overset{r}{\sum }}}{j}^p{\theta }_j-\frac{1}{\left(p-1\right)!}{\displaystyle \underset{j=1}{\overset{r}{\sum }}}{j}^{p-1}{\gamma }_j\right]$

(9) is of order p if

$\mathcal{l}\left[y\left(x\right):h\right]=o\left(h^{p+1}\right),c_0=c_1=\cdots =c_p\ne c_{p+1}=0$

$c_{p+1}$ is called the error constant and $c_{p+1}{h}^{p+1}{y}^{\left(p+1\right)}\left(x\right)$ is called the local truncation error (LTE).

For our scheme, the order is 8 with error constant given by

$\left[\begin{array}{c}\frac{h^9}{50803200}{u}^5\left(-15u^{3}v-15u^{3}w+12u^2{v}^2+12u^2{w}^2+42v^2{w}^2+5u^4-42uv{w}^2-42u{v}^{2}w+48u^{2}vw\right)\\ \frac{h^9}{50803200}{v}^5\left(-15u{v}^3-15v^{3}w+12u^2{v}^2+42u^2{w}^2+12v^2{w}^2+5v^4-42uv{w}^2+48u{v}^{2}w-42u^{2}vw\right)\\ \frac{h^9}{50803200}{w}^5\left(-15u{w}^3-15v{w}^3+42u^2{v}^2+12u^2{w}^2+12v^2{w}^2+5w^4+48uv{w}^2-42u{v}^{2}w-42u^{2}vw\right)\end{array}\right]$

3.2. Consistency

Definition 1 A block method is consistent if it has order $p\ge 1.$

3.3. Zero-Stability

Definition 2 A method is said to be Zero-stable if no root of the first characteristics polynomial has modulus greater than one, and if every root of modulus one has multiplicity not greater than one or is simple.

$\rho \left(\lambda \right)=\left|\lambda \left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]-\left[\begin{array}{ccc}0& 0& 1\\ 0& 0& 1\\ 0& 0& 1\end{array}\right]\right|$

The roots of the determinant gives $\lambda =0,0,1$ , hence the method is consistent.

3.4. Convergence

Definition 3 A method is said to be convergent if It is consistent and zero stable.

3.5. Region of Absolute Stability (RAS)

Definition 4 The Region of absolute stability (RAS) of a L M M is the set R = { $z=\lambda h$ : for z where the root of the stability polynomial are absolute less than one}.

Substituting the test equation $y^{\prime }=\lambda y,y^{″}={\lambda }^{2}y$ and writing $r={\text{e}}^{i\theta }$ , the region of absolute stability for our new methods are shown below:

In determining the region of absolute stability, we consider three cases in this paper as described below:

Definition 5 A-Stability

A method is said to be A-stable if its region of absolute stability contains the whole of the complex left hand-half plane $\mathrm{Re}\left(h\lambda \right)<0$ . Alternatively, a numerical method is called A-Stable if the solution tend to zero as $n\to \infty$ when the method is applied with fixed h to any differential equation of the form

$\frac{\partial y}{\partial x}=\lambda y$ , where l is a complex constant with negative real part. Hence, for our cases I, II, & III; The regions of absolute stabilities are given in Figures 1-3.

We conclude that the three cases presented in this paper are A-stable.

4. Numerical Experiments

In this section we considered four examples to test the efficiency of the method. We compared the results of the cases in order to conclude on the best way to fix u,v and w. The following notations are used to show the results $\psi \psi e-\upsilon \upsilon =\psi \psi \ast {10}^{-\upsilon \upsilon }.$

Example 1 We consider a system in the range $0\le x\le 20$

$y^{\prime }=\left(\begin{array}{cc}-1& 95\\ -1& -97\end{array}\right)y,y\left(x\right)=(11)$

with the exact solution

$y\left(x\right)=\frac{1}{47}\left(\begin{array}{c}95{\text{e}}^{-2x}-48{\text{e}}^{-95x}\\ 48{\text{e}}^{-96x}-{\text{e}}^{-2x}\end{array}\right)$

The eigenvalues of the Jacobian matrix are ${\lambda }_1=-2,{\lambda }_2=-96$ with the stiffness ratio 1:48. Source :Abdulimen [10] .

Example 2 We consider a simple nonlinear stiff system

$\begin{array}{l}{y^{\prime }}_1=-y_1,y_1\left(0\right)=5\\ {y^{\prime }}_2=y_1^2-2y_2,y_2\left(0\right)=5\end{array}$

$x\in \left[\mathrm{0,20}\right]$ , with the exact solution

$\left(\begin{array}{c}{y}_1\\ y_2\end{array}\right)=\left(\begin{array}{c}5exp\left(-x\right)\\ 5exp\left(-2x\right)\left(\mathrm{1\; <mo>\; \_\; </mo>\; 5}x\right)\end{array}\right)$

The eigenvalues of the Jacobian matrix are ${\lambda }_1=-1,{\lambda }_2=-2$ . Source: Yakubu and Markus [9] .

Example 3 In the this example we consider stiff nonlinear system of two dimensional Kaps problem with corresponding initial conditions

$\left[\begin{array}{c}{y^{\prime }}_1\left(x\right)\\ {y^{\prime }}_2\left(x\right)\end{array}\right]=\left[\begin{array}{c}-1002y_1\left(x\right)+1000y_2{\left(x\right)}^2\\ y_1\left(x\right)-y\left(x\right)\left(1+y_2\left(x\right)\right)\end{array}\right],\left[\begin{array}{c}{y}_1\left(0\right)\\ y_2\left(0\right)\end{array}\right]=\left[\begin{array}{c}1\\ 1\end{array}\right]$

The exact solution is

$\left[\begin{array}{c}{y}_1\left(x\right)\\ y_2\left(x\right)\end{array}\right]=\left[\begin{array}{c}exp\left(-2x\right)\\ exp\left(-x\right)\end{array}\right]$

Source: Yakubu and Markus [9] .

Example 4 We consider a four dimensional problems

$\left(\begin{array}{c}{y^{\prime }}_1\left(x\right)\\ {y^{\prime }}_2\left(x\right)\\ {y^{\prime }}_3\left(x\right)\\ {y^{\prime }}_4\left(x\right)\end{array}\right)=\left(\begin{array}{c}-{10}^4{y}_1\left(x\right)+100y_2\left(x\right)-10y_3\left(x\right)+y_4\left(x\right)\\ -1000y_2\left(x\right)+10y_3\left(x\right)-10y_4\left(x\right)\\ -y_3\left(x\right)+10y_4\left(x\right)\\ -0.1y_4\left(x\right)\end{array}\right)=\left(\begin{array}{c}{y}_1\left(0\right)\\ y_2\left(0\right)\\ y_3\left(0\right)\\ y_4\left(0\right)\end{array}\right)=\left(\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right)$

Within the range $0\le x\le 20$ . The eigenvalues of the Jacobian matrix ${\lambda }_1=-0.1$ , ${\lambda }_2=-10$ , ${\lambda }_3=-1000$ and ${\lambda }_4=-10000$ . The exact solution is given as

$\begin{array}{c}{y}_1\left(x\right)=-\frac{89990090}{8999010009}{\text{e}}^{-0.1x}+\frac{818090}{89901009}{\text{e}}^{-x}\\ +\frac{998911}{899010090}{\text{e}}^{-1000x}+\frac{89071119179}{89990100090}{\text{e}}^{-10000x}\end{array}$

$y_2\left(x\right)=\frac{9100}{8991}{\text{e}}^{-0.1x}-\frac{910}{8991}{\text{e}}^{-x}+\frac{9989911}{9989001}{\text{e}}^{-1000x}$

$y_3\left(x\right)=\frac{100}{9}{\text{e}}^{-0.1x}-\frac{91}{9}{\text{e}}^{-x}$

$y_4\left(x\right)={\text{e}}^{-0.1x}$

Source: Ehigie, Okunuga, and Sofoluwe [11] .

5. Conclusion

In this paper, we introduced three new second derivative linear multistep methods for the numerical solution of stiff initial value problems. Four numerical examples were considered, the results justified the proficiency of the second derivative method which is cheaper to implement since it does not require starting values and particularly the new methods show that lower step method gives better accuracy than higher step methods of the same order. We are able to achieve the aim of this paper which is to develop a class of second derivative linear multistep methods that are A-stable with large region of absolute stabilityas shown in Figures 1-3. The results from the high-order methods are very encourging (see Tables 1-4), therefore, we recommend further investigation of the second-derivative.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Ezzeddine, A.K. and Hojjati, G. (2012) Third Derivative Multistep Methods for Stiff Systems. International Journal of Non-Linear Sciences, 14, 443-450.
[2]	Enright, W.H. (1974) Second Derivative Multistep Methods for Stiff Differential Equations. SIAM J. Numerical Analysis, 11, 321-331. https://doi.org/10.1137/0711029
[3]	Enright, W.H. and Priyce, J.O. (1983) Two FORTRAN Packages of Assessing IV Methods. Technical Reports No. 16183, Development of Computer Sciences University of Toronto, Canada.
[4]	Brown, T.I. (1977) Some Characteristics of Implicit Multistep Multi-Derivative Integration Formulas. SIAM J. Numer. Math., 34, 59-84.
[5]	Cash, J.R. (1981) On the Exponential Fitting of Composite Multi-Derivative Linear Multi Step Method. SIAM J. Numer. Anal, 18, 808-812. https://doi.org/10.1137/0718055
[6]	Okunuga, S.A. (1997) Fourth Order Composite Two Step Method for Stiff Problems. Int. J. Comput. Math., 2, 39-47.
[7]	Abhulimen, C.E. and Okunuga, S.A. (2008) Exponentially Fitted Second Derivative Multistep Method for Stiff Initial Value Problems for ODEs. J. Eng.Sci. Appl., 5, 36-47.
[8]	Ngwane, F.F. and Jator, S.N. (2012) Block Hybrid Second Derivative Method for Stiff Systems. International Journal of Pure and Applied Mathematics, 80, 543-559.
[9]	Yakubu, D.G. and Markus, S. (2016) The Efficiency of Efficiency of Second Derivative Multistep Methods for the Numerical Integration of Stiff Systems. Journal of the Nigeria Mathematics Society, 35, 107-127.
[10]	Abdulimen, C.E. (2014) Exponentially Fitted Third Derivative Three Step Method for the Numerical Integration of Stiff IVP. Applied Mathematics and Computation, 243, 446-453.
[11]	Ehigie, J.O., Okunuga, S.A. and Sofoluwe, A.B. (2013) A Class of Exponetially Fitted Second Derivative Fitted Second Derivative Extended Backward Differentiation Formula for Solving Stiff Problems. Fascicuci Mathematics, 1, NR 51.