1. Introduction

In this paper we consider numerical methods for solving elliptic and parabolic advection-diffusion-reaction (ADR) equations in one spatial dimension. We study the case in which the diffusion coefficient, advection coefficient, reaction coefficient and some terms are discontinuous across an interface. First, we consider a one-dimensional elliptic interface problem. Using the immersed interface method (IIM) ^[1] for finite difference approximations, then we demonstrate how to modify numerical methods constructed for homogeneous elliptic problems in the case in which the diffusion coefficient, advection coefficient, and reaction coefficient are discontinuous across the interface along with discontinuous sources. Next, we consider a parabolic problem in which the diffusion coefficient, advection coefficient, and reaction coefficient are discontinuous across the interface ^[3].

If numerical schemes that are used in case of the homogeneous diffusion, advection, and reaction coefficients are not modified around interfaces representing a discontinuity in the diffusion, advection, and reaction coefficients, the schemes lose their optimal order of accuracy. The IIM is a method that uses the undetermined coefficient method and correction terms to account for interfaces of discontinuity in the parameters of the partial differential equations (PDEs) which improves the accuracy of the schemes ^[4].

In section 2, we give an overview of the IIM. In section 3, we consider an elliptic interface problem with discontinuous diffusion, advection, and reaction coefficient. In this case, we modify a centered difference scheme at points near the interface (irregular points) by adding a correction term. The coefficients of the difference scheme are obtained using the undetermined coefficients method. The centered difference scheme is second order accurate, but it loses the optimal order at irregular points where it is first order accurate.

In section 4, we consider a parabolic interface problem. We show how to modify a scheme which is centered difference in space and Crank-Nicolson in time to obtain a second order scheme in time and space.

2. Numerical Methods for Interface Problems

We consider one dimensional advection diffusion reaction equations (elliptic and time independent problems) with discontinuities in the coefficients. Discontinuities in the solutions of these PDEs occur when, for example, when there are different materials that are separated by sharp interfaces. Across such interfaces there can be sharp discontinuities in the diffusion coefficient as well as the advection and reaction terms.

By an interface problem, we mean an initial boundary value problem in which the following conditions are satisfied ^[4].

•  There can be one or several fixed interfaces in the solution domain.

•  The diffusion coefficients as well as advection velocity and/or reaction terms could be discontinuous across these fixed interfaces.

•  The source term may have a finite jump along interfaces.

•  The solution to an interface problem may be discontinuous, or its gradient derivative is discontinuous across the interface.

•  We have some prior knowledge of the jump conditions or interface conditions of the solution and ux across interfaces.

Numerical methods constructed for initial boundary value problems involving constant coefficients work poorly when applied without modification to problems with interfaces. There are several different techniques available in the literature to construct numerical methods for interface problems. Here, we present a technique called the Immersed Interface Method (IIM), which is a technique that modifies standard numerical methods for problems without interfaces to capture the discontinuities in the solution and ux. As outlined in ^[4], a method is called an IIM if it satisfies the following properties:

1. Uniform Grids: The is method uses a uniform grid in Cartesian, polar or spherical geometries instead of a body fitted grid as in the finite element method.

2. Jump Conditions: The interior boundary conditions at interfaces which are known from the PDEs are the jump conditions or interface conditions. More interface relations can be derived from given jump conditions.

3. Standard Schemes at Regular Points: Standard numerical methods are used to discretize the scheme at grid points that are far away from the interface (regular points), but the numerical methods are modified by correction terms at the points near the interface (irregular points). We obtain these correction terms for modifying the standard numerical schemes through the derivation of the truncation error at each irregular point using the given jump conditions.

4. Pointwise Convergence: We will be concerned with the pointwise convergence of the resulting immersed interface methods through the computation of the L₁ error instead of the average L₁ or L₂ error. This is because, the L₁ error will take into account the error near the interface, i.e. at the irregular points. The IIM generally has the same global order of accuracy as when the method is applied to regular problems without interfaces of discontinuity. However, the local truncation error may have a lower order of accuracy at irregular grid points than at regular grid points.

We will first consider an elliptic interface problem and then consider the time dependent case (parabolic problem).

3. Elliptic Interface Problem

We begin with a one-dimensional elliptic problem with a discontinuity at x =ξ which falls between irregular points x_j and x_j+1 for some j ∈ N so that ξ ∈ [x_j, x_j+1] (see Figure 1). The elliptic interface problem that we consider here is finding the solution u = u(x) of the initial-boundary value problem in domain Ω = [a, b] such that:

(3.1a)

(3.1b)

(3.1c)

(3.1d)

Figure 1. Diagram of ID-interface problem with one interface ξ

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View next figure

(3.1e)

(3.1f)

In the above equations, the coefficients β > 0,σ,α and the source term f are known a-priori. We assume that α and σ are piecewise constant throughout the domain and (3.1b) and (3.1c) are the interior boundary conditions or jump conditions at the fixed interface at x = ξ. For any function v, the notation and . Boundary conditions specifying the solution at the boundaries of are given in (3.1d). There is a jump at x = ξ in the coefficients β,α,σ and the source term f as well as the solution u. However, away from the interface β,α,σ and f are smooth.

Below we outline the procedure for constructing an IIM scheme for the interface problem (3.1).

Step 1: Construction of a Uniform Cartesian Grid

We first construct an uniform Cartesian grid. The discrete mesh h consists of grid points

Here h > 0 is the mesh step size. We also define half-points . The conventional finite difference schemes can be used on most grid points except at the irregular points. In Figure 1, x_j and x_j+1 are irregular points.

Step 2: Derivation of a Centered Finite Difference Scheme

Next, we derive the centered finite difference scheme for equation (3.1a). Since the solution is smooth away from the interface points, we use standard finite approximations at regular grid points.

Scheme at Regular Grid Points:

Using a second order centered difference scheme to discretize the spatial derivatives at regular points x_i, we have

where and

Thus, the fully discrete scheme at a regular grid point x = x_i is:

(3.2)

where u_i ≈u(x_i), and f_i = f(x_i).

At regular grid points away from the interface, scheme (3.2) is second order accurate in space. We define for 1 ≤ i ≤ N-1, i≠j, j + 1 the parameters:

(3.3a)

(3.3b)

(3.3c)

with definition (3.3), scheme (3.2) the scheme can be written as:

(3.4)

for 1 ≤ i ≤ N-1, i≠j, j + 1.

Scheme at Irregular Grid Points:

To determine the finite difference equations at irregular grid points we use the method of undetermined coefficients ^[2]. To keep a unified notation, we assume that the scheme at the irregular points x_j and x_j+1 is of the form

(3.5a)

(3.5b)

where we add correction terms C_j, and C_j+1 on the right hand side of the equations in (3.5a) and (3.5b) respectively in order to let the scheme achieve the optimal order of accuracy at irregular points, x_j and x_j+1. The correction terms are obtained by constructing the Local Truncation Error (LTE) at the irregular points, using Taylor expansions of the solutions in which the interior boundary conditions (3.1b) and (3.1c) are used. We will derive the coefficients in the next step as well as the correction terms C_j and C_j+1.

Step 3: Solution of a Linear System

We can write the scheme at all the grid points (regular and irregular) in vector form as

where

and

is a tridiagonal matrix. The solution of this linear system gives us the approximate solution to the elliptic interface problem (3.1).

We compute the jump condition for in order to use it to derive the local truncation error as well as obtain the correction terms. Consider (3.1a) at interface x = ξ and substitute the interior boundary conditions. We note that the coefficients σ and α are constant and the source term f has a jump at the interface x = ξ. Thus,

By definition and from the jump conditions (3.1b), (3.1c) we have

(3.6)

Substituting the jump conditions (3.1b) and (3.1c) in (3.6) we have

(3.7)

Moreover, from (3.6), (3.1b) and (3.1c) we have in terms of as

(3.8)

In the above β,α,σ and the source term f are the only discontinuous coefficients.

The coefficients at regular points are given in (3.3). In this step, we aim to find the coefficients where k = 1,2,3 in the scheme at irregular points. At regular points the solution is smooth and the truncation error is O(h²) since the centered difference scheme has second order accuracy, but the smoothness is poor near irregular grid points. We derive the local truncation error at irregular points x_j and x_j+1 and we expect to lose the optimal order of this error. Let the local truncation error at x = xi be denoted as T_i, 1 ≤ i ≤ N-1.

The local truncation error at x = x_j:

The scheme at the irregular point x_j is

(3.9)

Since x_j-1, x_j fall on the left of the interface x = ξ and x_j+1 falls on the right of interface (see Figure 1) we compute Taylor expansions of the solution u(x) of the scheme (3.9) around the interface x = ξ and rewrite terms in u^-(ξ). The local truncation error T_j is defined as

(3.10)

Taylor expansions for u(x_j-1); u(x_j) and u(x_j+1) around x = ξ, with respect to the interior conditions (3.1b), (3.1c), the PDE (3.1a) and the additional jump condition defined in (3.7) are computed as:

(3.11a)

(3.11b)

(3.11c)

We define

(3.12)

Substituting (3.11a), (3.11b), and (3.11c) in (3.10), using the PDE (3.1a) we have

(3.13)

Comparing (3.13) with the PDE (3.1a), which in approaching ξ from the left can be written as

We obtain the following system of linear equations in the coefficients

(3.14a)

(3.14b)

(3.14c)

We define

(3.15)

The local truncation error at x = x_j+1:

The scheme at the irregular point x_j+1 is

(3.16)

Since x_j falls on the left of the interface x = ξ and x_j+1, x_j+2 fall on the right of interface (see Figure 1) we compute Taylor expansions of the solution u(x) in the scheme (3.16) around the interface x = ξ and rewrite terms in u⁺(ξ). The local truncation error T_j+1 is defined as

(3.17)

Taylor expansion for u(x_j), u(x_j+1) and u(x_j+2) around x = ξ with respect to the interior conditions (3.1b), (3.1c), the PDE (3.1a) and the additional jump condition as defined in (3.8) are computed as

(3.18a)

(3.18b)

(3.18c)

We define h₄ = (x_j+2- ξ) and substitute (3.18a), (3.18b), and (3.18c) in (3.17) to get

(3.19)

When we compare (3.19) with the PDE (3.1a), which in approaching ξ from the right is given as

then, we obtain the following linear system in the coefficients and :

(3.20a)

(3.20b)

We define

The system of linear equations (3.20), (3.20b), (3.20b) can be solved easily using a computer program like MAPLE to obtain the coefficients of the centered finite difference scheme at the irregular points x_j and x_j+1.

We will measure the error using the L∞ norm defined as

(3.21)

where u_l is numerical solution and u(x_l) is the exact solution at grid point x_l.

Order of convergence is , where is the error corresponding to step size h_i, i = 1, 2.

Example 3.1: Consider the initial-boundary value problem with interface ξ = 0.5, and the jump conditions: [u] = 1/3 and [α-βu_x] =-1. In this problem σ = 0 and β, α, and the source term are given as

(3.22)

The true solution for this problem is

In Figure 2, we plot the numerical and exact solution for Example 3.1 with N=80 and h=1/80. Table 1 shows that the error converges as O(h²).

Table 1. Table of errors for Example 3.1

Download as

PowerPoint Slide

Larger image(png format)

Tables index

Veiw figure

View current table in a new window

View next table

In Figure 2, we plot the numerical solution and the exact solution for Example 3.1 with N=80 and h=1/80. Table 1 shows that the IIM reaches the optimal order of accuracy (i.e second order) in the case where the diffusion, advection, and reaction coefficient are discontinuous across the interfaces.

Figure 2. The plot of the numerical solution (the \o"-line) and the exact solution (the solid line) for Example

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

View next figure

Figure 3. The plot of Error for Example 3.1 with N=80, h=1/80, k=0.1 and t=0.5

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

View next figure

4. Time-Dependent Interface Problem

In this section, we will construct an IIM for the time-dependent (Parabolic) interface problem. Consider the following time-dependent initial-boundary value problem in one dimension with time dependent source term f and assume that the coefficients σ and α are piecewise constant. Find the solution u = u(x,t) on Ω(0,T] of the initial-boundary value problem:

(4.1a)

(4.1b)

(4.1c)

(4.1d)

(4.1e)

(4.1f)

(4.1g)

(4.1h)

In the above equations, the coefficients β > 0, σ, and α, the source term f, the jump conditions w and v, and the initial and boundary functions u₀(x) and g(x) are known a-priori. We assume β, α, σ, u, and f have a finite (time-dependent) discontinuity at interface x = ξ. We now present the construction of the IIM for problem (4.1).

Step 1: Construction of the Uniform Cartesian Grid for Space-Time Mesh

The discrete mesh consists of grid points (x_i, tⁿ) defined as

where h is the space step size, and k is the time step. We also define half-points as well as . Conventional finite difference schemes can be used on most grid points except at (see Figure 4), where x_j, x_j+1 are the spatial irregular points.

Figure 4. Stencil for CFD-CN scheme (4.8)

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

View next figure

Step 2: Derivation of a Centered Finite Difference Scheme in Space and Crank-Nicolson

Scheme in Time at Regular Grid Points We will derive the scheme for (4.1a) which is centered finite difference (CFD) in space and Crank-Nicolson (CN) in time. This scheme is unconditionally stable when applied to regular grid points if σ > 0 regardless of the jump, provided that β(x,t) has the same sign across the interface. Then the scheme discretized at is derived as he following:

In the above, . We define for 1 ≤ i ≤ N-1 the following parameters:

(4.2)

(4.3)

(4.4)

The coefficients (4.2), (4.3), and (4.4) correspond to the scheme at regular points. Thus, the standard Crank-Nicolson and centered spatial finite difference scheme at regular grid points can be written as:

where and .

Scheme at Irregular Grid Points:

To determine the finite difference equations at irregular grid points we use the method of undetermined coefficients ^[2]. We assume that the scheme at the irregular points and is of the form

(4.5)

(4.6)

where we add correction terms , , and , on the right hand side of the equations in (4.5) and (4.6) respectively in order to let the scheme reach the optimal order of accuracy (i.e second order) at irregular points x_j and x_j+1 for each time point tⁿ. We will derive the coefficients for k = 1, 2, 3; and l = j, j +1 in the next step. We also derive the correction terms and through computing the Local Truncation Error (LTE) at the irregular points, using Taylor expansions of the solutions at each time point tⁿ in which the interior boundary conditions (4.1b) and (4.1c) are used.

Step 3: Solution of a Linear System

This finite difference scheme for problem (4.1) on can be written in vector form as

(4.7)

where

with

where if for all

The solution of the linear system (4.7) of equations is the numerical solution

In this subsection we want to compute The source term has a jump across the interface at each time point

Then, we have

(4.8)

We define and use the given jump conditions (4.1b) and (4.1c) to get

(4.9)

We can also rewrite in terms of as

(4.10)

The coefficients of the scheme at regular points are given in (4.2), (4.3), and (4.4). The coefficients of the scheme at irregular points are derived by substituting Taylor expansions of the solution around interface in the local truncation error at and by using the jump condition (4.1b), (4.1c), and (4.9) on . The local truncation error and is slightly different than and computed in section 2 because of extra term .

The local truncation error at is

(4.11)

The scheme is discretized at as the following:

(4.12a)

(4.12b)

The Taylor expansions for and around are the same as above, but in terms of Also, we use he following expansion:

to obtain

(4.13)

where

(4.14)

The PDE that approaches from the left of the interface at time point is

(4.15)

where we substitute in (4.11) and compare the factors of , and we obtain the following linear system

The correction term corresponding to is given by

(4.16)

Also, we obtain the same linear system of equations and correction term for the time level n to compute and

The local truncation error at :

(4.17)

We use Taylor expansion of the solution around interface to obtain

(4.18a)

(4.18b)

(4.18c)

We will define The Taylor expansion for around is the same as above replacing by Also, we have

The local truncation error corresponding to is

(4.19)

The PDE that approaches from the right of the interface on time point is

(4.20)

From (4.19) and (4.20), we obtained the following linear system for which that the solutions are the coefficients of the scheme corresponding to :

The correction term corresponding to (x_j+1, tⁿ⁺¹) is

(4.21)

The linear system for computing the coefficients of the scheme at (x_j+1, tⁿ) and correction term are the same as in time level n + 1. These linear systems of equations can be solved using a computer program such as MAPLE.

Example 4.1: We will estimate the errors using the norm L∞ in space mesh at all grid point as

Also, the errors are estimated in time mesh using

(4.22)

Consider the initial-boundary value problem with interface: x = 0.5. and the jump conditions: [u] =1/2cos(t) as well as [α - βu_x] =1/3cos(t) Let T=[0,1]. Then

(4.23)

The true solution is

Table 2. Table of errors for Example 4.1 with k/h=1

Download as

PowerPoint Slide

Larger image(png format)

Tables index

Veiw figure

View current table in a new window

View previous table

View next table

Figure 5. The plot of the numerical solution (the \o"-line) and the exact solution (the solid line) for Example 4.1 with N=80, h=1/80, k=0.1 and t=0.5

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

View next figure

Figure 6. The plot of Error for Example 4.1 with N=80, h=1/80, k=0.1 and t=0.5

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

View next figure

In Figure 5, we plot the numerical solution for N = 80, t = 0.5, h = 1/80, k = 0.1. Table 2 shows the error converges to O(h²). Table 2 shows that the error converges to O(k²). Therefore, the truncation error is second order of accuracy in space as well as in time.

Example 4.2: Consider the initial-boundary value problem with interface: x = 0.5. and the jump conditions: [u] =1/2cos(t) as well as [α - βu_x] =1/3cos(t). Let T=[0,1]. then

(4.24)

The true solution is

Table 3. Table of errors for Example 4.1 with k/h=1

Download as

PowerPoint Slide

Larger image(png format)

Tables index

Veiw figure

View current table in a new window

View previous table

View next table

Figure 7. The plot of the numerical solution (the \o"-line) and the exact solution (the solid line) for Example 4.2 with N=80, h=1/80, k=0.1 and t=0.5

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

View next figure

In Figure 7, we plot the numerical solution for N = 80, t = 0.5, h = 1/80, k = 0.1. Table 3 shows the error converges to O(h²). Table 3 shows that the error converges to O(k²). Therefore, the truncation error is second order of accuracy in space as well as in time

Example 4.3: Consider the initial-boundary value problem with interface: x = 0.5. and the jump conditions: [u] =1/2cos(t) as well as [α - βu_x] =1/3cos(t) Let T=[0,1]. then

(4.25)

The true solution is

Figure 8. The plot of Error for Example 4.2 with N=80, h=1/80, k=0.1 and t=0.5

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

View next figure

In Figure 9, we plot the numerical solution for N = 80, t = 0.5, h = 1/80, k = 0.1. Table 4 shows the error converges to O(h²). Table 4 shows that the error converges to O(k²). Therefore, the truncation error has the second order of accuracy in space as well as in time.

Table 4. Table of errors for Example 4.1 with k/h=1

Download as

PowerPoint Slide

Larger image(png format)

Tables index

Veiw figure

View current table in a new window

View previous table

Figure 9. The plot of the numerical solution (the \o"-line) and the exact solution (the solid line) for Example 4.3 with N=80, h=1/80, k=0.1 and t=0.5

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

View next figure

Figure 10. The plot of Error for Example 4.3 with N=80, h=1/80, k=0.1 and t=0.5

Download as

PowerPoint Slide

Larger image(png format)

Figures index

Veiw figure

View current figure in a new window

View previous figure

5. Conclusions

In this paper, we have presented the Immersed Interface Method (IIM) to solve one-dimensional elliptic and parabolic interface problems with discontinuous diffusion, advection, and reaction coefficients. We presented numerical examples that compute the approximate solution across a fixed interface and demonstrate how to modify a centered difference scheme at irregular points by adding the correction terms. The correction terms were obtained by deriving the truncation error and using the given jump conditions. The method of undetermined coefficients was used to obtain the (modified) coefficients of the scheme at irregular points.

The result is that we obtained the second order of accuracy for the scheme modified by IIM.

References

[1]	X. Feng and Z. Li, Simplified Immersed Interface Methods for Elliptic Interface Problems with Straight Interfaces, Numerical Methods for Partial Differential Equations, 28 (2012), pp. 188-203.
	In article		CrossRef
[2]	R. J. LeVeque, Finite difference methods for ordinary and partial differential equations: Steady-state and time-dependent problems, (2007).
	In article
[3]	Zhilin Li, The immersed Interface method: A numerical approach for partial differential equation with interface, PhD thesis, Univerdity of Washington, 1994.
	In article
[4]	Z. Li and K. Ito, The immersed interface method: numerical solutions of PDEs involving interfaces and irregular domains, vol. 33, Society for Industrial Mathematics, 2006.
	In article		CrossRef