1. Introduction

In this paper, we consider the numerical solution of stiff and oscillatory first-order differential equations of the form,

(1)

where is the initial point, is the solution at the initial point and is assumed to satisfy Lipchitz condition stated below.

Theorem 1[13]: Let be defined and continuous for all points in the region defined by, finite, and let there exist a constant such that, for every such that are both in;

Then, if is any given number, there exists a unique solution of the initial value problem (1), where is continuous and differentiable for all .

Proof, see[11]

According to[18], equation (1) is used in simulating the growth of population, trajectory of particles, simple harmonic motion, deflection of a beam, etc. Few equations that are modeled in higher order differential equations are first reduced to systems of first-order before appropriate method of solution is applied. Most often, these problems do not have a closed form solutions, hence appropriate methods are adopted to solve such problems. Different methods have been proposed ranging from predictor-corrector methods to block methods. Despite the success recorded by the predictor-corrector methods, its major setback is that the predictors are in reducing order of accuracy especially when the value of the step-length is high and moreover the results are at overlapping interval[3]. Block methods which have advantage of being more efficient in terms of cost implementation, time of execution and accuracy were developed to handle the setbacks of predictor-corrector methods[19],[2] and[20].

Definition 1[13]: A differential equation is said to be stiff if , where is the eigen value of the differential equation.

Definition 2[6]: A nontrivial solution (function) of an ODE is called oscillating if it does not tend either to a finite limit or to infinity (i.e. if it has an infinite number of roots). The differential equation is called oscillating, if it has at least one oscillating solution.

In the quest for a method that gives better stability condition, [14] proposed an approximate solution which combined power series and exponential function. In this paper, we extend their work by developing a block integrator with step number k=5.

2. Derivation of the Extended Block Integrator

We consider an approximate solution that combines power series and exponential function of the form,

(2)

Interpolation and collocation procedures are used by choosing interpolation point at a grid point and collocation points at all points giving rise to system of equations whose coefficients are determined by using appropriate procedures. The first derivative of (2) is given by,

(3)

where for and is continuously differentiable. Let the solution of (1) be sought on the partition of the integration interval with a constant step-size , given by, , .

Then, substituting (3) in (1) gives,

(4)

Now, interpolating (2) at point and collocating (4) at points , leads to the following system of equations,

(5)

where

and

Solving (5), for and substituting back into (2) gives a continuous linear multistep method of the form,

(6)

where the coefficients of and are given by,

(7)

where . Evaluating (6) at gives a block scheme of the form,

(8)

where

,

,

,

,

3. Analysis of Basic Properties of the Extended Block Integrator

3.1. Order of the Extended Block Integrator

Let the linear operator associated with the block (8) be defined as,

(9)

Expanding (9) using Taylor series and comparing the coefficients of gives,

(10)

Definition 3[10]: The linear operator and the associated continuous linear multistep method (6) are said to be of order if

is called the error constant and the local truncation error is given by,

(11)

For extended block integrator,

(12)

Expanding (12) in Taylor series gives,

(13)

Equating the coefficients of the Taylor series expansion to zero yields, Therefore, the extended block integrator is of order six.

3.2. Zero Stability

Definition 4[10]: The block integrator (8) is said to be zero-stable, if the roots of the first characteristic polynomial defined by satisfies and every root satisfying have multiplicity not exceeding the order of the differential equation. Moreover, as where is the order of the differential equation, is the order of the matrices , see[5] for details.

For our integrator,

(14)

. Hence, the extended block integrator is zero-stable.

3.3. Consistency

The extended block integrator (8) is consistent since it has order .

3.4. Convergence

The extended block integrator is convergent by consequence of Dahlquist theorem stated below.

Theorem 2[8]: The necessary and sufficient conditions that a continuous LMM be convergent are that it be consistent and zero-stable.

3.5. Region of Absolute Stability

Definition 5[22]: Region of absolute stability is a region in the complex plane, where . It is defined as those values of such that the numerical solutions of satisfy for any initial condition.

We shall adopt the boundary locus method to determine the region of absolute stability of the extended block integrator. This is achieved by substituting the test equation,

(15)

into the block formula gives (8). This gives,

(16)

Figure 1. Showing the Stability Region of the Extended Block Integrator

Thus,

(17)

since is given by . Equation (17) is our characteristic/stability polynomial. For the extended block integrator, equation (17) is given by,

(18)

This gives the stability region shown in the figure below.

According to Fatunla (1988), stiff algorithms have unbounded RAS. Thus, from figures 1 above, the extended block integrator is suitable for solving stiff problems. Also, Lambert (1973) said that the stability region for L-stable schemes must encroach into the positive half of the complex plane. Thus, the extended block integrator is L-stable.

4. Numerical Experiments

We shall evaluate the performance of the block integrator on some challenging stiff and oscillatory problems which have appeared in literature and compare the results with solutions from some methods of similar derivation. The numerical results are obtained using MATLAB.

The following notations shall be used in the tables below;

ERR- Exact Solution-Computed Solution

ESSI- Error in[17]

ESYO- Error in[20]

Problem 1

Consider the highly stiff ODE

(19)

which has the exact solution

(20)

where is a positive constant and is a smooth slowly varying function. Equation (20) exhibits two widely different time scales: a rapidly changing term associated with and a slowly varying term associated with .

In[17], the authors considered a special case of (19) where . They solved the problem 1 by adopting an L-stable hybrid block Simpson’s method of order six. The authors in[20] also solved problem 1 using a block integrator with step number .

Problem 2

Consider the highly oscillatory ODE

(21)

whose exact solution is,

(22)

Though in[22], the authors did not solve this problem, he however observed that it has a solution that oscillates and grows exponentially in . He further stated that most numerical methods do not perform well on this problem. The authors in[20] solved problem 2 by adopting a block integrator with step number .

Table 1. Showing the results for stiff problem 1

Table 2. Showing the result for oscillatory problem 2

5. Conclusions

In this paper, we have presented an extended block integrator for the solution of stiff and oscillatory first-order ordinary differential equations. Our aim was to construct highly stable block integrator which is computationally more efficient than many of the existing numerical integrators for stiff and oscillatory problems. The approximate solution (basis function) adopted in this paper produced a block integrator with L-stable stability region. This made it possible for the block integrator to perform well on stiff and oscillatory problems. The extended block integrator proposed was also found to be zero-stable, consistent and convergent. The block integrator was also found to perform better than some existing methods in view of the numerical results presented.