ASYMPTOTIC-NUMERICAL METHOD FOR SINGULARLY PERTURBED DIFFERENTIAL DIFFERENCE EQUATIONS OF MIXED-TYPE

Abstract

A computational method for solving singularly perturbed boundary value problem of differential equation with shift arguments of mixed type is presented. When shift arguments are sufficiently small (o(ε)), most of the existing method in the literature used Taylor's expansion to approximate the shift term. This procedure may lead to a bad approximation when the delay argument is of O(ε). The main idea for this work is to deal with constant shift arguments, which are independent of ε. In the present method, we construct the formally asymptotic solution of the problem using the method of composite expansion. The reduced problem is solved numerically by using operator compact implicit method, and the second problem is solved analytically. Error estimate is derived by using the maximum norm. Numerical examples are provided to support the theoretical results and to show the efficiency of the proposed method.

1. Introduction

Many phenomena in real life and science may be modeled mathematically by delay differential or differential difference equations (DDEs). Equations of this type arise widely in scientific fields such as biology, medicine, ecology and physics, in which the time evolution depends not only on present states but also on states at or near a given time in the past [3]. If we restrict the class of DDEs to a class in which the highest derivative is multiplied by a small parameter, then we get a class of singularly perturbed differential difference equations (SPDDEs). These equations are used to model a large variety of practical phenomena, for instance, variational problems in control theory [4], description of the so-called human pupil-light reflex [10], evolutionary biology [19] and a variety of model for physiological processes or diseases [11].

In this paper, we consider the following SPDDE with negative as well as positive shifts on the interval Ω = (0, ℓ):

where 0 < ε ≪ 1 is the singular perturbation parameter and δ(ℓ < 2δ) is the shift parameter. We assume that the functions a, b, c, d, f, ϕ and ψ are sufficiently smooth such that

Under the above assumptions, the solution u of the problem (1) exhibits a single boundary layer at the left end of the interval = [0, ℓ]. When a(x) < 0, the solution exhibits boundary layer at the right end of the interval . If a(x) = 0, then the solution can exhibit layer or oscillatory behavior depending on the sign of reaction coefficient. Here, we deal with the two cases a(x) > 0 and a(x) < 0.

It is well known that the standard discretization methods for solving singular perturbation problems are not useful and fail to give accurate results when the perturbation parameter ε tends to zero. This motivates the need for other methods to solve this type of equations, whose accuracy does not depend on the perturbation parameter. For more details, one may refer to [2,12].

The study of boundary value problems for such equations was initiated by Lange and Miura [8,9], where the authors provided asymptotic approach to SPDDEs with small shifts in the case of layer behavior and rapid oscillations. In recent years, there has been growing interest to develop numerical methods for SPDDEs. In this area, we mention [5,6,7,13,17]. Most of these works focused mainly on the problems with very small shifts (of the order of ε). Kadalbajoo and Sharma [6] considered the problem where delay and advance parameter are of o(ε) and a convection coefficient is absent. Such problems have oscillatory behavior or boundary layers at both ends of the interval. Kadalbajoo et al. [5], Mohapatra and Natesan [13] and Salama and Al-Amery [17] presented numerical methods to solve SPDDEs with small shifts of mixed type when convection coefficient is non-vanishing. Recently, Subburayan and Ramanujam [18] have applied the initial value technique for singularly perturbed convection-diffusion problems with negative shift.

In this paper, we construct an asymptotic-numerical method to solve SPDDEs with shifts of mixed type, these shifts are fixed and do not depend on ε. In this method, an asymptotic expansion of the solution of this problem is constructed using the basic ideas of the method of matched asymptotic expansion [14,15]. Then, the initial-value problem (the reduced equation) is approximated by operator compact implicit (OCI) method [16], and the second problem, which retains the order of the original problem, is solved analytically. We show that the present method is useful for obtaining the numerical solution of the considered problem in both cases when boundary layer is at the left end as well as at the right end of the interval.

The remainder of the paper is organized as: In Section 2, a maximum principle and some important properties of the exact solution and its derivatives are established. An asymptotic expansion approximation is constructed in Section 3. The proposed numerical method is described in Section 4. Error estimate is derived in Section 5. In Section 6, numerical examples are presented, which validate the theoretical results. Finally, conclusion and discussion are indicated in Section 7.

 

2. The continuous problem

The problem (1) is equivalent to

where Ω1 = (0, ℓ−δ], Ω2 = [ℓ−δ, δ] and Ω3 = [δ, ℓ). Throughout this paper, C is a generic positive constant independent of ε and discretization parameter h, and we use the simple notation for the discrete maximum norm

Lemma 2.1 (Continuous minimum principle). Let Ψ(x) be a smooth function satisfying Ψ(0) ≥ 0 and Ψ(ℓ) ≥ 0. Then LεΨ(x) ≼ 0, ∀x ∈ Ω implies that Ψ(x) ≥ 0, ∀x ∈ .

Proof. Let x⋆ ∈ be such that Ψ(x⋆) =min{Ψ(x), x ∈ } and Ψ(x⋆) < 0. Clearly x⋆ ≠ 0, x⋆ ≠ ℓ and also Ψ′(x⋆) = 0 and Ψ′′(x⋆) ≥ 0.

We have the following:

which contradicts the hypothesis that LεΨ(x⋆) ≼ 0. Therefore Ψ(x⋆) ≥ 0. But since x⋆ was arbitrary point in , so that Ψ(x) ≥ 0, ∀x ∈ . □

Lemma 2.2. Let u be the solution of the problem (1), and let b(x)−c(x)−d(x) ≥ λ > 0, x ∈ Ω. Then

where C0 is a positive constant.

Proof. Consider the barrier function Ψ± as

Then, application of the minimum principle to the above barrier function, one can obtain the required result. □

Remark 2.1. Lemma 2.1 implies that the solution of the problem (1) is unique, and the existence of the solution is implied by its uniqueness and the linearity of the considered problem. Further, a bound on the solution is given in Lemma 2.2.

Theorem 2.3. If a, b, c, d and f are sufficiently smooth functions and a(x) ≥ α > 0, ∀x ∈ Ω. Then the derivatives of the solution u of the problem (1) satisfy

Proof. We rewrite (1) in the form

where

Multiplying both sides of (3) by and integrating over (0, x) we get

Integrating (4) from 0 to x, we have

Interchanging the order of the double integral, we get

Using the condition u(ℓ) = ψ(ℓ) in equation (6), we get

Since

and a bound on F(x) is obtained by using the bound on u given in Lemma 2.2.

Hence, we have

Now, from equation (7) we have

Using the above inequality in equation (4), we can conclude that

which proves (2) when k = 1. The estimates for higher order derivatives can be proved by using induction process. □

 

3. An asymptotic expansion

Since the differential difference equations have special properties, we construct piecewise formal solution in Ω1,Ω2 and Ω3. Moreover, we assume that the solution u satisfies the continuity conditions

where u(z−) and u(z+) denote the left and the right limits of u at x = z, respectively.

We represent the solution u(x) in the form

where . The coefficients ui and vi are called outer and boundary-layer correction coefficients, respectively. We refer the reader to [14,15] for complete description. From the boundary conditions, we obtain

Here, we consider the zero order asymptotic expansion approximation to the solution of the original problem (1) which is given by

where u0(x) is the solution of the following reduced problem

and v0(x) is the solution of the boundary value problem

which can be solved analytically as

Remark 3.1. In the case of right boundary layer, the stretching variables are used in the intervals [0, ℓ − δ], [ℓ − δ, δ] and [δ, ℓ], respectively, and one can follow the same procedure used for the case of left boundary layer.

 

4. Numerical method

In this section, we derive operator compact implicit (OCI) method for the initial value problem (11). Hence the numerical solution of the original problem (1) can be obtained by combining this solution with the solution (13).

4.1. Operator compact implicit method. In the following, we derive OCI method of order four for solving the initial value problem (11). To achieve that, we rewrite it in the form

where

We define the present methods in the form

where the operators R and Q are given by

It is worthwhile to mention that U0 and LNU0 are approximations of u0 and Lu0, respectively. The coefficients are to be functions of and h as we show later. Note that one of the coefficients in (16) must be taken as a multiplicative normalizing factor, so that we consider the following condition

For simplicity of notation we use and we also omit the subscript j in the quantities when convenient.

The local truncation error of (15) at is defined in the form

For u0 sufficiently smooth and using Taylor expansion, can be written in the form

where

The sufficient condition for the order of the local truncation error is given in the following lemma.

Lemma 4.1. For the implicit method (15) and (16), a sufficient condition for for k = 0, 1, 2, 3. Furthermore, if (15) and (16) are normalized according to (17). Then has formal order no larger than 4.

Proof. The proof of the this lemma is similar to that of Theorem 1.1 in [16] . □

Here, we construct the method (15) and (16) by the following conditions

Using (19) and the above conditions, we obtain

Now, we define as polynomials in h at each mesh point , j = 1, 2,…,N, and for simplifying the notations, the index j in will be dropped.

Substituting (23) into (22) and imposing (21), we get

Then, the method is defined in the form

where

4.2. A numerical solution for problem (11). In order to obtain the numerical solution for the problem (11), we apply the OCI method (15) and (16) on a specially designed mesh. The presence of shift parameter makes the problem (11) difficult. To overcome this difficulty, the mesh is designed in such a way that the terms containing shifts lie at the mesh points after the discretization. In what follows, we consider the following uniform mesh on Ω:

and we suppose that , where N0 is a positive integer, i.e., xN0 = δ. Thus, the OCI method (15) and (16) for the initial value problem (11) is defined as follows

where

and

β(n) is the value of (26) for U0N0 . The above process is to be repeated with suitable initial value β(0) until the profiles stabilize in both regions. For computational purposes, the iterative process stops at the nth iteration if any one of the following conditions is satisfied:

where ξ is a given tolerance.

4.3. A numerical solution for problem (1). A numerical solution for the original problem (1) is given by

 

5. Error estimate

Lemma 5.1. Let a, b, c and f be sufficiently smooth functions. Then the zeroorder asymptotic approximation ũ given in (10) of the solution u for the problem

(1) satisfies the following

where C is a positive constant independent of ε.

Proof. Following the method of proof that is done in [2] on the intervals [0, ℓ − δ], [ℓ − δ, δ] and [δ, ℓ] yields the desired estimate.

Denote by σ(x) the error of the piecewise cubic interpolation for u(x) defined as

Define

Theorem 5.2. Let a, b, c, d and f be sufficiently smooth functions and satisfy the condition

Then the numerical solution U0 given in (26) to the solution u0 of the problem (11) exists and the following error estimate holds:

In particular, the method (26)-(28) is convergent and the order of convergence is four.

Proof. Let e(xj) = u0(xj) − U0j . Then, for 0 ≼ j ≼ N − N0 and using (26), we get

It is clear that, if = xj + δ ≥ ℓ, then . Otherwise, ∈ (xi−1, xi] for some N0 + 1 ≼ i ≼ N. Thus, = xi − θh, 0 ≼ θ < 1. Using the cubic interpolation, we get

Thus

Hence

Since . Therefore, we get

Similarly, we have

From the condition (31) and using (33) and (34), we have

To obtain a bound on e(x), we observe that

Let , then the above inequality together with (35) complete the proof. □

Remark 5.1. Note that, in the error estimate (32), ∥τ∥ = O(h4) and ∥σ∥ = O(h4). Therefore, it follows that

Theorem 5.3. If u is the solution of the original problem (1), and U is the numerical solution given in (29). Then we have

where C is a positive constant independent of ε and h.

Proof. Let ũ = u0 + v0 be the asymptotic approximation to the solution of problem (1). Then using Lemma 5.1, we have

On the other hand, for 0 ≼ j ≼ N − N0, we have

where U0j is the numerical solution of problem (11) obtained by the OCI method (26)-(28), and V0j defined as

Hence, using triangle inequality, we get

Similarly, one can conclude that

Combining (38) with (39) and (40), we obtaine

□

 

6. Numerical results

In this section, we present several numerical examples to illustrate the applicability and efficiency of the proposed method, we consider the boundary value problem of singularly perturbed differential difference equations with the left end and the right end boundary layer.

For given value of N, the maximum absolute error is calculated as follows

where u is the exact solution and U is the numerical solution. In the case of the exact solution is not known, the maximum absolute error and the corresponding rate of convergence are evaluated using the double mesh principle [2]

Example 1 ([1]). Consider the constant coefficient problem with left boundary layer

subject to the interval and boundary conditions

The exact solution is given by

where

Example 2 ([18]). Consider the constant coefficient problem with right boundary layer

subject to the interval and boundary conditions

The exact solution is given by

where

Example 3 ([18]). Consider the variable coefficient problem with right boundary layer

subject to the interval and boundary conditions

The exact solution is not known.

Example 4 ([5]). Consider the variable coefficient problem with left boundary layer

with the interval conditions

The exact solution is not known.

 

7. Conclusion and discussion

Boundary value problem for second order singularly perturbed differential difference equations of mixed type having a boundary layer at one end (left or right) point is considered. To obtain an approximate solution for such type of problems, a mixed asymptotic-numerical method is proposed. In this method, the outer solution which corresponds to the reduced problem is solved numerically using OCI method while the inner solution is obtained analytically by making use of suitable stretching variable. Both cases, when the boundary layer occurs in the left and in the right side of interval are studied. The proposed method is analyzed for convergence, and the error estimate is also discussed.

The maximum absolute errors EN for different values of ε and N for Examples 1 and 2 are represented in Tables 1 and 2 respectively. The maximum absolute errors and the corresponding rates of convergence for different values of ε for Examples 3 and 4 are shown in Tables 3 and 4, respectively. The numerical results in these tables show that the proposed method is ε-uniformly convergent. Also, from Tables 1 and 2 the maximum absolute error between the numerical solution and exact solution stabilizes as ε → 0 for each N. While the maximum absolute error stabilizes and decreases rapidly with increasing N for each value of ε as indicated in Tables 3 and 4. Also, it can be observed that the rate of convergence is independent of the value of ε and δ. Comparisons in maximum absolute errors for problems in Examples 1 and 3 with the existing methods are presented in Tables 5 and 6. It is clear that the present method is robust with respect to the perturbation parameter and that our results are much better than the previous ones. The method presented here is very simple to implement, and with a little modification can be extended to high order SPDDEs and other types of differential difference equations like problems with discontinuous source terms and evolution equations.

Table 1.The maximum absolute errors EN for Example 1.

Table 2.The maximum absolute errors EN for Example 2.

Table 3.The values of for Example 3.

Table 4.The values of for Example 4 when δ = 0.5.

Table 5.Comparison of maximum absolute errors for Example 1 with [1].

Table 6.Comparison of maximum absolute errors for Example 3 for ε = 2−6 with [18].