Abstract

In this paper, an initial value method for solving a class of linear second-order singularly perturbed differential difference equation containing mixed shifts is proposed. In doing so, first, the given problem is modified in to an equivalent singularly perturbed problem by approximating the term containing the delay and advance parameters using Taylor series expansion. From the modified problem, two explicit initial value problems which are independent of the perturbation parameter are produced; namely, the reduced problem and the boundary layer correction problem. These problems are then solved analytically and/or numerically, and those solutions are combined to give an approximate solution to the original problem. An error estimate for this method is derived using maximum norm. Several test problems are considered to illustrate the theoretical results. It is observed that the present method approximates the exact solution very well.

1. Introduction

Singularly perturbed differential equations (SPDEs) refer to the study of classes of differential equations containing small parameter(s) multiplying the highest derivative(s). A well-known fact is that the solutions of such problems have a multiscale character, i.e., there are thin transition layer(s) where the solution varies very rapidly, while away from the layer(s), the solution behaves regularly and varies slowly. This leads to boundary and/or interior layer(s) in the solution of the problems [1]. Due to the presence of the layer regions, it has been shown that the classical numerical methods fails to produce good approximations for SPDEs. In fact, some numerical techniques employed for solving singularly perturbed boundary value problems (SPBVPs) are based on the idea of replacing these problems by suitable initial value problems (IVPs).The reason for this is that the numerical treatment of a boundary value problem is much more demanding than the treatment of the corresponding IVPs. There are different initial value methods in the literature of SPDEs developed for solving SPBVPs; for the detailed discussions of such methods, one can refer [1, 2] and the references therein.

Modeling automatic engines or physiological systems often involve the idea of control because feedback is used in order to maintain a stable state. However, much of this feedback require a finite time to sense information and react to it. A popular way to describe this process is to formulate delay differential equations or differential difference equations (DDEs) where the evolution of a dependent variable is a function of time which depends on not only current time but also earlier time [3]. A singularly perturbed differential difference equation is a differential equation in which the highest derivative is multiplied by a small parameter and which involves at least one shift term. Such type of equations arise frequently in the mathematical modeling of various practical phenomena, such as in optical bistable devices [4], in variational problems in control theory [5, 6], in the hydrodynamics of liquid helium [7], in description of the human pupil-light reflex [8], in microscale heat transfer [9], and in a variety of models for physiological processes or diseases [10, 11].

Singularly perturbed differential difference equations (SPDDEs) have not been satisfactorily discussed in the literature up to now; however, in the past few decades, there has been a growing interest in the numerical and/or asymptotic study of such problems. Lange and Miura initiate the study of boundary value problems for SPDDEs in a series of papers [1214] and discuss the case of small as well as large shift parameters. Tain [15] extended the concept of singular perturbation theory for ordinary differential equations to delay differential equations, applying O’Malley–Hoppensteadt technique to obtain approximate solutions. Kadalbajoo and Sharma [16] constructed an ε-uniform fitted mesh method for solving singularly perturbed differential-difference equations with mixed type of shifts. Kadalbajoo and Sharma [17] described a numerical approach based on finite differences to solve a mathematical model arising from neuronal variability. Patidar and Sharma [18] approximated the term containing delay by Taylor series expansion and then applied an ε-uniformly convergent nonstandard finite difference methods to SPDDEs with small delay. Kumar and Sharma [19] presented a numerical technique to approximate the solution of SPDDEs with delays as well as advances. In recent years, different scholars further developed numerical schemes for ordinary differential-difference equations with mixed shifts, to mention few [2024].

In this paper, we consider linear second-order singularly perturbed differential difference equation containing mixed shifts (i.e., both delay and advance parameters) on the reaction terms and propose an asymptotic numerical method, namely, initial value technique. This initial value method was first introduced by EL-Zahar and EL-Kabier [2]. In fact, they applied this method for solving singularly perturbed differential equations without shift parameters. In [25] this method was extended for singularly perturbed delay differential equations having a negative shift on the convection term. The proposed method here is similar in some respect to the methods in [2, 25] but differs in detail, and the model problem considered is completely different.

The rest part of this paper is organized as follows. In Section 2, the problem under consideration is stated. Section 3 is devoted on the analytic properties of the exact solution. The proposed method is described in Section 4, and the error analysis for the method is derived in Section 5. Several numerical examples are given in Section 6. Finally, discussions on the numerical results and conclusions are included in Section 7.

2. Statement of the Problem

Consider the following second-order linear SPDDE with a delay and advance term:

, under the following interval conditions:where , , , , , , and are bounded and continuously differentiable functions on , is the singular perturbation parameter, and and are the delay and advance parameters, respectively.

When the shifts are zero (i.e., ), the solution of the resulting problem exhibits layer behavior or turning point behavior depending on the coefficient of the convection term, i.e., whether does not change sign or changes sign on the interval . The layer behavior of the problem under consideration is maintained for and but only if they are sufficiently small. In general, the solution of problem (1)–(3) exhibits boundary layer behavior at one end of the interval depending on the sign of [16].

By using Taylor series expansion on the terms and in (1), we have

Using (4) and (5) in (1), we obtainwhere

Since (6) is an approximate version of (1)–(3), it is good to use different notation (say ) for the solution of this approximate differential equation. Thus, (6) can be rewritten aswhich differs from the original problem (1)–(3) by terms.

The approximation of (1)–(3) by (9) is acceptable, because of the condition that and are sufficiently small. This replacement is significant from the computational point of view. Further details on the validity of this transition can be found in Els’golts and Norkin [26]. Thus, the solution of (9) will provide a good approximation to the solution of (1)–(3). Further, we assume thatthroughout the interval , where M and θ are some positive constants. Under these assumptions, (9) has a unique solution which exhibits a boundary layer of width on the left side of the underlying interval [24].

Remark 1. In this paper, we consider only the case where there is one boundary layer at the left end of the interval. The case when the layer occurs at right end can be analyzed similarly. However, we present some test problems for the latter case.

3. Some a Priori Estimates

Throughout this paper, C (sometimes subscripted) denote generic positive constants independent of the singular perturbation parameter ε and in the case of discrete problems, also independent of the mesh parameter N. The maximum norm is used for studying the convergence of the approximate solution to the exact solution of the problem:

Now, we give the required bounds on the solution u of (9) that will be used to establish the error estimate. First, we consider the following property of the operator L.

Lemma 1. (continuous minimum principle).
Assume that is any sufficiently smooth function satisfying and . Then, for all implies that for all .

Proof. Let be an arbitrary point in such that and assume that . Clearly ; therefore, and . Moreover,which is a contradiction. It follows that our assumption is wrong. So, . Since is an arbitrary point, , .

Lemma 2. (stability result).
Let be the solution of the problem (9). Then

Proof. First we construct two barrier functions defined byThenTherefore, from Lemma 1, we obtain for all , which gives the required estimate.
Lemma 1 implies that the solution is unique and since the problem under consideration is linear, the existence of the solution is implied by its uniqueness. Furthermore, Lemma 2 gives the boundedness of the solution.

Lemma 3. The derivatives of the solution of the boundary value problem (9) satisfy the following estimates for :

Proof. For the proof of this lemma, the reader can refer to [17].
Lemma 4. The solution of (9) admits the following decomposition:where the regular (smooth) component satisfiesand the singular component satisfies

Proof. For the proof of this theorem, the reader can refer to [16].

4. Description of the Method

Here, we extend the initial value method proposed and used in the articles [2, 25] to approximate the solution of (9), when the solution of this problem exhibits a boundary layer near the left end (i.e., near at ). However, the procedure we follow here is somehow different.

First, we rewrite (9) equivalently aswhere

Let be the solution of the reduced problem (i.e., the problem which is obtained from (9) when ). So, (9) is reduced to an initial value problem of the following form:

If this problem is separable, then it can be integrated easily to give an exact solution, and if not, any initial value solver like fourth-order Runge–Kutta method can be used to approximate the solution. For the exact solution of the reduced problem, the following theorem gives an error bound.

Theorem 1. Let be the solution of (9) and be its reduced problem solution defined by (22). Then

Proof. First consider the following two barrier functions defined byThenfor an appropriate choice of C, i.e., .
Therefore, from Lemma 1, we obtain for all , which gives the required estimate.

Remark 2. From the above theorem, it is clear that the solution u of the boundary value problem (9) exhibits a strong boundary layer at and further away from the boundary layer region and in particular on , where , for sufficiently small values of ε, we haveFrom Theorem 1, it can be observed that the solution satisfies (9) on most part of the interval and away from . Thus, by replacing the solution by on the right part of (20), we obtain an asymptotically equivalent approximation aswhereIntegrating both sides of (27) with respect to x giveswhereUsing (22) in the above integral yieldsThen, substituting this in to (29) gives uswhere k is an integration constant. In order to determine k, we introduce the condition that the reduced equation of (31) should satisfy the boundary condition at . Thus, we get .
Hence, by substituting in (31), a first-order initial value problem which is asymptotically equivalent to the second-order boundary value problem (9) is obtained and written as follows:Over most of the domain , the solution of the reduced problem (22) behaves like the solution of (9). But, in the neighborhood of , there is a region in which this solution varies greatly from the solution of (9). To compensate the solution over this region (inner layer), a new inner variable is introduced by stretching the spatial coordinate x asUsing this stretching transformation in to (33), we obtainIn spite of this simplification, equation (34) remains a first-order differential equation and also regularly perturbed. For , it becomesThis is a differential equation for the solution of the layer region. Moreover, the solution of (36) is supposed to counteract for the fact that the solution of the reduced problem does not satisfy the boundary condition at and that this solution satisfiesFurther, using the substitution in to (36), we obtain the following boundary layer correction problem:This equation is a linear initial value problem with constant coefficient, which can easily be solved analytically, and givesFinally, from standard singular perturbation theory, it follows that the solution of the IVP (33) admits the representation in terms of the solutions of the reduced and boundary layer correction problems, which approximates the modified problem (9), that is,

Theorem 2. Let be the approximate solution of (9) given by (40), then it satisfies the following bound:where is the exact solution of (9).

Proof. First consider the following two barrier functions defined byThenfor an appropriate choice of C, i.e., , andfor an appropriate choice of C, i.e., .
Therefore, from Lemma 1, we obtain for all , which gives the required estimate.

5. Error Analysis

The numerical error of the present method has three sources: one from Taylor’s series approximation of the original problem (1)–(3), the second from the asymptotic approximation of the modified problem (9), and the last from the numerical approximation of the reduced problem (22). Let h be the mesh spacing of the domain of the problem.

5.1. Error on the Nonboundary Layer Domain

Let be the solution of the modified problem (9), be the exact solution of the reduced problem (22), and be the numerical solution of the reduced problem obtained from the fourth-order Runge–Kutta method. On the nonboundary layer domain, the error is

By using (26), (45) and triangle inequality, we conclude

Most of the time, the exact solution of the reduced problem can be easily obtained and the second term of the above error inequality is vanished.

5.2. Error on the Boundary Layer Domain

On the boundary layer domain, the asymptotic approximation error is generated from the reduction of order method and the numerical error from the numerical approximation of the outer solution , since the initial condition of the boundary layer correction problem (38) is affected by .

Let be the exact solution of the (37) and be the numerical solution of (37). On the boundary layer domain, the error is

The new method works well for singular perturbation problems since the singular perturbation parameter ε is extremely small.

6. Test Problems and Numerical Results

To demonstrate the applicability of the proposed method, we have implemented it on four boundary value problems of the form (1)–(3), exhibiting either right or left boundary layer. Three of them are constant coefficient problems, whereas the remaining one is of variable coefficients.

In case exact solution of the test problem is available for comparison, the pointwise error and maximum absolute error are calculated by using the following formula:where is the exact solution of the given problem and is the approximate solution obtained by using the proposed method, evaluated at N equally spaced mesh points.

In case the exact solution is not available, the pointwise error and maximum absolute error are calculated by using the double mesh principle given bywhere and denote the and components of the numerical solutions obtained by using N and meshes points, respectively. In addition, the corresponding rate of convergence is determined by

The exact solution of the BVPs (1)–(3) having constant coefficients (i.e., , , , , , and ), is given bywhere

Since both the reduced problem (22) and the layer correction problem (38) for constant coefficient BVPs are separable IVPs which can easily be solved analytically, and the proposed method gives the following formula for the approximate solution:

Following a similar derivation like that of Section 4 for constant coefficient problems exhibiting right boundary layer, we obtain

Example 1. Consider the following BVP with left boundary layer:withThe exact solution of this problem is given by (51)–(54) and the approximate solution using (55) becomesThe maximum pointwise errors of Example 1 for different values of ε, δ, and are given in Table 1. Plots of the approximate solution for different values of δ with fixed and are displayed in Figure 1.

Example 2. Consider the following nonhomogeneous BVP with left boundary layer:withThe exact solution of this problem is given by (51)–(54) and the approximate solution using (55) becomesThe maximum pointwise errors of Example 2 for different values of ε, δ, and η with are given in Table 2. Plots of the approximate solution for different values of δ with fixed , , and are displayed in Figure 2. And, Figure 3 shows the plots of the approximate solutions for , , and with different values of η.

Example 3. Consider the following BVP with right end boundary layer:withThe exact solution of this problem is given by (51)–(54) and the approximate solution using (56) becomesThe maximum pointwise errors of Example 3 for different values of ε, δ, and η with are given in Table 3. Plots of the approximate solution for different values of δ with fixed , , and are displayed in Figure 4. In Figure 5, the plots of the approximate solutions for , , and with different values of η are displayed.

Example 4. Consider the following variable coefficient BVP:withExact solution: not known.
Property: boundary layer near right end.
Solution: using Taylor’s series expansion, this problem can be written equivalently aswhereSince , , and , then and which implies the given problem has a boundary layer near right end. Then, the reduced problem becomesSince this problem is difficult to solve analytically, we apply classical fourth-order Runge–Kutta method to compute the approximate solution .
Next, the boundary layer correction problem becomeswhere , then solving this IVP analytically gives Finally, the approximate solution of the given problem becomesThe maximum errors and rate of convergences for Example 4 for different values of ε and N by taking are given in Table 4. The plots of the approximate solution for and with and are displayed in Figure 6.

7. Discussion

In this article, an initial value method for solving linear second-order singularly perturbed differential difference equation with both delay and advance parameters is considered. First, by applying Taylor’s series expansion on the term containing the delay and advance, the given problem is modified to an asymptotically equivalent singularly perturbed problem. Then, the solution of the modified problem is computed analytically and/or numerically by solving two initial value problems, namely, the reduced problem and the boundary layer correction problem, which are independent of the perturbation parameter ε. The method is simple to apply, very easy to implement on a computer, and offers a relatively simple tool for obtaining approximate solution.

To show the efficiency and applicability of the proposed method, four test problems are considered. For problems with exact solution, the maximum error is calculated and tabulated in Tables 13 for different values of ε, δ, and η. On the other hand, for the problem with no exact solution, the double mesh principle is used to compute the maximum error and rate of convergence, and the results are displayed in Table 4. From the results, it is observed that, for very small ε, the present method approximates the exact solution of the problems very well, and also it is superior to some of the existing methods in the literature such as [17, 20, 21, 24].

In addition, to examine the effect of the small shifts on the behavior of the solution, graphs of the solution for the test problems are displayed in Figures 15. We observe from Figures 13 that, when the solution exhibits left layer, the effect of both shifts on the solution in the layer region is negligible whereas that in the outer region is considerable. On the other hand, Figures 4 and 5 illustrate that, when the solution exhibits right layer, the changes in delay or advance affect the solution in layer region as well as outer region. But, the thickness of the layer increases as the size of the delay increases while it decreases as the size of the advance increases.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This research work was partially supported by the Bahir Dar University, College of Science.