Parameter-uniform approximations for a singularly perturbed convection-diffusion problem with a discontinuous initial condition

A singularly perturbed parabolic problem of convection-diffusion type with a discontinuous initial condition is examined. A particular complimentary error function is identified which matches the discontinuity in the initial condition. The difference between this analytical function and the solution of the parabolic problem is approximated numerically. A co-ordinate transformation is used so that a layer-adapted mesh can be aligned to the interior layer present in the solution. Numerical analysis is presented for the associated numerical method, which establishes that the numerical method is a parameter-uniform numerical method. Numerical results are presented to illustrate the pointwise error bounds established in the paper.


Introduction
In this paper, we examine a singularly perturbed convection-diffusion problem with a discontinuous initial condition of the form: Findû such that −εû ss +âû s +bû +û t =f , (s, t) ∈Q := (0, 1) × (0, T ]; (1a) u(s, 0) = φ(s) ∈ C 0 (0, 1);â > 0;b ≥ 0, with Dirichlet boundary conditions. As this is a parabolic problem, an interior layer emerges from the initial discontinuity, which is diffused over time if ε = O(1). However, when the parameter is small, the interior layer is convected along a characteristic curve associated with the reduced problem. In [8], we examined a related singularly perturbed reaction-diffusion problem (setâ ≡ 0 in (1)) with a discontinuous initial condition and we used an idea from [3] to first identify an analytical function which matched the discontinuity in the initial condition and also satisfied a constant coefficient version of the differential equation. A numerical method was then constructed to approximate the difference between the solution of the singularly perturbed reaction-diffusion problem and this analytical function. The numerical approximation involves approximating an interior layer function whose location, in the case of a reaction-diffusion problem, is fixed in time.
In the corresponding convection-diffusion problem, the location of the interior layer function moves in time and, from [5], we know that the numerical method needs to track this location. Shishkin [10] examined problem (1) in the case where the initial condition φ ∈ C 0 (0, 1) \ C 1 (0, 1). In [11, Chapter 10 and §14.2], Shishkin and Shishkina discuss the method of additive splitting of singularities for singularly perturbed problems with non-smooth data. We follow the same philosophy here.
When the convective coefficient depends solely on time (â(s, t) ≡â(t) > 0), the main singularity generated by the discontinuous initial condition can be explicitly identified by a particular complimentary error function. This error function tracks the location of the interior layer emanating from the discontinuity in the initial condition and it also satisfies the homogenous partial differential equation (1a) exactly. When this discontinuous error function is subtracted from the solutionû of (1), the remaining function (denoted below byŷ) contains no interior layer and it can be adequately approximated numerically by designing a numerical method which incorporates a Shishkin mesh in the vicinity of the boundary layer [6].
In this paper we deal with the more general case of the convective coefficient depending on both space and time. In this case, the situation is more complicated. The main singularity is again a particular complimentary error function which tracks the location of the interior layer, but when the coefficientâ in (1a) varies in space this complimentary error function does not satisfy the homogenous partial differential equation (1a). Moreover, when this discontinuous error function is subtracted from the solutionû of (1), the remaining functionŷ(s, t) contains its own interior layer. To generate an accurate numerical approximation to this remainderŷ, a coordinate transformation is first required in order that a mesh can be constructed to track the location of this internal layer. Hence the numerical method used to approximate the remainder (whenâ depends on space and time) is different to the numerical method used to approximate the remainder in the case of the convective coefficient solely depending on time. Needless to say, the more general method can also be applied to the case where the convective coefficientâ is independent of space. If the coordinate transformation is not used, in the numerical section we demonstrate that one does not generate a parameter-uniform approximation ifâ depends on the space variable.
In §2 we specify the continuous problem and deduce bounds on the partial derivatives of the solution. Some of the more technical details involved in the proofs of the bounds on the continuous solution are presented in the appendices. A piecewise-uniform mesh is constructed in §3, which is designed to be refined in the neighbourhood of the curve Γ * , which identifies the location of the interior layer at each time. To analyse the parameter-uniform convergence of the resulting numerical approximations on such a mesh, it is more convenient to perform the analysis in a transformed domain where the location of the interior layer is fixed in time. To simplify the discussion of the method and the associated numerical analysis, we discuss the case where there is no source term present in the problem in §2 and §3. In §4, we outline the modifications required when a source term is present. In §5, we present some numerical results to illustrate the performance of the method.
Notation: Throughout the paper, C denotes a generic constant that is independent of the singular perturbation parameter ε and all the discretization parameters. The L ∞ norm on the domain D will be denoted by · D and the subscript is omitted if D =Q. We also define the jump of a function at a point d by Functions defined in the computational domain will be denoted by f (x, t) and functions defined in the untransformed domain will be denoted byf (s, t).

Continuous problem
Consider the following convection-diffusion problem 1 : Findû such that In general, a moving interior layer and a boundary layer will appear in the solution. When the convective term depends on space then the path of the characteristic curveΓ * (associated with the reduced problem) is implicitly defined byΓ Since we have assumed thatâ > 0, the function d(t) is monotonically increasing. We restrict the size of the final time T so that the interior layer does not interact with the boundary layer. Thus, we limit the final time T 2 such that In the error analysis, we are required to impose a further restriction on the final time by assuming that 2T δ â s ≤ 1 − γ, 0 < γ < 1.
(2j) 1 As in [4], we define the space C 0+γ (D), where D ⊂ R 2 is an open set, as the set of all functions that are Hölder continuous of degree γ ∈ (0, 1) with respect to the metric · , where for all pi = (xi, ti), ∈ R 2 , i = 1, 2; p1 − p2 2 = (x1 − x2) 2 + |t1 − t2|. For f to be in C 0+γ (D) the following semi-norm needs to be finite f 0+γ,D := sup The space C n+γ (D) is defined by and · n+γ , · n+γ are the associated norms and semi-norms. The discontinuity in the initial condition generates an interior layer emanating from the point (d.0). By identifying the leading term 0.5[φ](d)ψ 0 in an asymptotic expansion of the solution, we can define the continuous function whereLŷ Note that in (2g) we impose the constraint [φ ](d) = 0 on the initial condition. This assumption permits us to complete the analysis of the numerical error. Based on the expansion (33) of the solution derived in the appendix, we note thatψ i ∈ C i−1 (Q), i ≥ 1, which implies (due to assumption (2g)) thatŷ(s, t) ∈ In addition, in (2g) we also assume thatâ s (d, 0) = 0, which results in the interior layer function (defined in (9)) being sufficiently regular to establish the bounds (12). The constraint (2j) is used in establishing the pointwise bound (11) on the interior layer function. This bound is used to determine the transition points in the Shishkin mesh around the interior layer. Finally, for sufficiently smooth and compatible boundary conditions at (0, 0) and (1,0), there is no loss in generality in assuming the constraints (2c), as the simple subtraction of the linear functionq(s, t) :=û(0, t)(1 − s) +û(1, t)s fromû leads us to problem (2) withf replaced byf 1 :=f −Lq.
Observe that the inhomogeneous term in (4) is continuous, but not in C 1 (Q) on the closed domain. The presence of this inhomogeneous term will induce an interior layer into the functionŷ. So if the convective coefficient a(s, t) depends on the space variable, we are required to transform the problem (2) so that the curveΓ * is transformed to a straight line, around which a piecewise-uniform Shishkin mesh is constructed.
Using this map the problem to solve numerically, transforms into the problem: Find y such that where Ly := −εy xx + κ(x, t)y x + g(x, t)y t , and the coefficients are Observe that g is a discontinuous function along x = d and [g](d, t) < 0 for all t > 0. In addition, for all t ≥ 0, |g − 1| ≤ C|d(t) − d| ≤ Ct and The transmission condition We associate the following differential operator with this transformed problem. For this operator L ε a comparison principle holds [5].
Using this comparison principle we see from (8) that The solution of problem (6) can be decomposed into the sum of a regular component v, a boundary layer component w, a weakly singular component and an interior layer z component: In Appendix B, the regular componentv ∈ C 4+γ (Q) and the boundary layer componentŵ ∈ C 4+γ (Q) are defined in the original variables (s, t).
The mapping X : (s, t) → (x, t) defined in (5) is not smooth along the interface x = d. Hence, in the transformed variables the regular component v is defined so that v ∈ (C 4+γ (Q + ) ∪ C 4+γ (Q − )) ∩ C 1 (Q) and satisfies the bounds Also, the boundary layer function and satisfies the bounds [5, bound in (9)] As y, v, w and ψ i , i = 2, 3, 4 are all bounded, then the interior layer function z is also bounded.
Remark 1. We note that ifâ(s, t) = a(t), then z ≡ 0 and the coordinate transformation is not needed for this problem class. Error estimates and extensive numerical results for this problem class are given in [6].
In addition, for x = d, Proof. The interior layer function z is decomposed into the sum of two subcomponents where z c satisfies the problem and z p satisfies the problem In Appendix C, the subcomponent z p is further decomposed into the sum (36) Bounds on the derivatives of the subcomponent z q are also given in (38). Moreover, it is established in (35) and (40) that . From (2j) and (7a), note the following Using a comparison principle seperately on each subdomain Q − and Q + , we can then obtain the bounds Combining this bound with the bounds on z q (x, t) (from (38) in the final Appendix C) we achieve the pointwise bound in (11).
We transform the problems Lz c (x, t) =: and now apply the standard argument from [9, pg.352] , separately on the subdomains Q − and Q + , to deduce the remaining bounds.

3
Numerical method in the transformed domain and associated error analysis We approximate the solution of problem (6) on a rectangular grid in the  (10) and (11) on the layer components, this mesh is defined with respect to the transition points which split the interval [0, 1] into the five subdomains The grid points are uniformly distributed within each subinterval in the ratio 3N 8 : N 8 : N 8 : N 4 : N 8 . We discretize problem (6) using an Euler method to approximate the time variable and an upwind finite difference operator to approximate in space. Hence the discrete problem 3 is: Find Y such that Associated with this discrete problem is the upwinded finite difference operator: For any mesh function U , define 3 We use the following notation for various finite difference operators: This discrete operator satisfies a discrete comparison principle [5] and we can then establish that Hence Y QN,M ≤ C. To perform the error analysis the discrete solution is decomposed into the sum where V and W are the discrete counterparts to v and w. Using a standard argument [2] one can establish that For the remainder of the numerical analysis we will assume that ε is sufficiently small so that When this is not the case, the argument is classical as then ε −1 ≤ C ln N . The additional terms Ψ i , i = 2, 3, 4; and Z are defined as follows: For i = 2, 3, 4 and By the discrete comparison principle, we have that Ψ m ≤ C( √ ε) m , m = 2, 3, 4 and we can examine the truncation error for ψ m , m = 2, 3, 4: Applying the argument from [13] (see [6, Theorem 1] for more details) we deduce that where we have used the bounds established in Appendix A for the singular functions ψ m , m = 2, 3, 4. From the proof of Theorem 2, we have the bounds Also, as Z ≤ C, we can use a discrete comparison separately on each subinterval to sharpen the bound on Z(x i , t j ).
Proof. (a) For 0 ≤ x i ≤ d, consider the following barrier function The parameter θ ≥ 1 is specified below and M and N are sufficiently large so that 0 < c ≤ 1 − θT ln N M and ln N ≥ 1 + 1 T .
So, it follows that, when κ(x i , t j ) ≥ 0 and for N sufficiently large We need a modification to the argument if at any mesh point κ(x i , t j ) < 0. From (8), For the fine mesh points, where d − σ 1 ≤ x i < d, Then, by choosing θ ≥ 1 + A, we get that On the coarse mesh where 0 < x i < d − σ 1 , then using the inequality nt ≤ (1 + t) n , t ≥ 0, Then, for sufficiently large N and 0 < x i < d − σ 1 , Finish using a discrete comparison principle with the barrier function B(x i , t j )+ Ct j N −1 ln N . (b) For x i ≥ d, consider the following barrier function and we further assume that Note that if 1 − σ < x i < 1, then Hence, for sufficiently large N and all the mesh points where d < x i < 1, we repeat the argument from part (a) to conclude that Proof. From Theorem 2 and using e −θs 2 ≤ e 1 4θ e −s , we deduce that In addition, from Theorem 3 it also follows that Then, using the triangular inequality estimate (22) follows when . Hence we only now need to consider the error in the internal fine mesh. Within the fine mesh |κ(x i , t j )| ≤ Cσ 1 and so for Consider the piecewise linear barrier function, B(x i ) defined by and then we deduce the error bound using the discrete barrier fuction . and the discrete maximum principle.
The main result of this paper can now be stated.
Theorem 5. For sufficiently large N and M ≥ O(ln(N )), If Y is the solution of (18) and y is the solution of (6). Then, the global approximation Y onQ generated by the values of Y onQ N,M and bilinear interpolation, satisfies Proof. By combining the bounds in (19)

Modifications when source term is present
Here we outline the modifications to the method and to the analysis when b > 0. The problem is (2), but the differential equation (2a) is replaced with − εû ss +âû s +bû +û t =f , (s, t) ∈ (0, 1) × (0, T ]. In addition to all of the constraints imposed in (2), we also assume that b ∈ C 2+γ (Q),b ≥ 0 and the additional constraintb s (d, 0) =b ss (d, 0) = 0. As before, Γ * is defined by d (t) =â(d(t), t), d(0) = d. The operatorL d , given in (26), is redefined aŝ and we introduce a new function ThenL d (Iψ i ) = 0, i = 0, 1, 2, 3, 4. We redefine the function (3) to bê The changes in the transformed problem (6) are: Find y such that The discrete problem is defined as in (18). In the proof of Theorem 2, the presence of the source term will only effect the discussion of the regularity of the component z p (x, t) in Appendix C. In addition, the component z R is in the space C 4+γ (Q − ) ∪ C 4+γ (Q + ), due to the additional constraint imposed onb, and then the bounds (40) are also satisfied. Consequently, the proof of Theorem 5 will still apply.

Numerical results
In In [6] it is proved that the co-ordinate transformation (5) is not needed in order to obtain a global approximation whenâ only depends on the variable t. Hence, we first examine if this transformation is needed ifâ =â(s, t). In  Table 1, we see that, without the mapping, the method is not parameteruniform. Example 1 is now approximated with the numerical scheme (18) proposed in this paper. The computed approximations to y andû are displayed in Figure 1 and the maximum two-mesh global differences are given in Table 2. These numerical results are in agreement with Theorem 5.
Note that the source term is present in this example and then problem (25) is approximated with the numerical method (18) on the Shishkin meshQ N,M . For this example, we have I(t) = (cos t − 0.1 sin t)e −t 2 /2 . In addition, observe thatâ s (d, 0) = 0 andb s (d, 0) = 0. In Table 3 we see that the numerical approximations converge with almost first order.  [6] J.L. Gracia and E. O'Riordan, Numerical approximations to a singularly perturbed convection-diffusion problem with a discontinuous initial condition, Arxiv.

Appendix A: A set of singular functions
In this appendix the singular functionsψ i , i = 0, 1, 2, 3, 4 are defined and bounds of their derivatives are given. These functions are the main terms in the regularity expansion (33) of the continuous solutionû(s, t). These bounds are used in the truncation error analysis of the interior layer component z.
It follows that Observe that the bounds on the time derivatives of these two functions do not depend adversely on the singular perturbation parameter ε. This contrasts with the bounds on the time derivatives of these functions in the original variables (s, t).
In the transformed variables, we see from (27) that The fact that Lψ i = 0, whenâ depends on the spatial variable, results in the functionŷ exhibiting an interior layer (see Remark 1.) The next singular function is and the subsequent three functions 4 are As the first space derivatives of these functions are involved in the analysis of the interior layer function, we explicitly record that For these singular functions 5 , we can establish the bounds The functions ψ1, ψ2 were defined earlier by Shishkin in [10, (4.8c)] and Bobisud in [1] 5 For n = 0, 1, 2, 3, 4, ψn(x, 0) = 2(d − x) n , x > d; ψn(x, 0) = 0, x < d. and ∂ ∂t ψ n (x, t) ≤ C, ∂ 2 ∂x 2 ψ n (x, t) ≤ C; n = 2, 3, 4; on the second time derivatives on the fourth space derivatives and on the third space derivatives These expressions will be used to deduce bounds for the component z p in the decomposition (13) of z.

Appendix B: Decomposition of the solution
In this appendix we decompose the solution of problem (2) into a regularv, boundary layerŵ and interior layerẑ components. Bounds for the derivatives ofv andŵ are established here and the bounds for the component z in Appendix C.
We have the following expansion for the solution of problem (2): u(s, t) = 0.5 Note that the smooth remainderR satisfies the singularly perturbed problem The function p d and the related functionp d satisfy By the definitions (15) of z and (14) of the subcomponent z c we have where the values of K 4 and K 5 can be obtained from (37). Hence, the function z R is sufficiently regular within each sub-domain to allow us use results from [9] to bound the derivatives of z R . That is, Moreover, using | √ g − 1| ≤ |g − 1| and |g − 1| ≤ Ct, we conclude that, In the stretched variable we have that