Level Set Function

In this approach, an interface is represented by the zero level set of function φ. The following is the definition of φ on whole domain:

φ(x) < 0 if x ∈ Ω⁻ φ(x) = 0 if x ∈ Γ φ(x) > 0 if x ∈ Ω⁺

We call φ is a signed distance function if φ(x) is the distance from x to the interface. We will use re-initialization process[10] to modified level set function into signed distance function so that use φ to define the outward unit normal vector n by

n = ∇φ

|∇φ|, (8)

and the curvature κ by

κ = ∇ · n = ∇ · µ ∇φ

|∇φ|

¶

. (9)

If X is on the interface but is not a grid point, we can still compute them by the interpolation method using specific points from the four corners of the rectangle that contains X.

3 One-Dimensional Elliptic Interface Prob-lems

In this section, we discuss with one-dimensional elliptic equation. Extending to two-dimensional elliptic equations is using dimension by dimension. The key idea of IIM is to avoid grid generation by correcting finite difference in the neighborhood of the interface. We only show the work about discrete form of (d²u/dx²)_j. The discrete form of (du/dx)_j can be easily obtained by the same way.

Assume the interface is located at x = α(Fig. 2). Let σ = α − xj

h . (10)

Clearly, we have 0 ≤ σ ≤ 1.

Figure 2: Uniform grid with interface located at x = α

The two jump conditions involving u and ux can be obtained from the original elliptic equation. A general jump conditions across the interface can be written as

[γu]_x = γ⁺u⁺− γ⁻u⁻= A (11) and

· βdu

dx

¸

x

= β⁺du dx

+

− β⁻du dx

−

= B. (12)

As mention in Sec. 2.2, the superscripts “+” and “−” represent the variables at the right and left side of the interface, respectively.

There are two kinds of grid points. One is called a regular point if the finite difference formula at this point only involves grid points on the same side of the interface. Otherwise, it is an irregular point. If the gird point i is a regular point, we use standard finite difference directly. For example:

µd²u dx²

¶

j

= u_j−1− 2u_j + u_j+1

h² + O(h²) (13)

or

µd²u dx²

¶

j

= −u_j−2+ 16u_j−1− 30u_j + 16u_j+1− u_j+2

12h² + O(h⁴). (14)

These two difference formulas can be easily obtained by the Taylor expansion.

3.1 Difference Formula for d

²

u/dx

²

at Irregular Points

Finite difference approximations for d²u/dx² at the irregular point j are considered by using a stencil of n points on the left side of interface and m points on the right(Fig. 2).

3.1.1 Difference Formulas with Four-point Stencil (n = m = 2)

Thus, we have to determine d_ks and the correction term on the right hand side of Eq. (15). By using the Taylor expansion at x = α⁻, we have Substituting Eqs. (11) and (12) into Eq. (20) and (21), we have

u_j+1 = A

Substituting Eqs. (18), (19), (22), and (23) into Eq. (15) leads to can conclude that the (d²u/dx²) at irregular points is O(h). To determine d_ks, it is necessary to solve the linear system equation as follows



Figure 3: New Cartesian grid with ex_k−1 = −x_j+2, ex_k = −x_j+1 , ex_k+1 = −x_j, and ex_k+2 = −x_j−1.

Clearly, dks are functions of σ and jump parameters: γ⁺, β⁺, cγ, and cβ. Finally, we’ve determined all d_ks. Eq. (15), together with Eq. (26), is an explicit difference formula for O(h) approximation to (d²u/dx²)_j. Moreover, it shows that the current formula at irregular points does not have singularity, even for the special cases of Γ coinciding with σ = 0 or σ = 1.

We can obtain (du/dx)_j by the same way. The general formulas for first and second derivatives terms are

µdu dx

¶

j

= (d₋₁− 2)u_j−1+ (d₀+ 2)u_j+ d₁u_j+1+ d₂u_j+2+ d_AA + hd_BB

2h + O(h),

µ (27) d²u dx²

¶

j

= d₋₁u_j−1+ d₀u_j + d₁u_j+1+ d₂u_j+2+ d_AA + hd_BB

h² + O(h), (28)

where d_ks are the same in Eq. (26).

3.1.2 Irregular Point Located at Right Side of Interface

The finite difference formulas at grid point j + 1 can be obtained from Eqs.

(27) and (28). Instead of using the method in previous section, we use the coordinate transformation (Fig. 3).

Let

v(x) = u(−x), e

γ(x) = γ(−x),

β(x) = β(−x),e (29)

[eγv]_x = eA, [ eβv⁰]_x = eB,

then the two jump conditions are

Consider the derivative of function v defined on the Cartesian grid ex_k. By the same thought in Sec. 3.1.1, we can easily obtain

µd²v

3.2 Difference Formulas with a General n + m Grid Stencil

In this section, we want to obtain arbitrary high order difference formula scheme. The difference formulas at irregular point j for n + m points are considered. In order to get uniform order, we often take n = m. The method we used is matched polynomial interpolation. We only provide the discrete form of d²u/dx² here.

The polynomial on the left side of Γ, interpolating through n grid points can be written as

P⁻(x) =

−n+1X

k=0

l_k(x)u_i+k+ a_nR(x), (33) where a_n is an undetermined coefficient to be decided, and

R(x) =

−n+1Y

k=0

(x − xi+k). (34)

l_k(x) is the Lagrange polynomial, i.e.

l_k(x) = Similarly, the polynomial on the right side of Γ, interpolating through m grid points can be written as

P⁺(x) = Xm k=1

r_k(x)u_i+k+ b_mQ(x), (36) where b_m is an undetermined coefficient to be decided, and

Q(x) = Ym k=1

(x − x_i+k). (37)

r_k(x) is the Lagrange polynomial, i.e.

r_k(x) = Our thought is that use the relation

µd²u

Thus, we only have to determine the unknown a_n. Once we find a_n, we can obtain the difference formula (d²u/dx²)_j, and then we get (d²u/dx²)_j+1 by the method in Sec. 3.1.2.

Substituting Eqs. (33) and (36) into Eq. (11) leads to

γ⁺

Rearrange the equation above we get

c₁₁a_n+ c₁₂b_m = β₁, (39)

Again, Substituting Eqs. (33) and (36) into Eq. (12) leads to

β⁺

Rearrange the equation above we get

c₂₁a_n+ c₂₂b_m = β₂, (41)

Solving Eqs. (39) and (41), we have

a_n= Xm

µ_ku_i+k+ ξ⁻_AA + ξ⁻_BB, (43)

where µ_k=

½ ₁

J{γ⁻β⁺Q⁰(α)l_k(α) − γ⁺β⁻Q(α)l⁰_k(α)} for k = −n + 1, · · · , 0

γ⁺β⁺

J {−Q⁰(α)r_k(α) + Q(α)r⁰_k(α)} for k = 1, · · · , m, ξ_A⁻= 1

J{β⁺Q⁰(α)}, (44)

ξ_B⁻= −1

J{γ⁺Q(α)},

J = γ⁺β⁻R⁰(α)Q(α) − γ⁻β⁺Q⁰(α)R(α).

In fact, if we take n = m = 2, we will have the same result as present in the previous section.

3.3 Numerical Results

We use four versions of current IIM tested in this thesis. In method C and D, grid points j − 1 and j + 2 are irregular points, but we can treat only j and j + 1 as irregular points and use fourth-order one-sided difference formula for j − 1 and j + 2.

Methods Order at regular Order at irregular Expected global order grid points grid points

Method A O(h²) O(h) O(h²)

Method B O(h²) O(h²) O(h²)

Method C O(h⁴) O(h) O(h²)

Method D O(h⁴) O(h²) O(h³)

Table 1: Four immersed interface method

Example 3.1 In this example, we use the IIM to solve the following prob-lem:

d²u

dx² + κu = βδ(x − α) , x ∈ (−0.5, 0.5), (45) where κ is discontinuous across the interface located at x = α:

κ(x) =

½ (α₁)² if − 0.5 < x ≤ α,

(α2)² if α < x < 0.5. (46)

The boundary condition is

u(−0.5) = u(0.5) = 0. (47)

An alternative way to state Eq. (45) requires that u(x) satisfies the equation d²u

dx² + κu = 0 , x ∈ (−0.5, α) ∪ (α, 0.5), (48) excluding the interface located at x = α, together with boundary conditions (47) and two natural jump conditions:

[u]_x = 0, (49)

The exact solution is

u_ex(x) =

Figure 4: (a)Comparison of the exact solution u_exand the numerical solution u(Method D with N = 80). (b) Numerical error(Method D with N = 640).

Table 2 shows the maximum-norm errors of four methods, the correspond-ing ratios and orders. Note that in order to compare the grid refinement results with the same conditions, all results are compared between grids N and N/4 because they have the same σ. So the ratio and order are defined as following:

Ratio = kE_N/4k_∞ kENk∞

, Order = log₄

µkE_N/4k_∞ kE_Nk_∞

¶ .

N σ Method A Method B

kE_nk_∞ Ratio Order kE_nk_∞ Ratio Order

20 1/3 9.795371 9.884993

40 2/3 9.537e-1 9.369e-1

80 1/3 2.135e-1 45.875 2.76 2.135e-1 46.79 2.77

160 2/3 5.277e-2 18.073 2.09 5.189e-2 18.05 2.08

320 1/3 1.297e-2 16.466 2.02 1.288e-2 16.39 2.01

640 2/3 3.275e-3 16.112 2.00 3.222e-3 16.10 2.00

N σ Method C Method D

kE_nk_∞ Ratio Order kE_nk_∞ Ratio Order

40 2/3 3.673e-1 3.762e-1

80 1/3 1.269e-2 8.797e-3

160 2/3 2.709e-3 135.5 3.54 1.016e-3 369.9 4.27

320 1/3 2.198e-4 57.72 2.92 2.521e-5 348.9 4.22

640 2/3 1.255e-4 21.57 2.21 1.112e-5 91.44 3.26

Table 2: Comparison of numerical errors

4 Two-Dimensional Elliptic Interface Prob-lems

We use a dimension by dimension approach to solve the two-dimensional problems. To compute two-dimensional problems, the grid points are classi-fied into four categories in x-direction:

1. Regular point;

2. Irregular point located on left side of interface;

3. Irregular point located on right side of interface;

4. Irregular point near two interface.

The definition in y-direction is similar to x-direction.

For regular point away from the interface, the derivatives with respect to x and y are approximated by standard central difference:

µd²u Remember that in 1D problems, we need two jump conditions: [γu]x and [βu_x]_x. But we can’t get these two necessary jump conditions from original equation directly, so we still have to make some effort to get these two jump conditions.

4.1 Poisson Equation with Interface

There are two natural jump conditions we can get from original poisson equation:

where s is a parameter of the interface. From Eq. (54), we obtain

·∂u

∂s

¸

= w⁰(s). (56)

Assume that n = (n₁, n₂) is the unit normal vector on Γ. Thus s = (−n₂, n₁) is the unit tangential vector on Γ. Hence, Eqs. (55) and (56) lead to:

[−n₂u_x+ n₁u_y] = w⁰(s), (57) [n₁u_x+ n₂u_y] = v(s). (58) By the definition, we have

−n₂[u_x] + n₁[u_y] = w⁰(s), (59) n₁[u_x] + n₂[u_y] = v(s). (60) Rearrange Eqs. (59) and (60) we have:

[u_x] = n₁v(s) − n₂w⁰(s), (61) [u_y] = n₂v(s) + n₁w⁰(s). (62) Since v(s) and w⁰(s) are known, we only have to calculate two unknowns, n1

and n₂, by level set method.

Note that all jump conditions for partial derivative of u are NOT in x or y-direction. So we have to use Eq. (7) to derive jump conditions for partial derivative in x or y-directions.

4.2 Elliptic Equation with Interface

The natural jump conditions for elliptic equation are

[u] = w(s), (63)

· β∂u

∂n

¸

= v(s). (64)

Again, form Eq. (63) we have

· β∂u

∂s

¸

= w⁰(s). (65)

By the same thought in Sec. 4.1, we have

£(βn²₁ + n²₂)u_x¤

= n₁v(s) − n₂w⁰(s) − [(β − 1)n₁n₂u_y], (66)

£(n²₁+ βn²₂)u_y¤

= n₁w⁰(s) − n₂v(s) − [(β − 1)n₁n₂u_x]. (67) For finite difference approximation of x derivatives at an irregular point, the jump condition (66) is used. So we have to decide the y derivative term on the right hand side of Eq.(66). We evaluate [uy] by one-sided difference at an order of accuracy which is consistent with the order of the overall calculations.

4.3 Numerical Results

4.3.1 Poisson Equation

Example 4.1 We use the example, which was used by Leveque and Li[3]

to test IIM in this thesis.

u_xx+ u_yy = Z

Γ

2δ(x − X(s))(y − Y (s))ds − 1 < x, y < 1, (68) where the interface Γ is a circle defined by x²+y² = 1/4. We can easily obtain unit normal vector n = (2x, 2y) where (x, y) ∈ Γ. The Dirichlet boundary condition is specified by using the exact solution:

uex(x, y) = The jump conditions at all points on Γ are

[u] = 0, (70)

where (x, y) is on the interface Γ. We also test delta function method in this example. Assume f has a delta function singularity along the interface Γ by discrete delta function. For example:

f (x, y) = Z

Γ

C(s)δ(x − X(s))δ(y − Y (s))ds, (74) where (X(s), Y (s)) is the arc-length parameterization of Γ. We use the dis-crete delta function

to calculate fi,j, and the form of which is fi,j =

Xm k=1

C(sk)dh(xi− Xk)dh(yj− Yj)∆sk. (76) Fig. 5 shows that the main error in the computations are originated from the interface. This demonstrates the importance of using higher-order method

Figure 5: Contour of numerical error

Methods DFM Method A Method B

kEnk∞ Order kEnk∞ Order kEnk∞ Order

30 × 30 0.031 2.204e-3 6.867e-4

60 × 60 0.015 0.99 2.873e-4 2.94 1.072e-4 2.68

120 × 120 0.008 0.96 5.312e-5 2.43 3.286e-5 1.70

240 × 240 0.004 0.93 1.225e-5 2.16 8.134e-6 2.01

Table 3: Comparison of numerical errors

Example 4.2 In this example, we consider the discontinuous Poisson prob-lem with the elliptic interface:

Γ : x² a² +y²

b² = 1, and we use the notations

4u_± = f_± in Ω_±,

[u] = w(s) on Γ, (77)

[u_n] = v(s) on Γ, u = u0 on ∂Ω.

We derive the jump conditions [u] and [βun] from the exact solution. Four different examples as shown in Table 4 are tested. Unlike in Example 4.1, the solution in this example is discontinuous.

Remember that in order to calculate the jump conditions (61) and (62), we need the unit normal vector n = (n₁, n₂). In this example, we use re-initialization process mentioned in Sec. 2.

Case 1 Case 2 Case 3 Case 4

u₋ 1 x²− y² e^xcos(y) sin(x) cos(y)

u+ 1 + log(2p

x²+ y²) 0 0 0

f− 0 0 0 −2 sin(x) cos(y)

f₊ 0 0 0 0

Table 4: Four test cases for Eq. (77).

In Case 2, since the exact solution is a polynomial, the maximum errors are machine errors. Note that for Case 4, we use the different grids since the Cartesian grid cannot fetch the interface behavior if we use the same grid for other cases.

N × N kE_nk_∞ Order N × N kE_nk_L² Order

Case 1 30 × 30 5.355e-4 30 × 30 5.004e-4

60 × 60 1.758e-4 1.61 60 × 60 1.153e-4 2.12 120 × 120 3.692e-5 2.25 120 × 120 3.036e-5 1.98 240 × 240 1.239e-5 1.58 240 × 240 6.709e-6 2.18

Case 2 30 × 30 3.608e-16 30 × 30 2.368e-16

60 × 60 4.441e-16 60 × 60 4.400e-16

120 × 120 1.443e-15 120 × 120 6.002e-16 240 × 240 1.221e-15 240 × 240 4.368e-16

Case 3 30 × 30 1.123e-4 30 × 30 9.544e-5

60 × 60 5.723e-5 0.97 60 × 60 4.102e-5 1.22 120 × 120 7.524e-6 2.93 120 × 120 6.330e-6 2.70 240 × 240 2.637e-6 1.51 240 × 240 2.170e-6 1.54

Case 4 40 × 40 1.789e-5 40 × 40 8.451e-6

80 × 80 2.918e-6 2.62 80 × 80 1.191e-6 2.82 160 × 160 6.865e-7 2.09 160 × 160 2.667e-7 2.16 320 × 320 1.028e-7 2.74 320 × 320 2.830e-8 3.24 Table 5: Comparison of numerical errors by method A with result for a = 0.6, b = 0.4.

4.3.2 Elliptic Equation

Example 4.3 We use the example, which was used by Leveque and Li[3]

to test IIM in this thesis. An elliptic equation with a delta function source term and with a discontinuous coefficient β as follows:

∇ · (β∇u) = f (x, y) + C

For the current case, the jump conditions on Γ are

[u] = 0, (79)

Here we derive the necessary two jump conditions for partial derivative by the method in Sec. 4.2.

C = 0.1 40 × 40 80 × 80 160 × 160 320 × 320

b kE_nk_∞ kE_nk_∞ Order kE_nk_∞ Order kE_nk_∞ Order 10.0 8.15e-5 2.53e-5 1.69 5.76e-6 2.14 1.38e-6 2.06

5.0 1.58e-4 4.97e-5 1.67 1.12e-5 2.15 2.71e-6 2.06

1.0 6.99e-4 2.31e-4 1.59 5.04e-5 2.20 1.24e-5 2.02

0.01 0.05985 0.02128 1.49 0.00431 2.30 0.00112 1.94 0.005 0.11961 0.04255 1.49 0.00861 2.30 0.00225 1.94

Table 6: Comparison of numerical errors

According to the exact solution, as b decreases, the maximum magnitude of |u(x, y)| increases. Therefore, the computational errors will increase when the value of b decreases.

Example 4.4 In this example, we consider the discontinuous elliptic equa-tion with elliptic interface:

Γ : x² a² +y²

b² = 1

Case 1 Case 2 Case 3

u₋ x²+ y² x²+ y² x²− y²

u₊ ¹⁻⁸⁰¹₄⁻¹⁰¹ + ^(x2+y2)22 ₁₀^+x2+y2 +^{0.1 log(2}

√x²+y²)

10 0 sin(x) cos(y)

β₋ x²+ y²+ 1 e^x− y 1

β₊ 10 0.5 2

f₋ 8(x²+ y²) + 4 2e^x(x + 2) − 6y 0

f₊ 8(x²+ y²) + 4 0 −2 sin(x) cos(y)

Table 7: Three test cases.

N × N kE_nk_∞ Order N × N kE_nk_L² Order

Case 1 30 × 30 2.169e-4 30 × 30 2.552e-4

60 × 60 5.600e-5 1.95 60 × 60 6.527e-5 1.97 120 × 120 1.508e-5 1.89 120 × 120 1.736e-5 1.91 240 × 240 3.135e-6 2.27 240 × 240 3.706e-6 2.23

Case 2 30 × 30 2.776e-16 30 × 30 1.086e-16

60 × 60 6.661e-16 60 × 60 1.833e-16

120 × 120 6.661e-16 120 × 120 1.535e-16 240 × 240 1.332e-15 240 × 240 4.247e-16

Case 3 30 × 30 2.093e-5 30 × 30 2.214e-5

60 × 60 8.626e-6 1.23 60 × 60 8.084e-6 1.45 120 × 120 1.449e-6 2.57 120 × 120 1.572e-6 2.36 240 × 240 5.628e-7 1.36 240 × 240 6.195e-7 1.34 Table 8: Comparison of numerical errors by method A with result for a = 0.8, b = 0.2.

In Case 2, because the exact solution is a polynomial, the maximum errors

4.3.3 Application of Heat Equation

Example 4.5 We use the example used by Shen and Li[6], the heat equa-tion with an interface as follows

ut = ∇ · (β∇u), Ω = [−1, 1] × [−1, 1], t ∈ [0, ∞] (81) where

β(x, y) =

( β⁻ if (x, y) ∈ Ω⁻,

β⁺ if (x, y) ∈ Ω⁺. (82) We give the initial condition

u(x, y, 1) = exp µ

−x²+ y² 4β

¶

and the Dirichlet boundary condition when (x, y) ∈ ∂Ω and two natural jump conditions [u] and [βu_n] by the exact solution:

u(x, y, t) = 1 t exp

µ

−x²+ y² 4βt

¶

. (83)

We use the Crank-Nicolson method in this example. The steps of our algo-rithm can be outlined as follows:

Step 1. Reinitialize φ to be an exact signed distance function by solving the equation, φ_t = sgn(φ₀)(1 − |∇φ|) to steady state.

Step 2. Compute outward normal vector n by Eq. (8). We compute this value at grid points neighboring the interface, then we interpolate its value on the interface whenever it is needed. Eq. (8) is numerically solved using center difference approximations to the partial derivatives of φ.

Step 3. Use φ to determine a flag matrix F . For example, let F (i, j) = 0 if (x_i, y_j) is regular point;

F (i, j) = 1 if (x_i, y_j) is irregular point.

Then we can use F to determine what kind the grid points are. Also, we compute the distance form irregular (xi, yj) to the interface in either the vertical or horizontal direction.

Step 4. Use the Crank-Nicolson method as following:

uⁿ⁺¹_i,j − uⁿ_i,j

4t = 1

2(β_i,j4uⁿ⁺¹_i,j + β_i,j4uⁿ_i,j).

Rearrange the equation above, we have

4t(β_i,j4uⁿ⁺¹_i,j ) − 2uⁿ⁺¹_i,j = −4t(β_i,j4uⁿ_i,j) − 2uⁿ_i,j (84) Use F to distinguish grid points. For regular point, we use five-point Laplacian to discretize 4uⁿ⁺¹_i,j on the left hand side of Eq. (84). For irregular point, we discretize 4uⁿ⁺¹_i,j term by IIM in this thesis. This scheme is always numerically stable and convergent but usually more numerically intensive as it requires solving a linear system AU = b of numerical equations on each time step.

Step 5. Since our interface is fixed, the matrix A is fixed. We repeat Step 4 to get the solution u at next time step.

Case 1.

N × N β⁺ β⁻ kEnk∞ Order kEnk_L² Order

30 × 30 1000 1 3.928e-5 4.413e-5

60 × 60 1000 1 1.029e-5 1.93 1.166e-5 1.92

120 × 120 1000 1 2.669e-6 1.95 3.038e-6 1.94 240 × 240 1000 1 6.996e-7 1.93 7.974e-7 1.93

30 × 30 1 1000 3.532e-5 4.253e-5

60 × 60 1 1000 9.832e-6 1.85 1.189e-5 1.84

120 × 120 1 1000 2.522e-6 1.96 3.015e-6 1.98 240 × 240 1 1000 6.304e-7 2.00 7.540e-7 2.00

30 × 30 5 1 4.500e-5 5.076e-5

60 × 60 5 1 1.192e-5 1.92 1.352e-5 1.91

120 × 120 5 1 3.084e-6 1.95 3.510e-6 1.95

240 × 240 5 1 8.012e-7 1.94 9.131e-7 1.94

Table 9: Numerical results with Γ : x²+ y² = 1/4, ∆t = h, T = 2.

Case 2.

N × N β⁺ β⁻ kE_nk_∞ Order kE_nk_L² Order

30 × 30 1000 1 5.220e-5 4.855e-5

60 × 60 1000 1 1.352e-5 1.95 1.287e-5 1.92

120 × 120 1000 1 3.478e-6 1.96 3.336e-6 1.95 240 × 240 1000 1 8.674e-7 2.00 8.280e-7 2.01

30 × 30 1 1000 5.023e-5 5.092e-5

60 × 60 1 1000 1.363e-5 1.88 1.363e-5 1.90

120 × 120 1 1000 3.492e-6 1.88 3.439e-6 1.99 240 × 240 1 1000 8.963e-7 1.96 8.834e-7 1.96

30 × 30 5 1 5.884e-5 5.544e-5

60 × 60 5 1 1.537e-5 1.94 1.475e-5 1.91

120 × 120 5 1 3.947e-6 1.96 3.806e-6 1.95

240 × 240 5 1 9.903e-7 1.99 9.528e-7 2.00

Table 10: Numerical results with Γ : ¡ _x

0.8

¢₂ +¡ _y

0.2

¢₂

= 1, ∆t = h, T = 2.

5 Conclusion

The advantage of the IIM in this thesis is that the finite difference formulas at irregular points are expressed in an explicit form, so they can be applied to difference problems without modifications. But there are still hard work when we use the natural jump condition [βu_n], so our method will fail if the problem with concave interface.

Our future work is to calculate [u_n] by the relation [u_n] = [βu_n] − [β]u⁻_n

β⁺ .

Thus, we have to solve the unknown term u⁻_n. Li[4] offered the method by a similar method introduced in original IIM.

There are two challenges. First one is modify our IIM by the method introduced above. Second one is use the modified IIM to solve the Stefan problems[11, 12, 13].

Acknowledgments The method and most notations we used in this thesis are referred by Zhong[1]. Especially we thank Tsai for providing the code of re-initialization process.

References

[1] X. Zhong. (2007). A new high-order immersed interface method for solv-ing elliptic equations with imbedded interface of discontinuity, Journal of Computational Physis. 225, 1066-1099.

[2] Z. Li, K. Ito. (2006). The Immersed Interface Method : Numerical So-lutions of PDEs Involving Interfaces and Irregular Domains, SIAM.

[3] R.J. LeVeque. (1994). Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J.

Numer. Anal. 31, 1001-1025.

[4] Z. Li. (1998). A Fast Iterative Algorithm for Elliptic Interface Problems, Journal on Numerical Analysis. 35, 230-254.

[5] Z. Li. (1997). Immersed interface methods for moving interface problems, Numerical Algorithms. 14, 269-293.

[6] Y. Shen, Z. Li, A numerical method for solving heat equations involving interfaces, Electron. J. Diff. Eqns., Conf. 03 (1999) 100-108.

[7] Y.C. Zhou, S. Zhao, M. Feig, G.W. Wei. (2006). High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources, J. Comput. Phys. 213, 1-30.

[8] A. Wiegmann, Kenneth P. Bube. (2000). The Explicit-Jump Immersed Interface Method: Finite Difference Methods for PDES with Piecewise Smooth Solutions, SIAM , Journal on Numerical Analysis, 37, 827-862.

[9] A. Mayo. (1984). The fast solution of Poissons and biharmonic equations on irregular regions, SIAM J. Numer. Anal. 21, 285-299.

[10] S. Osher, R. Fedkiw. (2003). Level Set Methods and Dynamic Implicit Surfaces, Springer.

[11] S. Chen, B. Merriman, S. Osher and P. Smereka. (1997). A simple level set method for solving Stefan problems, Journal of Computational Ph-ysis. 135, 8-29.

[12] Z. Li. (1997). Immersed interface methods for moving interface problems, Numerical Algorithms. 14, 269-293.

[13] Z. Li, B. Soni. (1999). Fast and accurate numerical approaches for stefan problems and crystal growth, Numerical Heat Transfer, Part B, 35,

在文檔中 橢圓介面方程的數值研究 (頁 10-0)