Fourth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 03, 1999, pp 100-108. http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)

A numerical method for solving heat equations involving interfaces ^{∗}

Zhilin Li & Yun-Qiu Shen

Abstract

In 1993, Li and Mayo [3] gave a finite-difference method with second order accuracy for solving the heat equations involving interfaces with constant coefficients and discontinuous sources. In this paper, we expand their result by presenting a finite-difference method which allows each coefficient to take different values in different sub-regions of the interface. Our method is useful in physical applications, and has also second order accuracy.

1 Introduction

Consider the heat equation

u_{t} = (βu_{x})_{x} + (βu_{y})_{y} + f(x, y, t) ,  (1)

where t is in [a, ∞) and (x, y) in an open region Ω ⊂ R^{2} which is divided into two sub-regions Ω^{+}, Ω^{−} by an irregular interface Γ.

Let Ω = Ω^{+} ∪ Γ ∪ Ω^{−} be an open rectangular region in R^{2}, and let

β(x, y) = \{ \begin{matrix}β^{+}, & if (x, y) ∈ Ω^{+}, \\ β^{−}, & if (x, y) ∈ Ω^{−}, \end{matrix}

where β^{+} and β^{−} are positive constants. Assume that f is continuous in each sub-region, and may have jumps of discontinuty across Γ. The solution u(x, y, t) and its normal derivative u_{n}(x, y, t) crossing the curve Γ have prescribed jumps

[u] ≡ u^{+}(x(s), y(s), t) − u^{−}(x(s), y(s), t) = ω(s, t),  (2)
[u_{n}] ≡ u_{n}^{+}(x(s), y(s), t) − u_{n}^{−}(x(s), y(s), t) = g(s, t),  (3)

where s is a parameter of Γ, the superscripts + and − denote the limiting values of a function from one side in Ω^{+} and another side in Ω^{−} respectively. Throughout this paper, we use [G] = G^{+} − G^{−} to denote the difference of the limiting values of a function G from different sides of the interface as in (2), (3). We also assume that the initial condition at t = a and the boundary condition on ∂B are known. In 1993, Li and Mayo [3] gave a finite-difference method with second accuracy for solving (1), assuming that f is discontinuous but β is constant. In this paper, we present a finite-difference method for solving (1), which allows β to be taken different values in different sub-regions of the interface. This feature is useful in physical applications.

∗Mathematics Subject Classifications: 65M06, 65M12, 65M15.

Key words: Finite difference, Second order accuracy, Interface, Local coordinates.

\begin{matrix}c \\ ◯ \end{matrix}2000 Southwest Texas State University and University of North Texas.

Published July 10, 2000.

We organize this paper as follows: In Section 2, we establish local coordinate systems around the interface. In Section 3, we give the correction terms for the finite-difference method. In Section 4, we show that our method is second order accurate. Finally, in Section 5, we give some numerical examples, in which the actual solutions are known, to confirm the theoretical result.

2 Local Coordinate Systems

We first give local coordinate systems along the interface Γ, as in [2, 3]. When a point (x_{0}, y_{0}) on the interface is fixed for the origin, we use the normal direction as the ξ-direction, which has an angle θ with the x-axis. Rotating the ξ-direction by ninety degrees counter-clockwise, we obtain the η-direction. Now we express the curve Γ as a function of the independent variable η locally.

ξ = χ(η).  (4)

We express (2),(3) locally by using the η coordinate.

[u] ≡ u^{+} − u^{−} = ω(η, t),  (5)
[u_{n}] ≡ u_{n}^{+} − u_{n}^{−} = g(η, t),  (6)

where ω and g are known in advance. From (4) and (6), we have

g = [u_{n}] = \frac{−[u_{η}]χ_{η} + [u_{ξ}]}{\sqrt{1 + χ_{η}^{2}}}.

Using differentiation with respect to η, noting χ_{η}(0) = 0, we have

[u_{η}] = ω_{η},  (7)
[u_{ξ}] = g,  (8)
[u_{ηη}] = −gχ_{ηη} + ω_{ηη},  (9)
[u_{ξη}] = ω_{η}χ_{ηη} + g_{η}.  (10)

Changing (1) locally, using

ξ = (x − x_{0})cosθ + (y − y_{0})sinθ,  (11)
η = −(x − x_{0})sinθ + (y − y_{0})cosθ,  (12)

we have

u_{xx} + u_{yy} = u_{ξξ} + u_{ηη},

therefore,

u_{t} = βu_{ξξ} + βu_{ηη} + f̃(ξ, η, t),  (13)

where f̃(ξ, η, t) = f(x(ξ, η), y(ξ, η), t). (13) implies

[u_{ξξ}] = [u_{ξξ}]_{1} + [\frac{u_{t}}{β} ] ,  (14)

where

[u_{ξξ}]_{1} = gχ_{ηη} − ω_{ηη} − [\frac{f̃}{β}] .  (15)

In (7)-(10)and (14)-(15), all functions in the right-hand sides are known except for [u_{t}/β] in (14) which will be explored in the next section.

3 Correction Terms

We first discretize in both of the x-direction and the y-direction with a mesh of size h.

u_{xx} ≈ δ_{x}u_{i,j} ≡ \frac{1}{h^{2}} (u_{i−1,j} − 2u_{i,j} + u_{i+1,j}),  (16)
u_{yy} ≈ δ_{y}u_{i,j} ≡ \frac{1}{h^{2}} (u_{i,j−1} − 2u_{i,j} + u_{i,j+1}).  (17)

We group all the grid points in Ω into two sets. The set S_{reg} consists the regular points, each point in one sub-region has no neighbor point in the another subregion. The set S_{irr} consists the irregular points, each point in one sub-region has at least one neighbor point in the other sub-region. For a regular grid point, the local truncation error of (16),(17) from βu_{xx} + β_{yy} + f is O(h^{2}).

For an irregular grid point, we need add some correction terms in (16),(17) such that the local truncation error is O(h), therefore the global error of the solution for solving the heat equation is O(h^{2}) after the discretization of time t in certain way. At first, we relate the jumps with respect to x and y to the jumps with respect to ξ and η by (11) and (12).

[u_{x}] = [u_{ξ}] cos θ − [u_{η}] sin θ,
[u_{y}] = [u_{ξ}] sin θ + [u_{η}] cos θ,  (18)
[u_{xx}] = [u_{xx}]_{1} + [\frac{u_{t}}{β} ] cos^{2} θ,

where

[u_{xx}]_{1} = [u_{ξξ}]_{1} cos^{2} θ − 2[u_{ξη}] cos θ sin θ + [u_{ηη}] sin^{2} θ,
[u_{yy}] = [u_{yy}]_{1} + [\frac{u_{t}}{β} ] sin^{2} θ,  (19)

where

[u_{yy}]_{1} = [u_{ξξ}]_{1} sin^{2} θ + 2[u_{ξη}] cos θ sin θ + [u_{ηη}] cos^{2} θ.

By (7)-(10) and (14)-(15), all the terms in the right-hand sides of (18)-(19) are known except for [u_{t}/β].

Now we consider an irregular grid point in the x-direction, there are four cases:

(a) (x_{i}, y_{j}) ∈ Ω^{−}, (x_{i+1}, y_{j}) ∈ Ω^{+}

(b) (x_{i}, y_{j}) ∈ Ω^{+}, (x_{i−1}, y_{j}) ∈ Ω^{−}

(c) (x_{i}, y_{j}) ∈ Ω^{+}, (x_{i+1}, y_{j}) ∈ Ω^{−}

(d) (x_{i}, y_{j}) ∈ Ω^{−}, (x_{i−1}, y_{j}) ∈ Ω^{+}

For case (a), let the intersection of the line segment connecting (x_{i}, y_{j}), (x_{i+1}, y_{j}) and Γ be (x^{∗}, y_{j}). Using Taylor’s series around x^{∗}, we have

\frac{1}{h^{2}} (u(x_{i−1}, y_{j}) − 2u(x_{i}, y_{j}) + u(x_{i+1}, y_{j}))

= \frac{1}{h^{2}} \{[u] + [u_{x}](x_{i+1} − x^{∗}) + \frac{[u_{xx}]}{2!}(x_{i+1} − x^{∗})^{2}\} + u_{xx}^{−} + O(h),

which implies

u_{xx} = δ_{x}u_{i,j} + C_{x}u_{i,j} − \frac{[\frac{u_{t}}{β}](x_{i+1} − x^{∗})^{2}}{2h^{2}} + O(h),  (20)

where

C_{x}u_{i,j} = −\frac{1}{h^{2}} \{[u] + [u_{x}](x_{i+1} − x^{∗}) + \frac{[u_{xx}]_{1}}{2!} (x_{i+1} − x^{∗})^{2}\} .

Similarly, for case (b), we have

u_{xx} = δ_{x}u_{i,j} + C_{x}u_{i,j} + \frac{[\frac{u_{t}}{β}](x_{i−1} − x^{∗})^{2}}{2h^{2}} + O(h),  (21)

where

C_{x}u_{i,j} = \frac{1}{h^{2}} \{[u] + [u_{x}](x_{i−1} − x^{∗}) + \frac{[u_{xx}]_{1}}{2!} (x_{i−1} − x^{∗})^{2}\} .

For case (c), we have

u_{xx} = δ_{x}u_{i,j} + C_{x}u_{i,j} + \frac{[\frac{u_{t}}{β}](x_{i+1} − x^{∗})^{2}}{2h^{2}} + O(h),  (22)

where

C_{x}u_{i,j} = \frac{1}{h^{2}} \{[u] + [u_{x}](x_{i+1} − x^{∗}) + \frac{[u_{xx}]_{1}}{2!} (x_{i+1} − x^{∗})^{2}\} .

For case (d), we have

u_{xx} = δ_{x}u_{i,j} + C_{x}u_{i,j} − \frac{[\frac{u_{t}}{β}](x_{i−1} − x^{∗})^{2}}{2h^{2}} + O(h),  (23)

where

C_{x}u_{i,j} = −\frac{1}{h^{2}} \{[u] + [u_{x}](x_{i−1} − x^{∗}) + \frac{[u_{xx}]_{1}}{2!} (x_{i−1} − x^{∗})^{2}\} .

Analogously, in the y-direction, we also have four cases. We add correction terms such that the local truncation is O(h).

Now using these correction terms in both x-and y-directions, we obtain a system of ordinary differential equations

(u_{i,j})_{t} = β(δ_{x}u_{i,j} + δ_{y}u_{i,j}) + f_{i,j}
+β ∑_{(x_{i_{0}} ,y^{∗})∈S_{irr}} \{C_{x}u_{i,j} + [\frac{u_{t}}{β} ] \frac{τ_{x_{0}}(x_{i_{0}} − x^{∗})^{2} cos^{2} θ}{2h^{2}} \}  (24)

+β ∑_{(x^{∗},y_{j_{0}} )∈S_{irr}} \{C_{y}u_{i,j} + [\frac{u_{t}}{β} ] \frac{τ_{y_{0}}(y_{j_{0}} − y^{∗})^{2} sin^{2} θ}{2h^{2}} \} + O(h),

where i_{0} = i − 1 or i + 1, j_{0} = j − 1 or j + 1. τ_{x_{0}}, τ_{y_{0}} = 1 or −1 according to (20)-(23). We have

[\frac{u_{t}}{β} ] = u_{t}^{−} [\frac{1}{β} ] + \frac{[u_{t}]}{β^{+}} = u_{t}^{+} [\frac{1}{β} ] + \frac{[u_{t}]}{β^{−}},

which imply

[\frac{u_{t}}{β} ] = (u_{i,j})_{t} [\frac{1}{β} ] + \frac{ω_{t}}{β̃} + O(h),  (25)

where β̃ = β^{−} if β = β^{+}, β̃ = β^{+} if β = β^{−}, and ω is defined in (5). Using (25), we have the following system of ordinary differential equations

(u_{i,j})_{t} = F (u_{i,j−1}, u_{i−1,j}, u_{i,j}, u_{i+1,j}, u_{i,j+1})

with the right-hand side F is equal to

\frac{β(δ_{x}u_{i,j} + δ_{y}u_{i,j} + ∑_{(x_{i_{0}} ,y^{∗})∈S_{irr}} C̃_{x}u_{i,j} + ∑_{(x^{∗},y_{j_{0}} )∈S_{irr}} C̃_{y}u_{i,j}) + f_{i,j}}{1 − ∑_{(x_{i_{0}} ,y^{∗})∈S_{irr}} β[\frac{1}{β}]\frac{τ_{x_{i_{0}}} (x_{i_{0}} −x^{∗})^{2} cos^{2} θ}{2h^{2}} − ∑_{(x^{∗},y_{j_{0}} )∈S_{irr}} β[\frac{1}{β}]\frac{τ_{y_{j_{0}}} (y_{j_{0}} −y^{∗})^{2} sin^{2} θ}{2h^{2}}} ,

(26)

where

C̃_{x}u_{ij} = C_{x}u_{ij} + \frac{ω_{t}}{β̃} \frac{τ_{x_{i_{0}}} (x_{i_{0}} − x^{∗})^{2} cos^{2} θ}{2h^{2}}

and

C̃_{y}u_{ij} = C_{y}u_{ij} + \frac{ω_{t}}{β̃} \frac{τ_{y_{j_{0}}} (y_{j_{0}} − y^{∗})^{2} sin^{2} θ}{2h^{2}} .

At a regular grid point, the local truncation error of the right-hand side of (26) from βu_{xx} + β_{yy} + f is O(h^{2}). In the next section, we will show that at an irregular grid point, the local truncation error is O(h).

Finally, we discretize time t by choosing Δt = h. We use Crank-Nicholson method.

u_{i,j,k+1} = u_{i,j,k} + Δt(0.5F (u_{i,j−1,k}, u_{i−1,j,k}, u_{i,j,k}, u_{i+1,j,k}, u_{i,j+1,k}) (27)
+0.5F (u_{i,j−1,k+1}, u_{i−1,j,k+1}, u_{i,j,k+1}, u_{i+1,j,k+1}, u_{i,j+1,k+1})),

which implies the local truncation error for discretizing on t is O((Δt)^{2}), [1, 4]. To solve for u_{i,j,k+1}, from (27), we use S.O.R. iteration with a certain parameter λ. Set

u_{i,j,k+1}^{(0)} = u_{i,j,k},

and for n = 1, 2, …, set

u_{i,j,k+1}^{(n+1)} = (1 − λ)u_{i,j,k+1}^{(n)} + λ(u_{i,j,k}
+0.5F (u_{i,j−1,k}, u_{i−1,j,k}, u_{i,j,k}, u_{i+1,j,k}, u_{i,j+1,k})  (28)

+0.5F (u_{i,j−1,k+1}^{(n+1)}, u_{i−1,j,k+1}^{(n+1)}, u_{i,j,k+1}^{(n)}, u_{i+1,j,k+1}^{(n)}, u_{i,j+1,k+1}^{(n)})) .

4 Accuracy Analysis

We first show that the denominator of the right-hand side of (26) is bounded below and above by positive constants. Let

D_{i,j} ≡ 1 − ∑_{(x_{i_{0}} ,y^{∗})∈S_{irr}} β [\frac{1}{β} ] \frac{τ_{x_{i_{0}}} (x_{i_{0}} − x^{∗})^{2} cos^{2} θ}{2h^{2}}  (29)

− ∑_{(x^{∗},y_{j_{0}} )∈S_{irr}} β [\frac{1}{β} ] \frac{τ_{y_{j_{0}}} (y_{j_{0}} − y^{∗})^{2} sin^{2} θ}{2h^{2}} .

and

E_{i,j} ≡ 1 + ∑_{(x_{i_{0}} ,y^{∗})∈S_{irr}} |β [\frac{1}{β} ] \frac{τ_{x_{i_{0}}} (x_{i_{0}} − x^{∗})^{2} cos^{2} θ}{2h^{2}} |  (30)

+ ∑_{(x^{∗},y_{j_{0}} )∈S_{irr}} |β [\frac{1}{β} ] \frac{τ_{y_{j_{0}}} (y_{j_{0}} − y^{∗})^{2} sin^{2} θ}{2h^{2}} | .

Lemma 4.1 Let D_{i,j} and E_{i,j} be as defined above. Then

D_{i,j} ≥ \frac{min(β^{+}, β^{−})}{max(β^{+}, β^{−})} and E_{i,j} ≤ 1 + max(β^{+}, β^{−}) |[\frac{1}{β} ]| ,  (31)

where i_{0} = i − 1 or i + 1, j_{0} = j − 1 or j + 1. τ_{x_{i_{0}}} , τ_{y_{j_{0}}} = 1 or −1 according to (20)-(23).

Proof: At first we prove the lower bound. Look at the first summation. In x-direction, we have four cases.

(a) (x_{i}, y_{j}) ∈ Ω^{−}, (x_{i+1}, y_{j}) ∈ Ω^{+}. Then x_{i_{0}} = x_{i+1}, τ_{x_{i_{0}}} = −1, β = β^{−}, and

β [\frac{1}{β} ] τ_{x_{i_{0}}} = 1 − \frac{β^{−}}{β^{+}} which is positive iff β^{−} < β^{+}.

(b) (x_{i}, y_{j}) ∈ Ω^{+}, (x_{i−1}, y_{j}) ∈ Ω^{−}. Then x_{i_{0}} = x_{i−1}, τ_{x_{i_{0}}} = 1, β = β^{+}, and

β [\frac{1}{β} ] τ_{x_{i_{0}}} = 1 − \frac{β^{+}}{β^{−}} which is positive iff β^{+} < β^{−}.

(c) (x_{i}, y_{j}) ∈ Ω^{+}, (x_{i+1}, y_{j}) ∈ Ω^{−}. Then x_{i_{0}} = x_{i+1}, τ_{x_{i_{0}}} = 1, β = β^{+}, and

β [\frac{1}{β} ] τ_{x_{i_{0}}} = 1 − \frac{β^{+}}{β^{−}} which is positive iff β^{+} < β^{−}.

(d) (x_{i}, y_{j}) ∈ Ω^{−}, (x_{i+1}, y_{j}) ∈ Ω^{+}. Then x_{i_{0}} = x_{i−1}, τ_{x_{i_{0}}} = −1, β = β^{−}, and

β [\frac{1}{β} ] τ_{x_{i_{0}}} = 1 − \frac{β^{−}}{β^{+}} which is positive iff β^{−} < β^{+}.

For all cases, a term in the first summation is positive iff β = min{β^{+}, β^{−}}. Similarly we can show that a term in the second summation is positive iff β = min{β^{+}, β^{−}} too. Only the positive terms will reduce the lower bound of the left-hand side of (29). So the left-hand side is bounded below by

D_{i,j} ≥ 1 − 0.5(1 − \frac{min(β^{+}, β^{−})}{max(β^{+}, β^{−})}) − 0.5(1 − \frac{min(β^{+}, β^{−})}{max(β^{+}, β^{−})}) = \frac{min(β^{+}, β^{−})}{max(β^{+}, β^{−})}.

Now we turn to the upper bound which can be proved directly from (30).

E_{i,j} ≤ 1 + 0.5β |[\frac{1}{β} ]| + 0.5β |[\frac{1}{β} ]| ≤ 1 + max(β^{+}, β^{−}) |[\frac{1}{β} ]| .

The lower bound (31) is useful for the stability of the numerical scheme and the upper bound (31) will be used in the proof of the following theorem.

Theorem 4.2 At irregular grid points (x_{i}, y_{j}) in Ω,

F (u_{i,j−1}, u_{i−1,j}, u_{i,j}, u_{i+1,j}, u_{i,j+1}) − (βu_{xx} + βu_{yy} + f)(x_{i}, y_{j}) = O(h) . (32)

Proof: From (24), (25), (26), and (31), we have

F (u_{i,j−1}, u_{i−1,j}, u_{i,j}, u_{i+1,j}, u_{i,j+1}) = (u_{i,j})_{t}
= β(δ_{x}u_{i,j} + δ_{y}u_{i,j} + ∑_{(x_{i_{0}} ,y^{∗})∈S_{irr}} C̃_{x}u_{i,j} + ∑_{(x^{∗},y_{j_{0}} )∈S_{irr}} C̃_{y}u_{i,j})

+f_{i,j} + (u_{i,j})_{t}\{ ∑_{(x_{i_{0}} ,y^{∗})∈S_{irr}} β [\frac{1}{β} ] \frac{τ_{x_{i_{0}}} (x_{i_{0}} − x^{∗})^{2} cos^{2} θ}{2h^{2}}

+ ∑_{(x^{∗},y_{j_{0}} )∈S_{irr}} β [\frac{1}{β} ] \frac{τ_{y_{j_{0}}} (y_{j_{0}} − y^{∗})^{2} sin^{2} θ}{2h^{2}} \}

= β(δ_{x}u_{i,j} + δ_{y}u_{i,j}) + f_{i,j}

+β ∑_{(x_{i_{0}} ,y^{∗})∈S_{irr}} \{C_{x}u_{i,j} + [\frac{u_{t}}{β} ] \frac{τ_{x_{i_{0}}} (x_{i_{0}} − x^{∗})^{2} cos^{2} θ}{2h^{2}} \}

+β ∑_{(x^{∗},y_{j_{0}} )∈S_{irr}} \{C_{y}u_{i,j} + [\frac{u_{t}}{β} ] \frac{τ_{y_{j_{0}}} (y_{j_{0}} − y^{∗})^{2} sin^{2} θ}{2h^{2}} \} + O(h)

= (βu_{xx} + βu_{yy} + f)(x_{i}, y_{j}) + O(h),

which proves (32).

At a regular grid point, the local truncation error is O(h^{2}). At an irregular grid point, the local truncation error is O(h) by Theorem 4.2. The local truncation error from the discretization of time is O((Δt)^{2}) = O(h^{2}). All these imply that the numerical solution has global error O(h^{2}).

5 Numerical Examples

We choose some examples in which the actual solutions are known; therefore, numerical error computations can be obtained to confirm the theoretical result of our method. We choose

Ω = (−1, 1) × (−1, 1), t ∈ [1, ∞),
Ω^{+} = {(x, y) ∈ Ω |x^{2} + y^{2} > 1/4},
Ω^{−} = {(x, y) ∈ Ω |x^{2} + y^{2} < 1/4},
Γ = {(x, y) ∈ Ω |x^{2} + y^{2} = 1/4}.

The exact solution for f = 0 is

u(x, y, t) = \frac{1}{t}e^{−(x^{2}+y^{2})/(4βt)}.  (33)

We give the initial condition when t = 1, and the Dirichlet boundary condition when x and y are equal to 1 or −1. We choose λ=1.75 in (28). For different pairs β^{+} and β^{−}, in t from 1.0 to 1.5, we obtain the data shown in Table 1.

The error is computed using the infinity norm. It shows that when h is divided by 2, the error becomes approximately a quarter of its previous value, which confirms that the numerical solution has second order accuracy.

Acknowledgements. The second author wants to thank Professor H. T. Banks for his advice and encouragement during his visit to the Center for Research in Scientific Computation. The second author also wants to thank Professor X.-B. Lin for his hospitality and suggestions. The first author was partially supported by an NSF grant DMS-96-26703, North Carolina State University FR&PD grant, and the ARO under grant number 39676-MA.

h  β^{+} β^{−}  error  ratio
0.100 1000 1 2.227782D-04 –
0.050 1000 1 5.391984D-05 4.13
0.025 1000 1 1.307318D-05 4.12
0.100 1 1000 2.392997D-05 –
0.050 1 1000 6.703319D-06 3.57
0.025 1 1000 1.766744D-06 3.79
0.100 5  1 2.629529D-04 –
0.050 5  1 6.351060D-05 4.14
0.025 5  1 1.550294D-05 4.10
0.100 1  5 5.461059D-05 –
0.050 1  5 1.309861D-05 4.17
0.025 1  5 3.263757D-06 4.01

Table 1: Accuracy of computations

References

[1] N. Gershenfeld, The Nature of Mathematical Modeling, Cambridge Univ. Press, Cambridge, 1999.

[2] R. J. LeVeque and Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31 (1994), 1019 - 1044.

[3] Z. Li and A. Mayo, ADI method for heat equations with discontinuities along an arbitrary interface, in Proc. Symp. Appl. Math. Vol. 48 (W. Gautschi ed.), AMS, 311-315, 1993.

[4] K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge Univ. Press, Cambridge, 1994.

Zhilin Li Center for Research in Scientific Computation and Department of Mathematics North Carolina State University Raleigh, NC 27695-8205, USA

Yun-Qiu Shen Department of Mathematics, Western Washington University Bellingham, WA 98225-9063, USA e-mail: yqshen@cc.wwu.edu