Stability Analysis of a Numerical Integrator for Solving First Order Ordinary Differential Equation

Samuel Olukayode Ayinde^1*, Adesoji Abraham Obayomi¹, Funmilayo Sarah Adebayo²

¹Department of Mathematics, Faculty of Science, Ekiti State University, Ado Ekiti, Nigeria.
²Federal Polytechnic, Ado Ekiti, Nigeria.

DOI: 10.4236/jamp.2017.511179                

Abstract

In this paper, we used an interpolation function to derive a Numerical Integrator that can be used for solving first order Initial Value Problems in Ordinary Differential Equation. The numerical quality of the Integrator has been analyzed to authenticate the reliability of the new method. The numerical test showed that the finite difference methods developed possess the same monotonic properties with the analytic solution of the sampled Initial Value Problems.

Keywords

Numerical Integrator, Autonomous and Non-Autonomous, Ordinary Differential Equation, Initial Value Problems, Stability Analysis

1. Introduction

Many Scholars have derived various Numerical Integrators using various techniques including interpolating functions that include the work of [1] [2] [3] [4] among others. All these authors have employed some analytically continuous functions to create numerically stable Integrators that can be used for ordinary differential equations. In this work we use an analytically differentiable interpolating function to create a one-step Finite Difference scheme for solving Initial Value Problems of first order Ordinary Differential Equations, and we are considering the concept of Nature, of Solutions of first order Ordinary Differential Equations to assume a theoretical solution and use that assumption to derive a discrete model that can be applied to some Ordinary differential equations.

Definition 1 [5]

Consider the nth-order ordinary differential equation

$F\left(x,y,y^1,\cdots ,y^n\right)=0$ (1)

where F is a real function of its (n + 2) arguments $x,y,y^1,\cdots ,y^n.$

1) Let f be a real function defined for all x in a real interval I and having an nth derivative (and hence also all lower ordered derivatives) for all $x\in \text{I}$ . The function f is called an explicit solution of the differential Equation (1) on interval I if it fulfills the following two requirements.

$F\left(x,f\left(x\right),f^1\left(x\right),\cdots ,f^n\left(x\right)\right)$ (2)

is defined for all $x\in \text{I}$ , and

$F\left(x,f\left(x\right),f^1\left(x\right),\cdots ,f^n\left(x\right)\right)=0$ (3)

for all $x\in \text{I}$ .

That is, the substitution of $f\left(x\right)$ and its various derivatives for y and its corresponding derivatives.

2) A relation $g\left(x,y\right)=0$ is called an implicit solution of (1) if this relation defines at least one real function f of the variable x on an internal I such that this function is an explicit solution of (1) on this interval.

3) Both explicit solutions and implicit solutions will usually be called simply Solutions.

We now consider the geometric significance of differential equations and their solutions. We first recall that a real function F(x) may be represented geometrically by a Curve $y=F\left(x\right)$ in the xy plane and that the value of the derivative of F at x, $F^{\prime }\left(x\right)$ , may be interpreted as the slope of the curve $y=F\left(x\right)$ at x.

2. Formulation of the Interpolating Function

Consider the initial value problem of the IVP

$y^{\prime }\left(x\right)=f\left(x,y\right),y\left(x_0\right)=\eta ,$ (4)

where $\eta$ is a discrete variables in the interval $\left[x_n,x_{n+1}\right]$ . In this we consider the method based on local representation of the theoretical solution $y\left(x\right)$ .

Let us assume that the theoretical solution $y\left(x\right)$ to the initial value problem 4) can be locally represented in the interval $\left[x_n,x_{n+1}\right],n\ge 1$ by the non-polynomial interpolating function given by:

$F\left(x\right)={\alpha }_1{\text{e}}^{-2x}+{\alpha }_2{x}^2+{\alpha }_{3}x+{\alpha }_4$ (5)

where ${\alpha }_1,{\alpha }_2$ and ${\alpha }_3$ are real undetermined coefficients, and ${\alpha }_4$ is a constant.

3. Derivation of the Integrator

We assumed that the theoretical solution $y\left(x\right)$ to the initial value problem (5) can be locally represented in the interval $\left[x_n,x_{n+1}\right],n\ge 0$ by the non-polynomial interpolating function;

$F\left(x\right)={\alpha }_1{\text{e}}^{-2x}+{\alpha }_2{x}^2+{\alpha }_{3}x+{\alpha }_4$ (6)

where ${\alpha }_1,{\alpha }_2$ and ${\alpha }_3$ are real undetermined coefficients, and ${\alpha }_4$ is a constant.

We shall assume $y_n$ is a numerical estimate to the theoretical solution $y\left(x\right)$ and $f_n=f\left(x_n,y_n\right)$ .

We define mesh points as follows:

$x_n=a+nh,n=0,1,2,\cdots$ (7)

We impose the following constraints on the interpolating function (6) in order to get the undetermined coefficients:

1) The interpolating function must coincide with the theoretical solution at $x=x_n$ and $x=x_{n+1}$ . Hence we required that

$F\left(x_n\right)={\alpha }_1{\text{e}}^{-2x_n}+{\alpha }_2{x}_n^2+{\alpha }_3{x}_n+{\alpha }_4$ (8)

$F\left(x_{n+1}\right)={\alpha }_1{\text{e}}^{-2x_{n+1}}+{\alpha }_2{x}_{n+1}^2+{\alpha }_3{x}_{n+1}+{\alpha }_4$ (9)

2) The derivatives of the interpolating function are required to coincide with the differential equation as well as its first, second, and third derivatives with respect to x at $x=x_n$ .

We denote the i-th derivatives of $f\left(x,y\right)$ with respect to x with $f^{\left(i\right)}$ such that

$F^1\left(x_n\right)=f_n,F^2\left(x_n\right)=f_n^1,F^3\left(x_n\right)=f_n^2,$ (10)

This implies that,

$f_n=-2{\alpha }_1{\text{e}}^{-2x_n}+2{\alpha }_2{x}_n+{\alpha }_3$ (11)

$f_n^1=4{\alpha }_1{\text{e}}^{-2x_n}+2{\alpha }_2$ (12)

$f_n^2=-8{\alpha }_1{\text{e}}^{-2x_n}$ (13)

Solving for ${\alpha }_1,{\alpha }_2$ and ${\alpha }_3$ from Equations (11) (12) and (13), we have

${\alpha }_1=-\frac{1}{8}{f}_n^2{\text{e}}^{2x_n}$ (14)

${\alpha }_2=\frac{1}{2}\left(f_n^1+\frac{1}{2}{f}_n^2\right)$ (15)

and

${\alpha }_3=\left(f_n-\frac{1}{4}{f}_n^2\right)-\left(f_n^1+\frac{1}{2}{f}_n^2\right)x_n$ (16)

Since $F\left(x_{n+1}\right)=y\left(x_{n+1}\right)$ and $F\left(x_n\right)=y(xn)$

Implies that $y\left(x_{n+1}\right)=y_{n+1}$ and $y\left(x_n\right)=y_n$

$F\left(x_{n+1}\right)-F\left(x_n\right)=y_{n+1}-y_n$ (17)

Then we shall have from (8) and (9) into (17)

$y_{n+1}-y_n={\alpha }_1\left[{\text{e}}^{-2x_{n+1}}-{\text{e}}^{-2x_n}\right]+{\alpha }_2\left[x_{n+1}^2-x_n^2\right]+{\alpha }_3\left[x_{n+1}-x_n\right]$ (18)

Recall that $x_n=a+nh$ , $x_{n+1}=a+\left(n+1\right)h$ with $n=0,1,2,\cdots$

Substitute (14) (15) (16), into (18), and simplify we have the integrator

$y_{n+1}=y_n-\frac{1}{8}{f}_n^2\left({\text{e}}^{-2h}-1\right)+\frac{1}{2}\left(f_n^1+\frac{1}{2}{f}_n^2\right)h^2+\left(f_n-\frac{1}{4}{f}_n^2\right)h$ (19)

for solution of the first order differential equation.

4. Properties of the Integration Method

4.1. Qualitative Properties of the Scheme

4.1.1. Definition 2 [6]

Define any algorithm for solving differential equations in which the approximation $y_{n+1}$ to the solution at the point $x_{n+1}$ can be calculated if only $x_n,y_n$ and h are known as one-step method. It is a common practice to write the functional dependence, $y_{n+1}$ , on the quantities $x_n,y_n$ and h in the form:

$y_{n+1}=y_n+h\varnothing \left(x_n,y_n;h\right)$ (20)

where $\varnothing \left(x_n,y_n;h\right)$ is the increment function.

The numerical integrator can be expressed as a one-step method in the form (20) above thus:

From (19) i.e. $y_{n+1}=y_n-\frac{1}{8}{f}_n^2\left({\text{e}}^{-2h}-1\right)+\frac{1}{2}\left(f_n^1+\frac{1}{2}{f}_n^2\right)h^2+\left(f_n-\frac{1}{4}{f}_n^2\right)h$

Expanding ${\text{e}}^{-2h}$ into the fourth term, we have

${\text{e}}^{-2h}={\displaystyle {\sum }_{r=0}^{\infty }\frac{{\left(-2h\right)}^r}{r!}}=1-2h+\frac{{\left(2h\right)}^2}{2!}-\frac{{\left(2h\right)}^3}{3!}+\cdots$ (21)

Put (21) into (19), then expand

$y_{n+1}=y_n+f_{n}h+f_n^1{x}_{n}h+\frac{1}{2}{f}_n^1{h}^2+\frac{1}{2}{f}_n^2{x}_{n}h+\frac{1}{6}{f}_n^2{h}^3$ (22)

$=y_n+h\left\{f_n+f_n^1\left(x_n+\frac{1}{2}h\right)+f_n^2\left(x_n+\frac{1}{6}{h}^2\right)\right\}$ (23)

Let $A=x_n+\frac{1}{2}h$ and $B=x_n+\frac{1}{6}{h}^2$ (24)

Thus our integrator (19) can be written compactly as

$y_{n+1}=y_n+h\left\{f_n+A{f}_n^1+B{f}_n^2\right\}$ (25)

Which is in the form

$y_{n+1}=y_n+h\varnothing \left(x_n,y_n;h\right)$ (26)

where $\varnothing \left(x_n,y_n;h\right)=\left\{f_n+A{f}_n^1+B{f}_n^2\right\}$ (27)

4.1.2. Theorem 1. [7]

Let the increment function of the method defined by (25) be continuous as a function of its arguments in the region defined by

$x\in \left[a,b\right],y\in \left(-\infty ,\infty \right);0\le h\le h_o,$

where $h_o>0$ , and let there exists a constant $L$ such that

$\left|\varnothing \left(x_n,y_n^{*};h\right)-\varnothing \left(x_n,y_n;h\right)\right|\le L\left|y_n^{*}-y_n\right|$ (28)

for all $\left(x_n,y_n;h\right)$ and $\left(x_n,y_n^{*};h\right)$ in the region just defined. Then the relation (28) is the Lipscitz condition and it is the necessary and sufficient condition for the convergence of our method (19).

We shall proof that (19) satisfies (28) in line with the established Fatunla’s theorem.

4.1.3. Proof of Convegence of the Integrator

The increment function $\varnothing \left(x_n,y_n;h\right)$ can be written in the form

$\varnothing \left(x_n,y_n;h\right)=\left\{f\left(x_n,y_n\right)+A{f}^{\left(1\right)}\left(x_n,y_n\right)+B{f}^{\left(2\right)}\left(x_n,y_n\right)\right\}$ (29)

where A and B are constants defined below.

$A=x_n+\frac{1}{2}h$

and

$B=x_n+\frac{1}{6}{h}^2$

Consider Equation (29), we can also write

$\varnothing \left(x_n,y_n^{*};h\right)=\left\{f\left(x_n,y_n^{*}\right)+A{f}^{\left(1\right)}\left(x_n,y_n^{*}\right)+B{f}^{\left(2\right)}\left(x_n,y_n^{*}\right)\right\}$

$\begin{array}{l}\varnothing \left(x_n,y_n^{*};h\right)-\varnothing \left(x_n,y_n;h\right)\\ =f\left(x_n,y_n^{*}\right)-f\left(x_n,y_n\right)+A\left[f^{\left(1\right)}\left(x_n,y_n^{*}\right)-f^{\left(1\right)}\left(x_n,y_n\right)\right]\\ +B\left[f^{\left(2\right)}\left(x_n,y_n^{*}\right)-f^{\left(2\right)}\left(x_n,y_n^{*}\right)\right]\end{array}$ (30)

Let $\bar{y}$ be defined as a point in the interior of the interval whose points are $y$ and $y^{*}$ , applying mean value theorem, we have

$\begin{array}{l}f\left(x_n,y_n^{*}\right)-f\left(x_n,y_n\right)=\frac{\partial f\left(x_n,\bar{y}\right)}{\partial y_n}\left(y_n^{*}-y_n\right)\\ f^{\left(1\right)}\left(x_n,y_n^{*}\right)-f^{\left(1\right)}\left(x_n,y_n\right)=\frac{\partial f^{\left(1\right)}\left(x_n,\bar{y}\right)}{\partial y_n}\left(y_n^{*}-y_n\right)\\ \text{and}{f}^{\left(2\right)}\left(x_n,y_n^{*}\right)-f^{\left(2\right)}\left(x_n,y_n\right)=\frac{\partial f^{\left(2\right)}\left(x_n,\bar{y}\right)}{\partial y_n}\left(y_n^{*}-y_n\right)\end{array}\}$ (31)

We define

$\begin{array}{l}L={\sup }_{\left(x_n,y_n\right)\in D}\frac{\partial f\left(x_n,y_n\right)}{\partial y_n}\\ L_1={\sup }_{\left(x_n,y_n\right)\in D}\frac{\partial f^{\left(1\right)}\left(x_n,y_n\right)}{\partial y_n}\\ \text{and}{L}_2={\sup }_{\left(x_n,y_n\right)\in D}\frac{\partial f^{\left(2\right)}\left(x_n,y_n\right)}{\partial y_n}\end{array}\}$

Therefore

$\begin{array}{l}\varnothing \left(x_n,y_n^{*};h\right)-\varnothing \left(x_n,y_n;h\right)\\ =\frac{\partial f\left(x_n,\bar{y}\right)}{\partial y_n}\left(y_n^{*},y_n\right)+A\left\{\frac{\partial f^{\left(1\right)}\left(x_n,\bar{y}\right)}{\partial y_n}\left(y_n^{*},y_n\right)\right\}\\ +B\left\{\frac{\partial f^2\left(x_n,\bar{y}\right)}{\partial y_n}\left(y_n^{*},y_n\right)\right\}\\ =L\left(y_n^{*}-y_n\right)+A{L}_1\left(y_n^{*}-y_n\right)+B{L}_2\left(y_n^{*}-y_n\right)\end{array}$ (32)

Taking the absolute value of both sides

$\begin{array}{l}\left|\varnothing \left(x_n,y_n^{*};h\right)-\varnothing \left(x_n,y_n;h\right)\right|\\ \le |L\left(y_n^{*}-y_n\right)+A{L}_1\left(y_n^{*}-y_n\right)+B{L}_2\left(y_n^{*}-y_n\right)|\\ \le \left|L+A{L}_1+B{L}_2\right|\left|y^{*}-y\right|\end{array}$ (33)

If we let $M=\left|L+A{L}_1+B{L}_2\right|$

then our Equation (33) turns to

$\left|\varnothing \left(x_n,y_n^{*};h\right)-\varnothing \left(x_n,y_n;h\right)\right|\le M\left|y^{*}-y\right|$ (34)

which is the condition for convergence.

4.2. Consistence of the Integrator

Definition 3 [8]

The integration scheme: $y_{n+1}=y_n+h\left(x_n,y_n;h\right)$ is said to be consistent with the initial-value problem $y^{\prime }\left(x\right)=f\left(x,y\left(x\right)\right)$ , $y\left(a\right)=y_o$ , $x\in \left[a,b\right]$ , $y\in R$ provided the increment function $\varnothing \left(x,y;h\right)$ satisfies the following relationship

$\varnothing \left(x,y;h\right)=f\left(x,y\right)$ (35)

The significance of the consistency of a formula is that it ensures that the method approximates the ordinary differential equation in its place.

Therefore from

$y_{n+1}=y_n+h\left\{f_n+f_n^1\left(x_n+\frac{1}{2}h\right)+f_n^2\left(x_n+\frac{1}{6}{h}^2\right)\right\}$ (36)

where $y_{n+1}=y_n+\theta \left(x_n,y_n;h\right)$ then

$\theta \left(x_n,y_n;h\right)=h\left\{f\left(x_n,y_n\right)+A{f}^{\left(1\right)}\left(x_n,y_n\right)+B{f}^{\left(2\right)}\left(x_n,y_n\right)\right\}$

and

$A=x_n+\frac{1}{2}h,B=x_n+\frac{1}{6}{h}^2$

If $h=0$ , then (36) reduced to $y_{n+1}=y_n$

$\Rightarrow \theta \left(x_n,y_n;0\right)=f\left(x,y\right)$ (37)

It is a known fact that a consistent method has order of at least one [9] . Therefore, the new numerical integrator is consistent since Equation (36) can be reduced to (37) when $h=0$ .

4.3. Stability Analysis of the Integration Method

We shall establish the stability analysis of the integrator by considering the theorem established by Lambert 1972.

Let $y_n=y\left(x_n\right)$ and $P_n=P\left(x_n\right)$ denote two different numerical solutions of initial value problem of ordinary differential Equation (35) with the initial conditions specified as $y\left(x_o\right)=\eta$ and $p\left(x_o\right)={\eta }^{*}$ respectively, such that $\left|\eta -{\eta }^{*}\right|<\epsilon$ , $\epsilon >0$ . If the two numerical estimates are generated by the integrator (19). From the increment function (26), we have

$y_{n+1}=y_n+h\varnothing \left(x_n,y_n;h\right)$ (38)

$P_{n+1}=P_n+h\varnothing \left(x_n,p_n;h\right)$ (39)

The condition that

$\left|y_{n+1}-P_{n+1}\right|\le K\left|\eta -{\eta }^{*}\right|$ (40)

is the necessary and sufficient condition that our new method (19) be stable and convergent.

Proof

From (27) we have

$y_{n+1}=y_n+h\left\{f_n+A{f}_n^1+B{f}_n^2\right\}$ (41)

Then let

$y_{n+1}=y_n+h\left\{f\left(x_n,y_n\right)+A{f}^1\left(x_n,y_n\right)+B{f}^2\left(x_n,y_n\right)\right\}$ (42)

and

$p_{n+1}=p_n+h\left\{f\left(x_n,p_n\right)+A{f}^1\left(x_n,p_n\right)+B{f}^2\left(x_n,p_n\right)\right\}$ (43)

Therefore,

$\begin{array}{c}{y}_{n+1}-p_{n+1}=y_n-p_n+h\{f\left(x_n,y_n\right)-f\left(x_n,p_n\right)+A\left[f^1\left(x_n,y_n\right)-f^1\left(x_n,p_n\right)\right]\\ +B\left[f^2\left(x_n,y_n\right)-f^2\left(x_n,p_n\right)\right]\}\end{array}$ (44)

Applying the mean value theorem as before, we have

$\begin{array}{c}{y}_{n+1}-p_{n+1}=y_n-p_n+h\{\frac{\delta f\left(x_n,p_n\right)}{\delta p_n}\left(x_n-p_n\right)+A\left[\frac{\delta f^1\left(x_n,p_n\right)}{\delta P_n}\left(x_n-p_n\right)\right]\\ +B\left[\frac{\delta f^2\left(x_n,y_n\right)}{\delta p_n}\left(x_n-p_n\right)\right]\}\end{array}$ (45)

$\begin{array}{c}{y}_{n+1}-p_{n+1}=y_n-p_n+h\{{\text{sup}}_{\left(x_n,p_n\right)\in D}\frac{\delta f\left(x_n,p_n\right)}{\delta p_n}\left(x_n-p_n\right)\\ +A\underset{\left(x_n,p_n\right)\in D}{\text{sup}}\frac{\delta f^1\left(x_n,p_n\right)}{\delta p_n}\left(x_n-p_n\right)\\ +B\underset{\left(x_n,p_n\right)\in D}{\text{sup}}\frac{\delta f^2\left(x_n,p_n\right)}{\delta p_n}\left(x_n-p_n\right)\}\end{array}$ (46)

$y_{n+1}-p_{n+1}=y_n-p_n+h\left\{L\left(x_n,p_n\right)+A{L}_1\left(x_n,p_n\right)+B{L}_2\left(x_n,p_n\right)\right\}$ (47)

Taking absolute value of both sides of (47) gives

$\left|y_{n+1}-p_{n+1}\right|\le \left|y_n-p_n\right|+h\left|L+A{L}_1+B{L}_2\right|\left|x_n-p_n\right|$ (48)

Let $N=h\left|L+A{L}_1+B{L}_2\right|$ and $y\left(x_o\right)=\eta$ , $P\left(x_o\right)={\eta }^{*}$ , given $\epsilon >0$ , then

$\left|y_{n+1}-p_{n+1}\right|\le N\left|y_n-p_n\right|$ (49)

and

$\left|y_{n+1}-p_{n+1}\right|\le N\left|\eta -{\eta }^{*}\right|<\epsilon$ , for every $\epsilon >0$ (50)

Then we conclude that our method (19) is stable and hence convergent.

5. The Implementation of the Integrator

Example 1

Using the Integrator (19) to solve the initial value problem

$y^{\prime }=2x^2-y,y\left(0\right)=-1,\text{intheinterval}0\le x\le 1$

The analytical solution $y\left(x\right)=-5{\text{e}}^{-x}+2x^2-4x+4,h=0.1$

Example 2

Consider the initial value problem

$y^{\prime }=2x-y,y\left(0\right)=1,\text{intheinterval}0\le x\le 1$

The analytical solution $y\left(x\right)=3{\text{e}}^{-x}-2\left(x+1\right),h=0.1$

6. Summary and Conclusion

In this paper, we have proposed a new integration for the solution of standard initial value problem of first order ordinary differential equations. The new method was found to be convergence, consistence, and stable.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Ayinde, S.O., et al. (2015) A New Numerical Method for Solving First Order Differential Equations. American Journal of Applied Mathematics and Statistics, USA, 3, 156-160.
[2]	Fatunla, S.O. (1988) Numerical Methods for Initial Value Problems in Ordinary Differential. Academic Press, San Diego, USA.
[3]	Ibijola, E.A., et al. (2010) On a New Numerical Scheme for the Solution of Initial Value Problems in ODE. Australian Journal of Basic and Applied Sciences, 4, 5277-5282.
[4]	Ogunrinde, R.B. (2010) A New Numerical Scheme for the Solution of Initial Value Problem (IVP) in Ordinary Differential Equations. Ph.D. Thesis, Ekiti State University, Ado Ekiti.
[5]	Zwillinger, D. (1997) Handbook of Differential Equations. 3rd Ed., Academic Press, Boston, MA.
[6]	Stoer, J. and Bulirsh, R. (1966) Numerical Treatment of ODEs by Extrapolation Methods. Numerische Mathematical, 8, 93-104.
[7]	Henrici, P. (1962) Discrete Variable Methods in ODEs. New York, John Wiley and Sons., U.S.A.
[8]	Lambert, J.D. (1972) Introductory Mathematics for Scientists and Engineers. John Wiley & Sons, New York.
[9]	Shepley, L.R. (1984) Differential Equations. Third Edition, John Wiley & Sons Inc., Canada.