SPLINE DIFFERENCE SCHEME FOR TWO-PARAMETER SINGULARLY PERTURBED PARTIAL DIFFERENTIAL EQUATIONS

Abstract

In this paper, we construct a numerical method to solve singularly perturbed one-dimensional parabolic convection-diffusion problems. We use Euler method with uniform step size for temporal discretization and exponential-spline scheme on spatial uniform mesh of Shishkin type for full discretization. We show that the resulting method is uniformly convergent with respect to diffusion parameter. An extensive amount of analysis has been carried out to prove the uniform convergence with respect to the singular perturbation parameter. The obtained numerical results show that the method is efficient, stable and reliable for solving convection-diffusion problem accurately even involving diffusion parameter.

1. Introduction

In the present paper, we are concerned with the numerical approximation of the partial differential equation:

Here QT = Ω × I, Ω ≡ (a, b), ∂Ω ≡ {a, b} , I ≡ (0, T), and a and b are real positive constants. We consider equation (1) associated with the Dirichlet boundary conditions

and the initial condition

where u(x, t) is the dependent variable, x and t are the independent variables and f , p, q and u0 are given functions of their arguments. It is usually assumed that 0 < εc ≪ 1 and 0 < εd ≪ 1.

This problem encloses both the reaction–diffusion problem when εc = 0 and the convection–diffusion problem when εc = 1. These problems occur naturally in various fields of science and engineering, such as, heat transfer with large Peclet numbers, combustion, nuclear engineering, control theory, elasticity, fluid mechanics, aerodynamics, quantum mechanics, optimal control, chemicalreactor theory, convection-diffusion process and geophysics; see [9,10,18].

There has been a lot of effort in developing numerical methods for the solution of singular perturbation problems that are uniformly convergent. Clavero et al.[1] and Kadalbajoo et al. [4] gives a uniform convergent numerical method with respect to the diffusion parameter to solve the one-dimensional time–dependent convection–diffusion problem. They used the implicit Euler method for the time discretization and the simple upwind finite difference scheme on a Shishkin mesh for the spatial discretization. Ramos [7] presented an exponentially fitted method for singularly perturbed, one-dimensional (convection–diffusion) parabolic problems, and showed its uniform convergence in the perturbation parameter. Surla and Jerkovic [13] considered a singularly perturbed boundary value problem using a spline collocation method. Sakai and Usmani [11,12] gave a new concept of B-spline in terms of hyperbolic and trigonometric splines which are different from the earlier ones. It is proved that the hyperbolic and trigonometric B-splines are characterized by a convolution of some special exponential functions and a characteristic function on the interval [0,1]. Sharma and Kaushik [14] presented a numerical approximation of a singularly perturbed time delayed convection diffusion problem on a rectangular domain. Finally, Zahra and El Mhlawy [18] developed a numerical method to solve two-parameter singularly perturbed semi-Linear boundary value problems. This method is based on exponential spline with a piecewise uniform Shishkin mesh which is shown to be uniformly convergent independent of mesh parameters and perturbation parameter εc and εd.

Our contribution in this paper is to introduce a full discretization scheme based on Rothe’s approximation in time discretization and on the exponential spline approximation in the spatial discretization for solving singularly perturbed, one-dimensional time-dependent convection–diffusion problem associated with homogeneous Dirichlet boundary conditions. It starts with the discretization in time by the 2-point Euler backward finite difference method in the time variable. After that we deal with the exponential spline method for the solution of the time discretized problem to compute the unknown function and the obtained system of algebraic equations is solved by iterative methods. Moreover, we obtain error estimates for the approximation process. The proposed technique is programmed using Matlab 7, and used in solving the proposed problem.

The rest of the paper is organized as follows: In section 2, we give some notations, assumptions and definitions, In §3, a time discretization scheme for the continuous problem is proposed in addition to the derivation of some a priori estimates. In §4, we prove the convergence results and error estimates of the semidiscretization scheme described in §3. We describe in §5 a full discretization scheme based on Rothe- exponential spline methods and the error estimates of the full discretized solution are derived. Finally, numerical results and discussions are presented and comparisons are made with other solutions in §6.

 

2. Notations, assumptions and definitions

In the sequel, we will denote by (. , .) either the standard inner product in L2 = L2(Ω) or the pairing between and V∗ ≡ H−1 (see e.g. [11]). we use the symbols |.|, ∥.∥ and ∥.∥∗ as the norms in L2(Ω), V, V∗, respectively. By we mean the strong and weak convergence. Also, we introduce some notations concerning the time discretization of our problem.

for any given family The following elementary relations will be used in the following analysis:

and Young’s inequality

where ε is a small constant. We will assume, throughout this work, the following hypotheses on the given data.

(H1) The function f : Ω × I → R is Lipschitz continuous in the sense of

with f(0, 0) = 0.

(H2) The functions p(x) and q(x) are smooth functions and satisfy

(H3)

Under these assumptions, we can define the weak solution of problem (1)-(2).

Definition 2.1. The measurable function u ∈ C(I; L2(Ω))∩L2(I; V ) with ∂tu ∈ L2(I; V∗) and u(x, 0) = u0(x) is said to be a weak (variational) solution of (1)−(3) if and only if the integral identity

holds for all φ ∈ V and a.e. t ∈ I.

Remark. It is evident that for any functions α, β ∈ H1(Ω) and γ ∈ L2(Ω)

 

3. The temporal discretization scheme, A prior estimate

Our main goal is to approximate (1)−(3) from a numerical point of view and to prove its convergence. The scheme starts with the discretization in time by the 2-point Euler backward finite difference method in the time variable, so that the problem is converted into a linear system of differential equations that easily solved numerically at each subsequent time level. Let m be a positive integer and subdivide the time interval I by the points tj, where tj= jτ , τ = T/m, j = 0, 1, · · · , m. The suggested discretization scheme of problem (8) consists of the following problem (in the weak sense):

where ∂x denotes to the derivative with respect to x and fj= f(x, tj).

The existence of a weak solution wj ∈ V is proved in [2]. By means of wj, (j = 0, 1, · · · , m) determined by the proposed scheme 10 -11 in each time step tj, we introduce the following piecewise linear functions (Rothe functions)

and the corresponding step function

Using the notation of Rothe function and its corresponding step function, a piecewise constant interpolation of equation 11 over I for each φ ∈ V yields

where tj − 1 ≤ t ≤ tj and with

In order to show the stability of the discrete solution and prove the convergence results, we shall derive some a priori estimates.

Lemma 3.1. Under the assumptions (H1)-(H3), there exists a positive constant C such that for any s,

Proof. Let us choose φ = τ δwj in 11 and summing over j from 1 to s, we obtain

Taking into consideration (4), the second and the fourth terms are estimated by

With the aid of 11, Younge’s inequality (5) and (9), we estimate the third term by

The last term is estimated by the use of (6) as follows

Collecting (16)-(20), using (H3), and choosing ε sufficiently small, we get

Applying Gronwall’s inequality, we conclude the proof.

Lemma 3.2. Uniformly with respect to n one has

Proof. The estimates (22) and (23) are a consequence of (15). By the use of the identity

the estimates (22)3, and (23) are a consequence of (15)2, (15)3 and the definitions of wm and and thus the proof completes.

 

4. Convergence results

This section is devoted to proving the convergence of the proposed scheme and estimating its accuracy. Before we are able to prove convergence we need to prove the compactness of in L2(I; L2(Ω)) which is a consequence of the following assertion.

Lemma 4.1. The estimate

holds uniformly for 0 < z < z0 and n.

Proof. [3] We sum up (12) for i = s+1, · · · , s+k considering φ = ws+k−ws)τ. Then we sum it up for s = 1, · · · , m − k and obtain the estimate

Hence for kτ ≤ z ≤ (k + 1)τ we conclude the desired estimate.

Theorem 4.1. There exists u ∈ L2(I; L2(Ω)) ∩ H1(I; L2(Ω)) such that

(in the sense of subsequences). Moreover, we have

Proof. The estimate (22)2 implies

Hence, from lemma 4.1, is compact in L2(I; L2(Ω)) because of Kolmogorov’s compactness argument. So we can conclude that wn → u and in L2(I; L2(Ω)) and also pointwise in Q. Also by the fact that we obtain that is weakly convergent in L2(I; V); i.e.

and this implies that

It remains to prove that ∂twm = ∂tu. For each t ∈ I, by lemma 3.3, ∂twm is uniformly bounded in the reflexive Banach space L2(I; L2(Ω)) and hence has a subsequence which converges weakly to an element ξ ∈ V∗ (Eberlin-Smulian theorem [6]). Thus

Using Fubini theorem, we get

Taking the limit as n → ∞, we obtain

which implies

Therefore, we have ∂tu = ξ. Using the above discussions and the fact that and hence and passing with n → ∞ in (14) the proof is complete.

To obtain an error estimate, let us start by introducing the following additional notation We now take the difference between (8) and (14) we easily obtain the equality for every φ ∈ V

Take φ = eu and write the new equation as I + II + III + IV = V. To begin with, we split the parabolic term into two parts I = I1 + I2, where

The second term of equation (37) is bounded by

The third and fourth terms of equation (37) are bounded by

In view of (H1) and Younge’s inequality we can estimate the right hand by

Collecting all the previous bounds, choosing ξ, ε, and κ sufficiently small and applying Gronwall lemma we obtain (28).

 

5. A Full discretization scheme

In this section, we use the exponential spline for the solution of discretization of problem (1-3) which is written in the form (denoting y(x) ≡ wj(x)):

with boundary conditions

in the time step tj, j = 1, 2, · · · , n, wj−1 is known and we have to solve the two parameter singularly perturbed boundary value problem

associated with the boundary condition (43). Here we use G(x) and H(x) to denote (q(x) + 1/τ ) and (fj + wj−1/τ ) respectively. To this end we consider a uniform mesh Δ with nodal point xi on [a, b] such that Δ : a = x1 < x2 < . . . < xm−1 < xm = b where

Let y(x) be the exact solution of the problem (44)-(43) and Si be an approximation to yi = y(xi) obtained by the segment Qi(x) passing through the points (xi, Si) and (xi+1, Si+1). Each mixed spline segment Qi(x) has the form, (for more details see [15-18]):

where ai,bi,ci and di are constants and k is a free parameter.

Following [18], we get the following relation for i = 1, 2, · · · , m − 1,

where:

θ = kh, α = (sinh(θ)−θ)/θ2 sinh(θ) and β = (2θ cosh(θ)−2 sinh(θ))/θ2 sinh(θ), when k → 0 that θ → 0 then and the relation defined by (47) reduces to an ordinary cubic spline relation:

at the point xi the proposed singularly perturbed problem may discretized by:

where, Gi = G(xi), Hi = H(xi) and pi = p(xi).

Substituting (49) into (47), we get the following linear equation for every i = 1, 2, · · · , m − 1

Equation (50) gives n−1 linear algebraic equations in n−1 unknowns Si,where Ai = α pi−1−β pi−3α pi+1,Di = 3 α pi−1+β pi−α pi+1,Ei = 4 α (pi+1− pi−1).

Now we investigate the convergence analysis of the suggested algorithm. The exponential spline solution of (1)-(3) is based on the linear equation given by (50) Let Y = (yi), S = (si), C = (ci), T = (ti) and E = (ei) = Y − S be m− 1 dimensional column vectors. It is easily seen that the system of equation (50) gives n − 1 linear algebraic equations in the n − 1 unknowns Si, i = 1, 2, · · · , m − 1. Then we can rewrite equation (50) in the standard matrix equation as:

Now the matrix B may be written as

where

and

We have the local truncation error are

From Eq. (52), then the row sum of the matrix B satisfies: for |pi| ≤ P ≥ 0,

for small values of εd and εc then the matrix B is irreducible and monotone and it follows that B−1 exist thus the problem (51) has a unique solution.

Lemma 5.1 ([5]). For any 0 < d < 1, we have up to a certain order q that it depends on the smoothness of the data

Shishkin Mesh Strategy

We form the piecewise-uniform grid in such a way that more points are generated in the boundary layer regions than outside of it. We divide the interval into three subintervals:

where transition parameters are given by

where

the real solutions of characteristic equation

with n to be the number of subdivision points of the interval [0,1] and we place n/4, n/2 and n/4 mesh points, respectively, in [0, γ1], [γ1, 1−γ2] and [1−γ2, 1]. Denote the step sizes in each subinterval by and respectively.

Where the resulting piecewise-uniform Shishkin mesh may be represented by:

where n is the number of discretization points and set of mesh points with:

Thus, a uniform mesh is placed on each of these subintervals.

Lemma 5.2 ([5]). The solution y(x) of (42)-(43) has the representation

Theorem 5.1. Let S(x) be the approximate solution of the solution y(x) of (42), then

where C is a constant independent of εd and εc.

Proof. The estimate is obtained on each subinterval Ωi = [0, 1] separately. Let any function z on Ωi, see[5]

and so it is obvious that, on Ωi

Taking maximum norm on both sides, we get

and by appropriate Taylor expansions it is easy to see that

From (58) and Lemma 5.1, on Ωi,

Also, using Lemma 5.3 and equation (58) on Ωi,

Case 1. The argument now depends on whether In this case ψ1 ≤ C ln n and ψ2 ≤ C ln n. Then the result follows at once from (58)

Case 2. When ln n and ln n. Suppose that i satisfies and Then and respectively. Now from (58)

If i satisfied Then γ1 ≤ xi and xi ≤ 1 − γ2 or γ2 ≤ 1 − xi and so

Using the above values in (59), we get the required result.

We are now in a position to formulate the main contribution of this work by combining (28) and (59) in the following theorem.

Theorem 5.2. Let the assumptions (H1)-(H2) be satisfied. Then uα → u in u ∈ H1(I; L2(Ω))∩L2(I; V ) for α → 0 where u is the variational solution of problem (1)-(3) and uα is constructed from the solution of problem (42)-(43). Moreover, we have the error estimates

 

6. Numerical Results

We present a case study to illustrate the performance of the error estimation procedures of section 5. For purposes of comparison, contrast and performance, two cases with known solutions were chosen.

Case study 1: Consider the following singularly perturbed parabolic partial differential equation, see [1].

ut − εduxx + ux = f(x, t), T = 1,in the interval [0, 1].

associated with the Dirichlet boundary conditions u(0, t) = 0, u(1, t) = 0. and the initial condition u0(x) = e−1/ε +(1− e−1/ε)x−e−(1−x)/ε, where f(x, t) = e−t(−c1 + c2(1 − x) + e−(1−x)/ε) . The exact solution is given by u(x, t) = e−t(c1 + c2x − e−(1−x)/ε), where c1 = e−1/ε, c2 = 1 − e−1/ε.

Case study 2: Consider the following singularly perturbed parabolic partial differential equation.

ut − εduxx + εcux = f(x, t), T = 1,in the interval [0, 1].

associated with the Dirichlet boundary conditions u(0, t) = 0, u(1, t) = 0. and the initial condition

where

TABLE 1.Maximum pointwise errors for Case study1

TABLE 2.Maximum pointwise errors and order of convergence for Case study1

The exact solution is given by

where

All the computations were performed by using MATLAB 7.

TABLE 3.Comparison of maximum pointwise errors for Case study1

TABLE 4.Maximum errors for Case study2 n = 128 and τ = 1.6/(n + 1)2

To compute the experimental rates of convergence Ordn for every fixedε, we use the rate of convergence from:

The numerical results displayed in Table 1-3 clearly indicate that the proposed method based on an exponential-spline function with Shishkin mesh is ε-uniformly convergent. Numerical solution profiles for case study are given in figure 1, respectively, for n = 128, τ = 0.02and different values of εd = 1, 2−6 and 2−12. For comparison, our results are better than the numerical solution presented in Clavero et al. [1]. Numerical solution profiles for case study2 are given in table 4 and also figure 2 and 3 for different values of εd and εc.

It is clearly that all solutions are in good harmony with each approximate and exact solution and the convergence obtained is parameter uniform and satisfies the theoretical predictions.

FIGURE 1.Exact and approximate solutions profiles of the case study1 for n = 128, τ = 0.02 and εd = 2−6and εd = 2−12 respectively

FIGURE 2.Exact and approximate solutions profiles of the case study2 for n = 128, 256, τ = 1.6/(n + 1)2, εc = 2−6, and εd = 2−6.

FIGURE 3.Exact and approximate solutions profiles of the case study2 for n = 64, 128, τ = 0.002, εc = 2−6, and εd = 2−6.

 

7. Conclusion

In this paper, we have investigated the application of Rothe’s approximation in time discretization and the exponential spline approximation in the spatial discretization to solve the convection–diffusion problems. The method is shown to be uniformly convergent i.e., independent of mesh parameters and perturbation parameters εd and εc. It has also been seen that the accuracy in the numerical results for the proposed scheme is comparable to that obtained by the uniformly convergent method with a uniform mesh [1]. It has been found that the proposed algorithm gives highly accurate numerical results and higher order of convergence than an upwind finite difference scheme. The present method provides a very simple and accurate method for solving convection–diffusion problems.