# An exponential fitting predictor-corrector formula for stiff systems of ordinary differential equations.

1. Introduction

We consider a system of first order initial value problem,

y' = f(x, y), y([x.sub.0]) = [y.sub.0] (1.1)

for y, f [member of] [R.sup.m], x [member of][a, b]

In addition, it is assumed that problem (1.1) satisfies the hypothesis of existence and uniqueness theorem. In this paper we focus our attention on initial value problem whose Lipschtz condition is of order one with respect to y.

i.e., L = [delta]f/[delta]y

Liniger and Willougbby (1967) introduced the concept of exponential fitting and suggested three new A-stable schemes with k = 1. However Cash (1981) derived an exponentially Multiderivative multistep method of orders up to 5 with step that all the step number k = 2 and a numerical investigation of these methods shows that all are A-stable for all choices of fitting parameters. Okunuga (1994, 1997) developed second derivative multistep methods of orders 2, 4, 5, and 6 for stiff problems in which the method is applicable. Otunta and Abhulimen (2005, 2006,) derived exponentially fitted third derivative multistep methods of various order for stiff initial value problems. The numerical result indicated that our methods are A-stable.

Definition 1: The initial value problem (1.1) is said to be stiff over the interval R if for every x[member of] R, the eigenvalues ([[lambda].sub.t](x), t = 1,...,n) of the Jacobian ([delta]f/[delta]y)satisfy the following conations,

(a) Re [[lambda].sub.t](x) < 0, t = 1,...,n (1.2)

(b) max [absolute value of Re [[lambda].sub.r](x)/Re [[lambda].sub.t](x)] >> 1, r, t = 1,...,n (1.3)

where [x.sub.n] = a + nh, n = 0,1,...N (1.4)

represent a uniform subdivision of the interval of integration R with the step length h given by h = (a - b)/N.

2. The Development of the Integration Formula

The development of our proposed third derivative integration formula involves a pair of formulae as follows.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.1)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.2)

Where [f.sub.n+i], [g.sub.n+i] and [v.sub.n+i] are respectively first, second and third derivative functions of y(x) evaluated at step n + i. thus equations (2.1) and (2.2) serves as the predictor and corrector respectively. When deriving exponentially fitted integrator, the approach is to allow both (2.1) and (2.2) to posses free parameter

For the purpose of efficient implementation, we cast the integration formula into predictor-corrector process.

Thus, the derivation of predictor-corrector integrations formula of order 5 involves two stages. Using Taylor series expansion,. we obtains five set of equations with 12 unknown parameters from equation (2.1.)

To obtain the predictor formula, we set [[bar.[alpha]].sub.1] = [[bar.[beta]].sub.1] = [[bar.[phi]].sub.1] = [[bar.[omega]].sub.1] = [[bar.[omega]].sub.2] = 0 and [[bar.[beta]].sub.2] = a, as free parameter, and [[alpha].sub.2] = + 1

we then obtain the following set of simultaneous equations from (2.115)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.3)

when we solve the above equations we obtain,

[[bar.[alpha]].sub.0] = -1 [[bar.[alpha]].sub.2] = 1, [[bar.[beta]].sub.0] = 2 - a, [[bar.[phi]].sub.0] = [5/3] - [4/3] a, [[bar.[phi]].sub.2] = [1/3] - [2/3] a, [[bar.[omega]].sub.0] = [2/3] - [2/3] a

when these values of the parameters are substituted into (2.10), we obtain the predictor formula as:

[[bar.y].sub.n+2] - [[bar.y].sub.n] = h[[ay'.sub.n+2] + (2 - a)[y'.sub.n]] + [h.sup.2] [([1/3] - [2/3] a)[y".sub.n+2] + (([5/3] - [4/3]a)[y".sub.n]] + [h.sup.3][(([2/3] - [2/3] a)[y".sub.n]] (2.4)

Now, for exponential fitting purpose, we apply (2.4) to scalar test problem to obtain,

y' = [lambda]y, y(0) = 1, Re([lambda]) < 0 and [lambda]h = u (2.5)

[y.sub.n+2][1 - au -([1/3] - [2/3]a)[u.sup.2]]= [y.sub.n][1 + (2 - a)u + ([5/3] - [4/3] a)[u.sup.2]([2/3] - [2/3] a)[u.sup.3]] (2.6)

[[[bar.y].sub.n+2]/[y.sub.n]] = 1 + (2 - a)u + ([5/3] - [4/3]a)[u.sup.2] + ([2/3] - [2/3]a)[u.sup.3]/1 - au - ([1/3] -[2/3]a)[u.sup.2] = R([bar.u]) (2.7)

For the purpose of stability analysis, we obtain the free parameter a from (2.7).

But [y.sub.n+2]/[y.sub.n] = [e.sup.2q] then equation (2.7) yields,

1 + 2u + [1/3][u.sup.2](5 + [e.sup.2q]) + [2/3][u.sup.3] - [e.sup.2q] = a[[2/3][u.sup.2]([e.sup.2q] + 2) + u(1 + [e.sup.2q]) + [2/3][u.sup.3]] (2.8)

From which we obtain,

a = 1 + 2u + [1/3][u.sup.2](5 + [e.sup.2q]) + [2/3][u.sup.3] - [e.sup.2q]/[2/3][u.sup.2](2 + [e.sup.2q]) + u(1 + [e.sup.2q]) + [2/3][u.sup.3] (2.9)

Again to obtain the corresponding order 5 corrector integration formula, we obtain six set of equations from (2.2) by using Taylor series expansion as it was done in the predictor.

We therefore impose the same condition as in predictor, and in addition, we let [[beta].sub.3] = b as free parameters and [[alpha].sub.2] = 1, the values of the unknown parameters are obtained as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.10)

When we solve the above equations we obtain

[[alpha].sub.0] = -1, [[alpha].sub.2] = 1, [[beta].sub.0] = [6/5] - [43/16]b, [[beta].sub.2] = [4/5] + [27/16]b, [[beta].sub.3] = b

[[phi].sub.0] = [8/5] + [33/16]b, [[phi].sub.2] = -[1/5] - [27/8]b, [[omega].sub.0] = [2/15] - [9/8]b.

When these values are substitute into (2.2), we obtain fifth order corrector formula as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(2.11) As before, we apply (2.11) to test function (2.5) to obtain

[[y.sub.n+2]/[y.sub.n]] = 1 + ([6/5] - [43/16]b)u + bu. [[y.sub.n+3]/[y.sup.n]] + ([8/5] - [33/16]b)[u.sup.2] + ([2/15] - [9/8]b)[u.sup.3]/1 - ([4/5] - [27/16]u + ([1/5] + [27/8]b)[u.sup.2] = R(u) (2.12)

We need to obtain the ratio [y.sub.n+3]/[y.sub.n] in (2.11) for the purpose of exponential fitting condition and uniting both the predictor and corrector formulas to form a single integrator formula.

But, [[y.sub.n+3]/[y.sub.n]] = [([e.sup.2u]).sup.3/2] = [[[[bar.y].sub.n+2]/[[bar.y].sub.n]].sup.3/2] = [[R([bar.u])].sup.3/2] = R * (u) (2.13)

Using (2.13), equation (2.12) becomes,

[[y.sub.n+2]/[y.sub.n]] = [1 + ([6/5] - [43/16]b)u + bu. R * ([bar.u]) + ([8/5] + [33/16]b)[u.sup.2] + ([2/15] - [9/8]b)[u.sup.3]/1 - ([4/5] + [27/16]b)u + ([1/5] + [27/8]b)[u.sup.2])] = R(u) (2.14)

Equation (2.14) now unites both the predictor and corrector formulas. So, equation (2.14) is now the exponentially third derivative integration formula of order 5. We then obtain the value of the free parameter b from (2.12) for stability purpose.

b = 1 + [1/5](6u - 8[u.sup.2]) + [2/15][u.sup.3] + [4/5]u[e.sup.2u] - [e.sup.2u](1 + [1/5][u.sup.2])/[27/8][u.sup.2][e.sup.2u] - [33/16][u.sup.2] + [9/8][u.sup.3] - [1/16](27u [e.sup.2u] - 43u) - u[e.sup.3u] (2.16)

3.0 Stability Consideration of the Method

In order to determine the stability of the method, the evaluation of the values of the free parameters a and b in the open left half plane (-[infinity], 0] become very necessary.

In Cash (1981), it is straightforward to find the conditions which a and b need to satisfy such that [absolute value of [y.sub.n+2]/[y.sub.n]] < 1 i.e., [absolute value of R(u)] < 1 for all u with R(u) < 0. Thus necessary and efficient for this inequality to hold are given by application of the maximum modulus theorem and are; (i) [absolute value of R(u)] < 1 (ii) R(u) analytic in R(u)< 0.

If condition (i) holds, it follows that R(u) is analytic as u [right arrow] - [infinity]. Thus, (i) and (ii) will guarantee A stability.

However, since equation (2.7) is contained in (2.14), then equation (2.14), must also satisfy [absolute value of R(u)] < 1, i.e., in equation (2.14), R(u) - 1 < 0 as u [right arrow] -[infinity].

Now from equation (2.14), we have,

[1 + ([6/5] - [43/16]b)u + bu [R.sup.*](u) + ([8/5] + [33/16]b)[u.sup.2] + ([2/15] - [9/8])[u.sub.3]/1 - ([4/5] + [27/16]b)u + ([1/5] + [27/8]b)[u.sup.2]] -1 < 0 (31)

so, as u [right arrow] -[infinity], equation (3.1) yields,

[[b/q] [R.sup.*]([bar.u]) + [89/16]b + [7/5] + ([2/15] - [9/8b])u/[1/5] + [27/8]b] <0 (3.2)

when we simplify equation (2.17) further, we obtain,

b > [-8/135] i.e b < 0. (3.3)

Furthermore, we also examine the stability function of the predictor, if it satisfies the inequality [bar.R](u) - 1 < 0.

so, from (2.7), we have,

[1 + (2 - a)u + ([5/3] - [4/3]a)[u.sup.2] +([2/3] - [2/3]a)[u.sup.3]/1 - au - ([1/3] - [2/3]a)[u.sup.2]] -1 < 0 (3.4)

So, we then have from (3.4);

[2 - 2a/[1/3] - [2/3]a] < 0 (35)

this gives, a < [1/2] or a > 1 (.3.6)

Thus, the inequality (3.3) and (3.6) will help to determine the interval of stability of this method.

Again if we take the case of R(u) > -1, we have a < 1/2 or a > 1.

In order to determine the interval of absolute stability of the method, we find finite limits of both a(u) and b(u) as u [right arrow] 0. We take the limits using L'Hospital rule. From

(2.9), we have; [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Similarly, from (2.16) we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Thus, as u decreases the values of a and b are monotonically increasing. However, we established from the algebraic work done on the condition that, R(u) is satisfy with the ranges of a and b, that within the ranges a[member of](1/2, 4/5) and b [member of]([8/135], -64/105), ur integration formula satisfy A-stability condition given by the maximum modulus.

4.0 Numerical Experiments

To show the effectiveness and validity of our newly derived integration formula, we present some numerical examples below. All numerical examples are coded in Fortran 77 and implemented on digital computer..

Problem 1 Jackson-Kenue (1974)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.1)

The eigenvalues of the Jacobian matrix of the system are [[lambda].sub.1] = -2 and [[lambda].sub.1] = -96 with stiffness ratio 48.

The exact solution is given as: y = (95[e-.sup.2x] -48[e.sup.-96x])/47 z = (48[e.sup.-96x] - [e.sup.-2x])/47

Now, for the purpose of comparison, we denote our newly third derivative exponential fitted predictor-corrector method of order five as AF5,

Jackson and Kenue (1974) as J-K, Okunuga (1994) method of order 6 as OK6, Abhulimen and Otunta (2007, 2008) methods of order 7 and 8 as AB7 and AB8 and Abhulimen (2009) method of order 9 as NM9 respectively.

From table (4.1) above at h = 0.03125 and 0.0625, we observed that our proposed method performed better than other existing methods. The level of accuracy of our method is further demonstrated by representing each of the step lengths in form of a chart.

[FIGURE 4.1 OMITTED]

Problem 2 (Enright 1972).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The error tolerance given by Enright (1972) is [10.sup.-6].

The eigenvalues of the system are the non-zero elements in the leading diagonal. So, the eigenvalues are [[lambda].sub.1] = -0.1, [[lambda].sub.2] = -10, while the exact solution is given as,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The stiffness ratio is given as [10.sup.2]. While the values of [A.sub.i] and [B.sub.i] are determined from the initial condition imposed on the derivations of y (x).

For the purpose of comparison of results, we denote the error of our fifth and eighth order integrators respectively by ER5 and ER8.

The numerical results of using the newly derived methods of orders 5 to AB8 on problem 2, are given in tables 4.2(a) in the range of x [member of] [0,10].

The error tolerance of the numerical solution of problem 2, given by Enright (1972) is [10.sup.-6].

However, from the results above, we observed that the errors in our newly derived methods of order 4 and 6, could be raised [10.sup.-15].

Problem 3 (Second order differential equation (Okunuga (1994))

[[d.sup.2]y/dx] + 1001 [dy/dx] + 1000y = 0,] y(0) = 1, [y.sup.1](0) = 1 (4.3)

The system (4.3) can be rewritten as a first order system.

[d[y.sub.1]/dx] = [y.sub.2] [y.sub.1](0) = 1 [d[y.sub.2]/dx] = -1001[y.sub.2] - 1000[y.sub.1] [y.sub.2] (0) = -1 (4.4)

Thus we obtain a 2x2 system of stiff IVP. The eigenvalues of the Jacobian matrix of (4.4) are [[lambda].sub.1] = -1 and [[lambda].sub.2] = -1000.

The eigen value is 1000

The general solution of (4.3) is y(x) = A[e.sup.-x] + B[e.sup.-1000x]. If we impose the initial conditions in 0 [less than or equal to] x [less than or equal to] 1 the exact solution is y(x) = [e.sup.-x].

Note: for purpose of comparison, A-K5 represents Abhulimen and Okunuga (2008).

The result of this problem using the newly derived method obtained at x = 1; is given in table (4.3) below. Note that A-K5 represents Abhulimen and Okunuga (2008).

It will be observed from table (4.3), that our newly derived method of order 5 has higher degree of accuracy than OK6 and A-K5. To illustrate the efficiency pf our method of order 5 over existing methods, we represent the Numerical results of problem 3 in form of charts, as shown in figures (44) and (4.5) below.

5.0 Conclusion

The numerical results so far obtained in this paper, show the efficiency of the newly derived integrator of order five. We also observed that for an exponentially fitted problems, our integrator do not use small step lengths, as may be required by many multistep methods before good accuracy is obtain. Finally, the new integrator derived in this paper is capable of handling stiff problems for which exponential fitting is applicable.

References

[1] C.E. Abhulimen (2006): 'On Exponentially Fitted Multiderivative Numerical Methods for Stiff Initial Value Problems in Ordinary Differential Equations'. Ph.D thesis, Department of Mathematics, University of Benin, Benin City, Nigeria.

[2] C.E. Abhulimen and F.O. Otunta (2006), "A Sixth Order Multiderivative Multistep Methods for Stiff System of Differential Equations". International Journal of Numerical Mathematics (IJNM) 2 No 1, 248-268.

[3] C.E Abhulimen and F.O Otunta (2007): A Family of Two Step Exponentially Fitted Multiderivative Methods for the Numerical Integration of Stiff IVPs in ODEs. International Journal of Numerical Mathematics IJNM 2 No. 1, 1-21.

[4] C.E Abhulimen and Okunuga S.A (2008) "Exponentially-Fitted Second Derivative Multistep Methods for Stiff Initial Value Problems in Ordinary Differential Equations. Journal of Engineering Science And Applications (JESA) Vol 5, Number 1&2, pp36-49

[5] C.E Abhulimen (2008) " A Class of Methods of Exponentially-Fitted Third Derivative Multistep Method for Solving Stiff Differential Equations." International Journal of Physical Sciences (IJPS) Vol 3 (8) pp 188-193. SSN 1992-1950@2008 Academic Journal.

[6] C.E. Abhulimen and Otunta F.O (2009): A new class of Exponential Fitting Numerical Integration for Initial Value Problems in Ordinary Differential Equations. Nigerian Mathematical society (NMS), vol 28, pp1-14.

[7] J.R. Cash (1981): "On Exponential Fitting of Composite Multiderivative Linear Methods", SIAM J. Multistep Numerical Anal. 18, No 5, 808-821

[8] S.C Chu and .M. Berman (1974): "An Exponential Method for the Solution of Systems of Ordinary Differential Equations". Comput. Mach 117 (12), 699-702.

[9] W.H Enright and J.D. Pryce (1983): "Two FORTRAN Packages for Assessing IV Methods". Technical Report No. 16/83. Department of Computer Science, University of Toronto, Canada

[10] Jackson L.W and Kenue S.K (1974): "A Fourth Oredr Exponential Fitted Method". SIAM Journal on Numerical Anal, Vol 6. No 4, pp. 965-978.

[11] W.S. Liniger & Willoughby (1970): "Efficient Integration Methods for Stiff System of ODEs". SIAM J. on Numer Anal. 7, 47-65.

[12] S.A. Okunuga (1994):"Composite Multiderivative Linear Multistep Methods for Stiff IVPS in ODES". Ph.D Thesis, Dept. of Mathematics University of Lagos, Nigeria.

[13] S.A. Okunuga (1997): "Fourth Order Composite Two-step Method for Stiff Problems". Inter J.Comp Maths 2, 39-47.

[14] F.O Otunta and C.E Abhulimen (2005): "A Fourth Order Exponentially Fitted Multiderivative Method for Stiff Initial value Problems". Nigerian Association of Mathematical Physics. 9, 295-306.

[15] D. Voss (1988): "A Fifth Order Exponentially Fitted Formula". SIAM J. Numer. Anal. 25, No 3.

C.E Abhulimen

Department of Mathematics, Ambrose Alli University, Ekpoma, Edo State, Nigeria

Email: cletusabhulimen@yahoo.co.uk
```Table 4.1: Comparison of Numerical results of Problem 1 at x = 1.0
for various step lengths.

Step Methods of y(1) z(1) x
length integration [10.sup.-2]

0.03125 AF5 0.27354958 -0.28794694
AB8 0.27356554 -0.28794325
NM9 0.273565455 -0.28794375
0.0625 AF5 0.23554556 -0.28794741
AB7 0.273545505 -0.28794751
AB8 0.273505505 -0.28774751
NM9 0.273505505 -0.287747
OK6 0.27354657 -0.28355004
J-K 0.27355005 -0.28794742
0.05 AF5 0.27354738 -0.28794461
OK6 0.27354864 -0.28794594
Exact solution 0.27355004 -0.28794110

Step Methods of Minimum Maximum
length integration error (y) error (z)

0.03125 AF5 4.5 x [10.sup.-7] 4.7 x [10.sup.-9]
AB8 3.7 x [10.sup.-3] 4.0 x [10.sup.-5]
NM9 3.7 x [10.sup.-3] 4.0 x [10.sup.-5]
0.0625 AF5 5.9 x [10.sup.-6] 6.2 x [10.sup.-8]
AB7 4.0 x [10.sup.-5] 6.0 x [10.sup.-5]
AB8 1.7 x [10.sup.-5] 1.8 x [10.sup.-5]
NM9 7.9 x [10.sup.-75] 8.3 x [10.sup.-7]
OK6 3.4 x [10.sup.-6] 3.7 x [10.sup.-8]
J-K 5 x [10.sup.-7] 4 x [10.sup.-7]
0.05 AF5 2.7 x [10.sup.-6] 1.4 x [10.sup.-8]
OK6 1.3 x [10.sup.-6] 1.4 x [10.sup.-8]
Exact solution

Table 4.2(a): Performance of the new integrator on problem
2, h = 0.5, x [member of] [0, 10].

Method [y.sub.1](10) [y.sub.2](10)

AF5 0.3678794352 0.77472 x 10
AB8 0.367879207 0.64136 x 10
Exact Solution 0.3678794357 0.3701 x 10

Errors

ER5 5.0 x [10.sup.-7] 1.3 x [10.sup.-15]
ER8 2.3 x [10.sup.-7] -6.4 x [10.sup.-11]

Method [y.sub.3](10) [y.sub.4](10)
-3.34435

AF5 -3.34475 x 10 -36.7879432
AB8 -3.34435 x 10 -36.7879207
Exact Solution -3.344358 -36.78794357
ER5 -4.4 x [10.sup.-7] -5.0 x [10.sup.-6]
ER8 -2.1 x [10.sup.-6] -2.3 x [10.sup.-5]

Table (4.3): Numerical results on second order ODE.

Step h Method y(1) Error (y)

0.1 AF5 0.367879435 5.2 x [10.sup.-8]
OK6 0.367879436 5.6 x [10.sup.-8]
A-K5 0.36787960 1.8 x [10.sup.-8]
0.125 AF5 0.367879441 2.7 x [10.sup.-8]
A-K5 0.36789500 1.4 x [10.sup.-8]
Exact Solution 0.367879435

Figure (4.2)

AF5 AB7 AB8

Series l 5.90E-10 4.00E-05 1.70E-05
Series 2 6.20E-08 6.00E-05 1.80E-05

NM9 OK6 J-K

Series l 7.90E-05 3.40E-06 5.00E-07
Series 2 8.30E-07 3.70E-08 4.00E-07

Figure (4.3)

1 2

AF5 2.70E-06 2.70E-08
OK6 1.30E-06 1.40E-08

Figure (4.4)

Absolute eeors of numerical values of y(x)
of methods of order 5 (0.1) and (0.125)
for problem 3 at h = 0.1

methods

AF5 OK6 A-K5

Series 1 5.20E-09 5.60E-08 1.80E-07

Figure (4.5)

for h = 0.125

methods

AF5 A-K5

Series 1 2.70E-08 1.40E-06
```
COPYRIGHT 2009 Research India Publications
No portion of this article can be reproduced without the express written permission from the copyright holder.