Numerical treatment of boundary value problems for second order singularly perturbed delay differential equations

Numerical treatment of boundary value problems for second order singularly perturbed delay differential equations

M.K. Kadalbajoo^I; K.K. Sharma^II

^IDepartment of Mathematics, IIT Kanpur, India, Professor, E-mail: [e-mail]

^IIDepartment of Mathematics, IIT Kanpur, India, Research Scholar, E-mail: [e-mail]

ABSTRACT

In this paper, the boundary value problems for second order singularly perturbed delay differential equations are treated. A generic numerical approach based on finite difference is presented to solve such boundary value problems. The stability and convergence analysis of the method is studied. The solution of the boundary value problems when delay is zero, exhibits layer behavior. Here, the study focuses on the effect of delay on the boundary layer behavior of the solution via numerical approach. The effect of the delay on the boundary layer behavior of the solution is shown by carrying out some numerical experiments.

Mathematical subject classification: G.1.7.

Key words: delay differential equation, singularly perturbed boundary value problem, boundary layer, oscillations.

1 Introduction

Delay differential equations play an important role in the mathematical modeling of various practical phenomena in the biosciences and control theory. Any system involving a feedback control will almost always involve time delays. These arise because a finite time is required to sense information and then react to it. A delay differential equation is of the retarded type if the delay argument does not occur in the highest order derivative term. If we restrict this class to a class in which the highest derivative term is multiplied by a small parameter, then it is said to be singularly perturbed delay differential equation of the retarded type. The boundary value problems for such a class of delay differential equations are ubiquitous in the modeling of several physical and biological phenomena like first exit time problem in modeling of activation of neuronal variability [8], in the study of bistable devices [1] and evolutionary biology [14], in a variety of models for physiological processes or diseases [11, 12, 14], to describe the human pupil-light reflex [10], variational problems in control theory [4, 5], and in describing the motion of the sunflower [13].

The literature on delay differential equations is mainly centered on first order initial value problems [3, 7]. But due to the versatility of boundary value problems for the second order delay differential equation in the mathematical modeling of processes in various kinds of application fields where they provide the best and sometimes the only realistic simulation of the observed phenomena, one cannot neglect this area, e.g., in Stein's model, the distribution representing inputs is taken as a Poisson process with exponential decay. If in addition, there are inputs that can be modeled as a Wiener process with variance parameter s and drift parameter µ, then the problem for expected first-exit time y, given initial membrane potential x Î (x₁, x₂), can be formulated as a general boundary value problem for a linear second order differential difference equation(DDE) [9]

where the values x = x₁ and x = x₂ correspond to the inhibitory reversal potential and to the threshold value of membrane potential for action potential generation, respectively. s and µ are variance and drift parameters, respectively, y is the expected first-exit time and the first order derivative term -xy¢(x) corresponds to exponential decay between synaptic inputs. The undifferentiated terms correspond to excitatory and inhibitory synaptic inputs, modeled as Poisson process with mean rates l_E and l_I, respectively, and produce jumps in the membrane potential of amounts a_E and a_I, respectively, which are small quantities and could be dependent on voltage. The boundary condition is

In this paper, we extend the numerical study of boundary value problems for second order singularly perturbed delay differential equations which was initiated in the paper [6] wherein the authors consider the case when delay d is o(e) and use Taylor's series to tackle the delay argument. But in the case when the delay d is of O(e), the method presented in the paper [6] fails due to the Taylor's series approximation for the term containing delay, which may lead to a bad approximation. To resolve that problem, we present here a numerical method composed of a standard upwind finite difference scheme on a special type of mesh to solve the boundary value problem for a singularly perturbed delay differential equation. To tackle the term containing delay, we construct a special type of mesh. We choose the mesh parameter h so that h = d/m, where m = p.mantissa of d, p > 1 is an integer. This method works nicely for either d = O(e) or d = o(e).

2 Statement of the problem

Consider a boundary value problem for the linear second order singularly perturbed differential equation with retarded argument

on 0 < x < 1, 0 < e << 1 and 0 < d < 1, subject to the interval and boundary conditions

where a(x), b(x), f(x) and f(x) are smooth functions and g is a constant. For d = 0, the solution of the boundary problem (2.1), (2.2) exhibits layer or oscillatory behavior depending on the sign of (a(x) + b(x)). If (a(x) + b(x)) < -M < 0, where M is a positive constant, the solution of the problem (2.1), (2.2) exhibits layer behavior and if (a(x) + b(x)) > M > 0, it exhibits oscillatory behavior. The boundary value problem considered here is of the reaction-diffusion type, therefore if the solution exhibits layer behavior, there will be two boundary layers which will be at both the end points {0,1}. In this paper, we study both the cases, i.e., when the solution of the problem exhibits layer as well as oscillatory behavior and show the effect of delay on the layer and oscillatory behavior. In particular, as delay increases then the layer behavior of the solution is destroyed and the solution begins to exhibit oscillatory behavior across the interval.

3 Numerical scheme

In this section, we discretize the boundary value problem (2.1), (2.2) using the method composed of a standard central finite difference operator on a special type of uniform mesh. To tackle the delay term, we choose the mesh parameter as h = d/m, where m = p.mantissa of d, p > 1 is an integer.

The difference scheme for the boundary value problem (2.1), (2.2) is given by

where D₊D_-y_i = (y_i-1 - 2y_i + y_i+1)/h². A simplification of Eq. (3.1) leads to a system of N - 1 difference equations in N - 1 unknowns

where

b_i = b(x_i), a_i = a(x_i) and f(x_i) = f_i for all i = 0, 1, ¼ N, and f(x_i) = f_i for all i = -m, -m + 1, ¼ 0.

4 Error estimate

4.1 Case I: Layer behavior

i.e., (a(x) + b(x)) < -M < 0, x Î .

Lemma 1 [Discrete Minimum Principle]. Suppose Y₀> 0 and Y_N> 0. Then L^N Y_i < 0 for all i = 1, 2, ¼ N - 1 implies that Y_i > 0 for all i = 0, 1, ¼ N.

Proof. Let k be such that Y_k = min_0<i<N Y_i and assume that Y_k < 0. Then we have Y_k - Y_k-1< 0, Y_k+1 - Y_k > 0 and for 1 < k < m

which contradicts the hypothesis that L^NY_i< 0 for all i = 1, 2, ¼ N - 1. Therefore our assumption that Y_k < 0 is wrong; hence Y_k> 0. We have chosen k fixed but arbitrary, so Y_i > 0 for all i, 0 < i < N.

Theorem 1. Under the assumptions that a(x)> 0 and (a(x)+b(x)) < -M < 0, where M is a positive constant, the solution of the difference scheme given by Eq. (3.1) with boundary conditions (3.2) exists, is unique, and satisfies

where C is a positive constant.

Proof. To prove the uniqueness and existence, suppose áu_iñ and áv_iñ be two solutions to the system of difference equations (3.1) with boundary conditions (3.2). Then suppose z_i = u_i - v_i is a mesh function. Clearly z₀ = 0 = z_N and for 1 < i < N - 1, we have

L^N z_i = L^N u_i - L^N v_i = 0, since u_i and v_i satisfy Eq. (3.1).

Thus the mesh function z_i satisfies the hypothesis of Lemma 1. Therefore an application of Lemma 1 to the mesh function z_i yields

Again if we set z_i = -(u_i - v_i), then once again z_i is a mesh function satisfying z₀ = 0 = z_N and L^N z_i = 0 for all i = 1, 2, ¼ N - 1. Thus, again an application of the discrete minimum principle for the mesh function z_i yields for i = 0, 1, ¼ N,

From Eqs. (4.1) and (4.2), we get z_i = 0 which implies the uniqueness of the solution to the difference scheme (3.1), (3.2). For linear equations, the existence is implied by uniqueness.

Now we shall prove the estimate. For that we introduce two barrier functions y^± defined by

where C > 1 is a constant.

Case i) For 1 < i < m, we have

Since a(x) and f(x) are bounded, so we choose the positive constant C so that the sum of the moduli of the first and second terms dominates the modulus of the third term in the above inequality. We then get

Case ii) For m < i < N - 1, we have

Combining both the cases, we obtain

An application of Lemma 1 to the mesh functions yields

which proves the required estimate.

4.2 Case II: oscillatory behavior

i.e., (a(x) + b(x)) >M > 0, x Î

Lemma 2 [Discrete Maximum Principle]. Suppose Y₀> 0 and Y_N > 0. Then L^NY_i > 0 for all i = 1, 2, ¼ N - 1 implies that Y_i> 0 for all i = 0, 1, ¼ N.

Proof. Let k be such that Y_k = max_0<i<N Y_i and assume Y_k < 0. Then we have Y_k - Y_{k -1}> 0, Y_{k +1} - Y_k< 0 and for 1 < k < m

which contradicts the hypothesis of Lemma 2. Therefore our assumption that Y_k < 0 is wrong; hence Y_k> 0. We have chosen k fixed but arbitrary, so Y_i> 0 for all i, 0 < i < N.

Theorem 2. Under the assumptions that a(x) > 0, b(x) > 0 and (a(x) + b(x)) > M > 0, where M is a positive constant, the solution of the difference scheme given by Eq. (3.1), (3.2) exists, is unique, and satisfies

where C > 1 is a constant.

Proof. To prove the uniqueness and existence, suppose áu_iñ and áv_iñ be two solution to the discrete problem (3.1), (3.2). Suppose z_i = u_i - v_i is a mesh function. Then

z₀ = 0 = z_N

and for 1 < i < N - 1, we have

L^N z_i = L^N u_i - L^N v_i = 0, since u_i and v_i satisfy Eq. (3.1).

Thus the mesh function z_i satisfies the hypothesis of Lemma 2; therefore an application of Lemma 2 to the mesh function z_i gives

If we set z_i = -(u_i - v_i), clearly z_i is a mesh function satisfying z₀ = 0 = z_N and L^N z_i = 0 for all i = 1, 2, ¼ N - 1. So again an application of the discrete maximum principle for the mesh function z_i yields

From Eqs. (4.6) and (4.7), we get u_i - v_i = 0, i = 0, 1, ¼ N which proves the uniqueness of the solution to the discrete problem (3.3), (3.2). For linear equations, the existence is implied by uniqueness.

Now the required bound on the solution remains to be proved and to obtain this, we use the barrier functions y^± defined by

where C is a constant and C > 1.

Case i) For 1 < i < m, we have

Using the inequality b_i M^-1> 1 in the above inequality (4.8) followed by a simplification yields

In the above inequality (4.9) the first and second terms are positive and a(x) and f(x) are bounded so we choose the constant C so that the total of the first two terms dominates the modulus of the third term which gives

Case ii) For m < i < N - 1, we have

Combining both the cases, we obtain

An application of Lemma 2 for yields

which proves the required estimate.

Theorems 1 and 2 imply that the solution to the discrete problem (3.1), (3.2) in both the cases (i.e., either when the solution of the problem exhibits layer or oscillatory behavior) is uniformly bounded, independently of the mesh parameter h and the parameter e, which proves that the difference scheme is stable for all mesh sizes.

Corollary 1. The error e_i= y(x_i) - y_i between the solution y(x) of the continuous problem (2.1), (2.2) and the solution áy_iñ of the discretized problem (3.1), (3.2) satisfies the estimates

where T_i satisfies

T_i< (h²e²/12)|y^iv(x)|.

Proof. The truncation error T_i is given by

T_i = e²[(y_i-1 - 2y_i + y_i+1 - y¢¢(x_i)].

Now using Taylor's series and after some simplifications, we obtain

T_i< |y^iv(x)|.

We have

L^Ne(x_i) = L^N y(x_i) - L^N(y_i) = T_i, i = 1, 2, ¼ N - 1

and e₀ = 0 = e_N. Then by using Theorem 1 and Theorem 2, we obtain the required error estimate on the error in both the cases (i.e., when the solution exhibits layer as well as oscillatory behavior).

5 Computational results

Example 1 (pp. 263, [9]) e²y¢¢(x) - 2y(x - d) - y(x) = 1, under the interval and boundary conditions

y(x) = 1, -d <x< 0, y(1) = 0.

Example 2 (pp. 265, [9]) e²y¢¢(x) + 0.25y(x - d) - y(x) = 1, under the interval and boundary conditions

y(x) = 1, -d <x< 0, y(1) = 0.

Example 3 (pp. 265, [9]) e²y¢¢(x) + 0.25y(x - d) + y(x) = 1, under the interval and boundary conditions

y(x) = 1, -d <x< 0, y(1) = 0.

Example 4 (pp. 265, [9]) e²y¢¢(x) + y(x - d) + 2y(x) = 1, under the interval and boundary conditions

y(x) = 1, -d <x< 0, y(1) = 0.

For d = 0, the solution to the boundary value problem (2.1), (2.2) exhibits layer or oscillatory behavior according to the sign of the coefficient of the reaction term. To demonstrate the effect of delay on the layer and oscillatory behavior of the solution and efficiency of the method, we consider several numerical examples and solve them using the proposed method. In examples 1 and 4, the coefficient of the delay term is of O(1) while that of o(1) in examples 2 and 3. We have plotted the graphs of the solutions of the problems for e = 0.01 with different values of d to show the effect of delay on the boundary layer or oscillatory behavior of the solution. The maximum absolute error for the considered examples is calculated using the double mesh principle [2] (as the exact solutions for the considered examples are not available), as shown in Tables 1 and 2.

6 Discussion

A numerical study of boundary value problems for second order singularly perturbed differential difference equations with delay as well as advance is initiated in paper [6]. In [6], the authors gave a numerical scheme to solve such boundary value problems in the case when the delay is of small order of e, i.e., d = o(e) but the method fails in the case when the delay is of capital order of e, i.e., d = O(e). In this paper, we present a generic numerical approach to solve the boundary value problems for second order delay differential equations which works nicely in both the cases, i.e., when d = o(e) and d = O(e).

To show the effect of delay on the boundary layer or oscillatory behavior of the solution, several numerical experiments are carried out in section 5. We observe that when the coefficient of the delay term in the equation is of O(1) and the delay term is of O(e), the layer behavior of the solution is no longer preserved and the solution exhibits oscillatory behavior. Not only is the layer behavior destroyed, but also oscillations previously confined to the layer region are extended throughout the entire interval [0, 1]. As the delay increases, the amplitude of the oscillations increases as shown in Figures 2-5 for example 1. If the order of the coefficient of the delay term in the equation is of o(1), the delay affects the boundary layer solution but maintains the layer behavior, although the delay is of O(e) as shown in Figures 6-9 for example 2. From Figure 1, we observe that when the delay is of o(e), the solution maintains layer behavior although the coefficients in the equation are of O(1) and as the delay increases, the thickness of the left boundary layer decreases while that of the right boundary layer increases.

To demonstrate the effect on the oscillatory behavior, we consider the examples 3 and 4 when the solution of the problem exhibits oscillatory behavior for delay equal to zero. We observe that if the coefficient of the delay term is of o(1), the amplitude of the oscillations increases slowly as the delay increases provided the delay is of o(e) (Figure 10) and increases exponentially as the delay increases if the delay is of O(e) (Figures 11, 12).

Received: 10/IV/03. Accepted: 31/VIII/04

#578/03.

• [1] M.W. Derstine, F.A.H.H.M. Gibbs and D.L. Kaplan, Bifurcation gap in a hybrid optical system, Phys. Rev. A 26 (1982), pp. 3720-3722.
• [2] E.P. Doolan, J.J.H. Miller and W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, (1980).
• [3] R.D. Driver, Ordinary and Delay Differential Equations, Belin-Heidelberg, New York, Springer, (1977).
• [4] V.Y. Glizer, Asymptotic solution of a singularly perturbed set of functional-differential equations of riccati type encountered in the optimal control theory, Differ. Equ. Appl., 5 (1988), pp. 491-515.
• [5] V.Y. Glizer, Asymptotic solution of a boundary-value problem for linear singularly-perturbed functional differential equations arising in optimal control theory, J. Optim. Theory Appl., 106 (2000), pp. 309-335.
• [6] M.K. Kadalbajoo and K.K. Sharma, Numerical analysis of boundary-value problems for singularly-perturbed differential-difference equations with small shifts of mixed type, J. optim. Theory Appl., 115 (2002), pp. 145-163.
• [7] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, INC, (1993).
• [8] C.G. Lange and R.M. Miura, Singular perturbation analysis of boundary-value problems for differential-difference equations. vi. small shifts with rapid oscillations, SIAM J. Appl. Math., 54 (1994), pp. 273-283.
• [9] C.G. Lange and R.M. Miura, Singular perturbation analysis of boundary-value problems for differential-difference equations. v. small shifts with layer behavior, SIAM J. Appl. Math., 54 (1994), pp. 249-272.
• [10] A. Longtin and J. Milton, Complex oscillations in the human pupil light reflex with mixed and delayed feedback, Math. Biosciences, 90 (1988), pp. 183-199.
• [11] B.J. MacCartin, Exponential fitting of the delayed recruitment/renewal equation. J. Comput. Appl. Math., 136 (2001), pp. 343-356.
• [12] M.C. Mackey and L. Glass, Oscillations and chos in physiological control systems, Science, 197 (1977), pp. 287-289.
• [13] M.L. Pena, Asymptotic expansion for the initial value problem of the sunflower equation, J. Math. Anal. Appl., 143 (1989), pp. 471-479.
• [14] M. Wazewska-Bzyzewska and A. Lasota, Mathematical models of the red cell system, Mat. Stos., 6 (1976), pp. 25-40.

Publication Dates

• Publication in this collection
  07 Nov 2005
• Date of issue
  Aug 2005

History

• Accepted
  31 Aug 2004
• Received
  10 Apr 2003