Tailored Finite Point Method for a Singular Perturbation Problem with ...
J Sci Comput DOI 10.1007/s10915-009-9292-2
Tailored Finite Point Method for a Singular Perturbation Problem with Variable Coefficients in Two Dimensions Houde Han · Zhongyi Huang
Received: 1 December 2008 / Revised: 11 March 2009 / Accepted: 24 March 2009 © Springer Science+Business Media, LLC 2009
Abstract In this paper, we propose a tailored-finite-point method for a type of linear singular perturbation problem in two dimensions. Our finite point method has been tailored to some particular properties of the problem. Therefore, our new method can achieve very high accuracy with very coarse mesh even for very small ε, i.e. the boundary layers and interior layers do not need to be resolved numerically. In our numerical implementation, we study the classification of all the singular points for the corresponding degenerate first order linear dynamic system. We also study some cases with nonlinear coefficients. Our tailored finite point method is very efficient in both linear and nonlinear coefficients cases. Keywords Tailored finite point method · Singular perturbation problem · Singular point
1 The Model Problem We are interested in the numerical simulation of singular perturbation problems with variable coefficients. In this paper, we shall focus on problems in two dimensions. Let us consider the following problem: ∂uε ∂uε (x) + q(x) (x) = 0, ∂x1 ∂x2 ∂uε = g, ∂n
−εuε (x) + p(x) uε |D = f,
x ∈ ,
N
H. Han was supported by the NSFC Project No. 10471073. Z. Huang was supported by the NSFC Project No. 10676017, the National Basic Research Program of China under the grant 2005CB321701. H. Han · Z. Huang () Dept. of Mathematical Sciences, Tsinghua University, Beijing 100084, China e-mail: [email protected] H. Han e-mail: [email protected]
(1.1) (1.2)
J Sci Comput
where ⊂ R2 is a bounded domain. ε is a small positive parameter, p and q are two piecewise smooth functions given in domain . When 1 the solution may have some boundary layers and some interior layers in downwind directions. These layers are characterized by rapid transitions in the solution, and are thus difficult to capture in a numerical approximation without using a large number of unknowns. Also, such layers tend to cause spurious oscillations in a numerical solution to the problem. Methods for the numerical solution of problems such as (1.1)–(1.2) in bounded or unbounded domains that attempt to deal with these difficulties have been developed by many mathematicians, see e.g., [1, 2, 4, 5, 9, 10, 12, 13, 16–18, 20–23]. If ε 1, in the ‘outer domain’, i.e. in the domain out of the boundary layers and the interior layers, the second order term ‘−εuε (x)’ can be dropped from the equation. One obtains the first order partial differential equation p(x)
∂uε ∂uε (x) + q(x) (x) = 0. ∂x1 ∂x2
The characteristic curves of this equation are given by the solution of the following ordinary differential system dx2 dx1 = p(x), = q(x). (1.3) dt dt It is well-known that there is a complete classification of the singular points of the ordinary differential system (1.3) when p and q are linear functions, p(x) = ax1 + bx2 and q(x) = cx1 + dx2 [11], i.e. dx2 dx1 = ax1 + bx2 , = cx1 + dx2 , dt dt with a, b, c, d ∈ R. The equation a − λ b D(λ) = = λ2 − (a + d)λ + ad − bc = 0 c d − λ
(1.4)
is called the characteristic equation of (1.4). Let α = −(a + d), β = ad − bc, then D(λ) = λ2 + αλ + β = 0, it has two roots −α ± α 2 − 4β λ1,2 = . 2 We have the following classification of the singular points: • • • • • • • • •
Saddle Point: β < 0; Stable Node: β > 0, α > 0, α 2 − 4β > 0; Stable Focus: β > 0, α > 0, α 2 − 4β < 0; Unstable Node: β > 0, α < 0, α 2 − 4β > 0; Unstable Focus: β > 0, α < 0, α 2 − 4β < 0; Critical Stable Node/Degenerate Stable Node: β > 0, α > 0, α 2 − 4β = 0; Critical Unstable Node/Degenerate unstable Node: β > 0, α < 0, α 2 − 4β = 0; Center: β > 0, α = 0; Infinite Singularities: β = 0.
In this paper, we propose an effective numerical method, called the “tailored finite point method”, to solve the problem (1.1)–(1.2) numerically. Especially, we shall use our method
J Sci Comput
to simulate the problem with different ‘singular points’. We will also show some numerical results for the case in which the coefficients p(x) and q(x) are nonlinear functions.
2 Finite Point Method In this section we present our new approach to construct a discrete scheme for the (1.1). The new scheme is called a “a tailored finite point method”, because the finite point method has been tailored to some particular properties or solutions of the problem. The finite point method [3, 14, 15, 19] is a development of finite difference method, in which the meshless technique is emphasized. For one dimensional singular perturbation problem, the method is very close to the method of “exponential fitting” discussed in [1, 10, 16] and [20] (Chap. I, Sect. 2). We have presented the tailored finite-point schemes for the Laplace equation and a singular perturbation problem with constant coefficients in [7], we also presented the TFPM for Helmholtz equation in heterogeneous media [6] and interface problems [8]. Here we shall construct a tailored-five-point scheme for the singular perturbation equation (1.1). First, we approximate the coefficients p(x) and q(x) by piecewise constant functions, i.e. we assume that p(x) ≈ p(x 0 ) ≡ p0 and q(x) ≈ q(x 0 ) ≡ q0 in the domain 0 (cf. Fig. 1). Let p0 x1 q0 x2 uε (x1 , x2 ) = v ε (x1 , x2 ) exp + . (2.1) 2ε 2ε Then v ε satisfies: −v ε +
p02 + q02 ε v = 0. 4ε 2
(2.2)
Let α0 =
p0 , 2ε
β0 =
q0 , 2ε
μ0 =
α02 + β02 .
(2.3)
The solution v ε of (2.2) can be expanded at x = x 0 as: v ε (x1 , x2 ) = a0 I0 (μ0 r) +
∞ n=1
Fig. 1 The location mesh of points x 0 , x 1 , . . . , x 9
In (μ0 r) (an cos nθ + bn sin nθ ) ,
(2.4)
J Sci Comput
Fig. 2 The numerical solutions of TFPM for Example 3.2 (saddle point): h = 1/32
where (r, θ ) is the ordered pair of polar coordinates of x = (x1 , x2 ) with the pole at x 0 and In is the n-th order modified Bessel function of first kind. Therefore the solution uε of (1.1) can be expanded at x = x 0 as: uε (x1 , x2 ) = b00 + er(α0 cos θ+β0 sin θ) ∞ 0
0 0 × a0 I0 (μ0 r) + In (μ0 r) an cos nθ + bn sin nθ .
(2.5)
n=1
We take the first few terms in (2.5) to construct an approximation u |0 ≈ ε
b00
+e
r(α0 cos θ+β0 sin θ)
a00 I0 (μ0 r) +
N n=1
In (μ0 r)
an0 cos nθ
+ bn0 sin nθ
.
(2.6)
J Sci Comput
Fig. 3 The numerical solutions of TFPM for Example 3.3 (stable node): h = 1/32
Now, we define the tailored 5-point scheme as follows. Set
W 4 = u(r, θ ) u = b00 + er(α0 cos θ+β0 sin θ) a00 I0 (μ0 r) + I1 (μ0 r) a10 cos θ + b10 sin θ , ∀an0 , bn0 ∈ R, n = 0, 1 . Consider μ0 > 0. First, suppose √ |p0 | ≥ ( 2 + 1)|q0 |,
or
√ |q0 | ≥ ( 2 + 1)|p0 |.
(2.7)
In this case, the tailored finite point equation uses the points x 0 , x 1 , x 2 , x 3 , x 4 . The equation at x 0 is 4 n=0
cn u(x n ) = 0
(2.8)
J Sci Comput
where the numbers cn (n = 0, 1, 2, 3, 4), satisfy 4
cn w(x n ) = 0,
∀w ∈ W 4 .
(2.9)
n=0
Let αh = α0 h, βh = β0 h, and
C1 (μ0 h) = 2I0 (μ0 h) cosh(αh ) − cosh(βh ) .
(2.10)
After some algebra we arrive at the formulas c1 = −
c0 e−αh I0 (μ0 h) − cosh(βh ) , C1 (μ0 h)
c0 e−βh cosh(αh ) − I0 (μ0 h) , C1 (μ0 h) c0 eαh I0 (μ0 h) − cosh(βh ) , c3 = − C1 (μ0 h) c2 = −
c4 = −
c0 eβh cosh(αh ) − I0 (μ0 h) . C1 (μ0 h)
(2.11) (2.12) (2.13) (2.14)
If (2.7) does not hold, the tailored finite point equation uses the points x 0 , x 5 , x 6 , x 7 , x 8 . The equation at x 0 in this case is c0 u(x 0 ) +
8
cn u(x n ) = 0
(2.15)
n=5
where the numbers cn (n = 0, 5, 6, 7, 8), satisfy c0 w(x 0 ) +
8
cn w(x n ) = 0,
∀w ∈ W 4 .
(2.16)
n=5
Let
√
C2 (μ0 h) = 2I0 ( 2μ0 h) cosh(αh + βh ) − cosh(αh − βh ) .
(2.17)
After some algebra we then arrive at the formulas e−αh −βh √ I0 ( 2μ0 h) − cosh(αh − βh ) , C2 (μ0 h) √ eαh −βh cosh(αh + βh ) − I0 ( 2μ0 h) , c6 = −c0 C2 (μ0 h) eαh +βh √ I0 ( 2μ0 h) − cosh(αh − βh ) , c7 = −c0 C2 (μ0 h) √ e−αh +βh cosh(αh + βh ) − I0 ( 2μ0 h) . c8 = −c0 C2 (μ0 h) c5 = −c0
(2.18) (2.19) (2.20) (2.21)
Equations (2.8) and (2.15) give the 5 point tailored finite point scheme at the point x 0 for the problem (1.1).
J Sci Comput
Fig. 4 The numerical solutions of TFPM for Example 3.4 (stable focus): h = 1/32
We can also obtain a tailored 9-point scheme in the similar way. Setting 3 0
0 r(α0 cos θ+β0 sin θ) 0 0 a0 I0 (μ0 r) + W = u(r, θ ) u = b0 + e In (μ0 r) an cos nθ + bn sin nθ ,
8
n=1
∀an0 , bn0 ∈ R, n = 0, 1, 2, 3 , we seek numbers cn (n = 0, . . . , 8), such that 8 n=0
cn w(x n ) = 0,
∀w ∈ W 8 .
(2.22)
J Sci Comput
Fig. 5 The numerical solutions of TFPM for Example 3.5 (center): h = 1/32
Let √
C3 (μ0 h) = 2I0 ( 2μ0 h) cosh(αh ) + cosh(βh ) − 2I0 (μ0 h) cosh(αh − βh ) + cosh(αh + βh ) . (2.23) After some algebra we arrive at the tailored 9-point scheme √ c0 e−αh 2I0 ( 2μ0 h) − cosh(αh + βh ) − cosh(αh − βh ) , 2C3 (μ0 h) √ c0 e−βh 2I0 ( 2μ0 h) − cosh(αh + βh ) − cosh(αh − βh ) , c2 = − 2C3 (μ0 h) √ c0 eαh 2I0 ( 2μ0 h) − cosh(αh + βh ) − cosh(αh − βh ) , c3 = − 2C3 (μ0 h) √ c0 eβh 2I0 ( 2μ0 h) − cosh(αh + βh ) − cosh(αh − βh ) , c4 = − 2C3 (μ0 h)
c1 = −
J Sci Comput
Fig. 6 The numerical solutions of TFPM for Example 3.6 (unstable node): h = 1/32
c5 = −
c0 e−αh −βh −2I0 (μ0 h) + cosh(αh ) + cosh(βh ) , 2C3 (μ0 h)
c6 = −
c0 eαh −βh −2I0 (μ0 h) + cosh(αh ) + cosh(βh ) , 2C3 (μ0 h)
c7 = −
c0 eαh +βh −2I0 (μ0 h) + cosh(αh ) + cosh(βh ) , 2C3 (μ0 h)
c8 = −
c0 e−αh +βh −2I0 (μ0 h) + cosh(αh ) + cosh(βh ) . 2C3 (μ0 h)
For μ0 = 0, i.e. p0 = q0 = 0, (1.1) is reduced to the Laplace equation. Then we should use the standard five-point finite difference scheme or nine-point finite difference scheme for Laplace equation, i.e. u1 + u2 + u3 + u4 − 4u0 = 0, h2
J Sci Comput
Fig. 7 The numerical solutions of TFPM for Example 3.7 (unstable focus): h = 1/32
or 4(u1 + u2 + u3 + u4 ) + (u5 + u6 + u7 + u8 ) − 20u0 = 0. h2
3 Numerical Implementation In this section, we give some numerical examples to show the efficiency of our new method. In the following tables and figures, the ‘exact’ solution u is obtained numerically on a very fine mesh. uFh P 5 is the solution solved by our five point scheme. Furthermore, we let
EhF P ,2 =
uFh P 5 − u L2 () , u L2 ()
EhF P ,ε =
uFh P 5 − u ε, , u ε,
F P 5/9 = the maximum nodal errors of 5-point TFPM or 9-point TFPM, E∞
J Sci Comput
Fig. 8 The numerical solutions of TFPM for Example 3.8 (critical stable node): h = 1/32
where u 2ε, = ε|u|21, + u 2L2 () , and
|u|21,
≡
|∇u(x)|2 dx1 dx2 ,
u 2L2 ()
≡
|u(x)|2 dx1 dx2 .
Example 3.1 (Accuracy test) First, we consider the following problem: −εu + pux1 + qux2 = 0, u|x1 =−1 = x2 + 1, u|x2 =1 = 1 − x1 ,
∀x ∈ ,
u|x2 =−1 = x1 + 1,
(3.1) u|x1 =1 = 1 − x2 ,
(3.2)
J Sci Comput Table 1 Example 3.1: Relative errors of the numerical solutions for TFPM with = 0.01 Mesh size h
1/8
1/16
1/32
1/64
EhF P ,2
3.561E−1
9.299E−2
2.389E−2
5.973E−3
1.9
2.0
2.0
1.214E−1
3.096E−2
7.815E−3
2.0
2.0
2.0
Convergence order EhF P ,ε
4.791E−1
Convergence order
Table 2 Example 3.1: The nodal errors of the numerical solutions by TFPM, h = 1/32
0.5
0.1
0.01
0.001
FP5 E∞
5.741E−2
5.742E−2
5.846E−2
6.024E−2
FP9 E∞
3.473E−2
3.467E−2
3.410E−2
3.439E−2
where
= x = (x1 , x2 ) ∈ R2 −1 < x1 < 1, −1 < x2 < 1 ,
(3.3)
and p(x) = x1 + x2 ,
q(x) = x1 − x2 .
The relative errors of our approximations are shown in Table 1. For our five-point TFPM, it seems that we have the second order convergence rate for L2 norm and energy norm. These convergent rates coincide with the corresponding results in 1D which were proved in [8]. From Table 2, we can find that our TFPM has a uniformly nodal errors as → 0. In what follows, we show numerical results for all kinds of singular points of the dynamic system (1.4). We consider different choices of the coefficients p and q, including cases when p and q are nonlinear functions of x. Example 3.2 (Saddle Point) We consider the problem (3.1)–(3.3) with p(x) = x2 ,
q(x) = x1 .
In this case, we have two real eigenvalues with opposite sign for the linear dynamic system (1.4). The results are shown in Fig. 2. We can see that, although the solution has very sharp boundary layers as → 0, our TFPM can capture them in very coarse mesh without any pseudo-oscillation. Example 3.3 (Stable Node) We consider the problem (3.1)–(3.3) with p(x) = −x1 + 0.5x2 ,
q(x) = 0.5x1 − x2 .
In this case, we have two negative real eigenvalues for the linear dynamic system (1.4). The results are shown in Fig. 3. In this case, the solution has four interior layers, and the origin is a singular point. Our TFPM also can detect these results in coarse mesh.
J Sci Comput
Fig. 9 The numerical solutions of TFPM for Example 3.9 (degenerate stable node): h = 1/32
Example 3.4 (Stable Focus) We consider the problem (3.1)–(3.3) with p(x) = −x1 − x2 ,
q(x) = x1 − x2 .
In this case, we have two conjugated complex eigenvalues with negative real parts for the linear dynamic system (1.4). The results are shown in Fig. 4. Here the solution has four curved interior layers. Our TFPM can get them very clearly. Example 3.5 (Center) We consider the problem (3.1)–(3.3) with p(x) = x2 ,
q(x) = −x1 .
In this case, we have two conjugated imaginary eigenvalues for the linear dynamic system (1.4). The results are shown in Fig. 5.
J Sci Comput
Fig. 10 The numerical solutions of TFPM for Example 3.10 (critical unstable node): h = 1/32
Example 3.6 (Unstable Node) We consider the problem (3.1)–(3.3) with p(x) = x1 + 0.5x2 ,
q(x) = 0.5x1 + x2 .
In this case, we have two positive real eigenvalues for the linear dynamic system (1.4). The results are shown in Fig. 6. Example 3.7 (Unstable Focus) We consider the problem (3.1)–(3.3) with p(x) = x1 − x2 ,
q(x) = x1 + x2 .
In this case, we have two conjugated complex eigenvalues with positive real parts for the linear dynamic system (1.4). The results are shown in Fig. 7. From Figs. 5–7, we can see that our TFPM can capture the sharp boundary layers in very coarse mesh without any artificial oscillation.
J Sci Comput
Fig. 11 The numerical solutions of TFPM for Example 3.11 (degenerate unstable node): h = 1/32
Example 3.8 (Critical Stable Node) We consider the problem (3.1)–(3.3) with p(x) = −x1 ,
q(x) = −x2 .
In this case, we have two double negative real eigenvalues with single elementary divisor for the linear dynamic system (1.4). The results are shown in Fig. 8. In this case and Example 3.9 (cf. Fig. 9), the solutions have four interior layers when → 0. Our TFPM can obtain them very robustly. Especially, the contour curves are very smooth even when h . Example 3.9 (Degenerate Stable Node) We consider the problem (3.1)–(3.3) with p(x) = −x1 ,
q(x) = x1 − x2 .
In this case, we have two double negative real eigenvalues with double elementary divisor for the linear dynamic system (1.4). The results are shown in Fig. 9.
J Sci Comput
Fig. 12 The numerical solutions of TFPM for Example 3.12: h = 1/32
Example 3.10 (Critical Unstable Node) We consider the problem (3.1)–(3.3) with p(x) = x1 ,
q(x) = x2 .
In this case, we have two double positive real eigenvalues with single elementary divisor for the linear dynamic system (1.4). The results are shown in Fig. 10. Example 3.11 (Degenerate Unstable Node) We consider the problem (3.1)–(3.3) with p(x) = −x1 ,
q(x) = x1 − x2 .
In this case, we have two double positive real eigenvalues with double elementary divisor for the linear dynamic system (1.4). The results are shown in Fig. 11. Example 3.12 We consider the problem (3.1)–(3.3) with p(x) = x1 + x2 ,
q(x) = 0.
J Sci Comput
Fig. 13 The numerical solutions of TFPM for Example 3.13: h = 1/32
In this case, we have a line x1 + x2 = 0 of singular points for the linear dynamic system (1.4). The results are shown in Fig. 12. Example 3.13 We consider the problem (3.1)–(3.3) with p(x) = x1 − x2 ,
q(x) = −x1 + x2 .
In this case, we have a line x1 − x2 = 0 of singular points for the linear dynamic system (1.4). The results are shown in Fig. 13. In Examples 3.10–3.13, the solutions have very sharp boundary layers. But our TFPM can get them without any oscillation in very coarse mesh. Example 3.14 We consider the problem (3.1)–(3.3) with nonlinear coefficients p(x) = x1 x2 ,
q(x) = x22 − x14 .
J Sci Comput
Fig. 14 The numerical solutions of TFPM for Example 3.14: h = 1/32
In this case, the origin O is a high order singular point for the nonlinear dynamic system (1.3). The results are shown in Fig. 14. Here the solution has some boundary layers and interior layers. Our TFPM can get all of them with high accuracy in coarse mesh. Example 3.15 We consider the problem (3.1)–(3.3) with nonlinear coefficients p(x) = x1 x2 − 3x13 ,
q(x) = x22 − 6x12 x2 + x14 .
In this case, the origin O is also a high order singular point. The results are shown in Fig. 15.
Example 3.16 We consider the problem (3.1)–(3.3) with nonlinear coefficients p(x) = x13 x2 − x1 x23 ,
q(x) = x12 x22 − x24 + x16 .
J Sci Comput
Fig. 15 The numerical solutions of TFPM for Example 3.15: h = 1/32
In this case, the origin O is also a high order singular point. The results are shown in Fig. 16. From Figs. 15–16, we can see that our TFPM has high accuracy in coarse mesh even the solution has very complex structures. 4 Conclusion In this paper, we present a tailored-finite-point method for a number of linear singular perturbation problems with variable coefficients in two dimensional case. First, we approximate the coefficients with piecewise constants. Then we use the eigenfunction expansion of the local problem with constant coefficients to get the tailored scheme. We consider all kinds of the singular points of the reduced first order linear differential equation, including the saddle point, the stable/unstable node, the stable/unstable focus, etc. Furthermore, we also consider the nonlinear cases. All of our numerical results show that our new method can achieve very high accuracy with very coarse mesh whenever ε is small or large, i.e. we can achieve high accuracy even when our mesh size is much larger than the thickness of the boundary layers and the interior layers. By contrast, as we know, the traditional finite element method
J Sci Comput
Fig. 16 The numerical solutions of TFPM for Example 3.16: h = 1/32
or finite difference method can not get satisfactory numerical results with the same mesh, especially for small ε. We will study the error estimates of our TFPM soon. Acknowledgement
The authors thank Prof. Bruce Kellogg for fruitful discussions on this work.
References 1. Berger, A.E., Han, H.D., Kellogg, R.B.: A priori estimates and analysis of a numerical method for a turning point problem. Math. Comput. 42, 465–492 (1984) 2. Brayanov, I., Dimitrova, I.: Uniformly convergent high-order schemes for a 2D elliptic reaction-diffusion problem with anisotropic coefficients. Lect. Notes Comput. Sci. 2542, 395–402 (2003) 3. Cheng, M., Liu, G.R.: A novel finite point method for flow simulation. Int. J. Numer. Methods Fluids 39, 1161–1178 (2002) 4. Ge, L., Zhang, J.: High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids. J. Comput. Phys. 171, 560–578 (2001) 5. Han, H.D.: The artificial boundary method—numerical solutions of partial differential equations on unbounded domains. In: Li, T., Zhang, P. (eds.) Frontiers and Prospects of Contemporary Applied Mathematics. Higher Education Press/World Scientific, Singapore (2005)
J Sci Comput 6. Han, H., Huang, Z.: A tailored finite point method for the Helmholtz equation with high wave numbers in heterogeneous medium. J. Comput. Math. 26, 728–739 (2008) 7. Han, H., Huang, Z., Kellogg, B.: A Tailored finite point method for a singular perturbation problem on an unbounded domain. J. Sci. Comput. 36, 243–261 (2008) 8. Huang, Z.: Tailored finite point method for the interface problem. Netw. Heterogeneous Media 4, 91–106 (2009) 9. Hemker, P.W.: A singularly perturbed model problem for numerical computation. J. Comput. Appl. Math. 76, 277–285 (1996) 10. Il’in, A.M.: Differencing scheme for a differential equation with a small parameter affecting the highest derivative. Math. Notes 6, 596–602 (1969) 11. Iooss, G., Joseph, D.D.: Elementary Stability and Bifurcation Theory. Springer, New York (1980) 12. Li, J., Chen, Y.: Uniform convergence analysis for singularly perturbed elliptic problems with parabolic layers. Numer. Math. Theor. Methods Appl. 1, 138–149 (2008) 13. Li, J., Navon, I.M.: Uniformly convergent finite element methods for singularly perturbed elliptic boundary value problems: convection-diffusion. Comput. Methods Appl. Mech. Eng. 162, 49–78 (1998) 14. Lin, H., Atluri, S.N.: The Meshless Local Petrov-Galerkin (MLPG) method for solving incompressible Navier-Stokes equations. CMES 2, 117–142 (2001) 15. Mendeza, B., Velazquez, A.: Finite point solver for the simulation of 2-D laminar incompressible unsteady flows. Comput. Methods Appl. Mech. Eng. 193, 825–848 (2004) 16. Miller, J.J.H.: On the convergence, uniformly in ε, of difference schemes for a two-point boundary singular perturbation problem. In: Hernker, P.W., Miller, J.J.H. (eds.) Numerical Analysis of Singular Perturbation Problems, pp. 467–474. Academic Press, San Diego (1979) 17. Morton, K.W.: Numerical Solution of Converction-Diffusion Problems. Applied Mathematics and Mathematical Computation, vol. 12. Chapman and Hall, London (1996) 18. Morton, K.W., Stynes, M., Süli, E.: Analysis of a cell-vertex finite volume method for convectiondiffusion problems. Math. Comput. 66, 1389–1406 (1997) 19. Oñate, E., Idelsohn, S., Zienkiewicz, O.C., Taylor, R.L.: A finite point method in computational mechanics. Applications to convective transport and fluid flow. Int. J. Numer. Methods Eng. 39, 3839–3866 (1996) 20. Roos, H.-G., Stynes, M., Tobiska, L.: Numerical Methods for Singularly Perturbed Differential Equations. Springer, New York (1996) 21. Shishkin, G.I.: A finite difference scheme on a priori adapted meshes for a singularly perturbed parabolic convection-diffusion equation. Numer. Math. Theor. Methods Appl. 1, 214–234 (2008) 22. Stynes, M.: Steady-state convection-diffusion problems. Acta Numer. 14, 445–508 (2005) 23. Wesseling, P.: Uniform convergence of discretization error for a singular perturbation problem. Numer. Methods Partial Differ. Equ. 12, 657–671 (1996)
Tailored Finite Point Method for a Singular Perturbation Problem with Variable Coefficients in Two Dimensions Houde Han · Zhongyi Huang
Received: 1 December 2008 / Revised: 11 March 2009 / Accepted: 24 March 2009 © Springer Science+Business Media, LLC 2009
Abstract In this paper, we propose a tailored-finite-point method for a type of linear singular perturbation problem in two dimensions. Our finite point method has been tailored to some particular properties of the problem. Therefore, our new method can achieve very high accuracy with very coarse mesh even for very small ε, i.e. the boundary layers and interior layers do not need to be resolved numerically. In our numerical implementation, we study the classification of all the singular points for the corresponding degenerate first order linear dynamic system. We also study some cases with nonlinear coefficients. Our tailored finite point method is very efficient in both linear and nonlinear coefficients cases. Keywords Tailored finite point method · Singular perturbation problem · Singular point
1 The Model Problem We are interested in the numerical simulation of singular perturbation problems with variable coefficients. In this paper, we shall focus on problems in two dimensions. Let us consider the following problem: ∂uε ∂uε (x) + q(x) (x) = 0, ∂x1 ∂x2 ∂uε = g, ∂n
−εuε (x) + p(x) uε |D = f,
x ∈ ,
N
H. Han was supported by the NSFC Project No. 10471073. Z. Huang was supported by the NSFC Project No. 10676017, the National Basic Research Program of China under the grant 2005CB321701. H. Han · Z. Huang () Dept. of Mathematical Sciences, Tsinghua University, Beijing 100084, China e-mail: [email protected] H. Han e-mail: [email protected]
(1.1) (1.2)
J Sci Comput
where ⊂ R2 is a bounded domain. ε is a small positive parameter, p and q are two piecewise smooth functions given in domain . When 1 the solution may have some boundary layers and some interior layers in downwind directions. These layers are characterized by rapid transitions in the solution, and are thus difficult to capture in a numerical approximation without using a large number of unknowns. Also, such layers tend to cause spurious oscillations in a numerical solution to the problem. Methods for the numerical solution of problems such as (1.1)–(1.2) in bounded or unbounded domains that attempt to deal with these difficulties have been developed by many mathematicians, see e.g., [1, 2, 4, 5, 9, 10, 12, 13, 16–18, 20–23]. If ε 1, in the ‘outer domain’, i.e. in the domain out of the boundary layers and the interior layers, the second order term ‘−εuε (x)’ can be dropped from the equation. One obtains the first order partial differential equation p(x)
∂uε ∂uε (x) + q(x) (x) = 0. ∂x1 ∂x2
The characteristic curves of this equation are given by the solution of the following ordinary differential system dx2 dx1 = p(x), = q(x). (1.3) dt dt It is well-known that there is a complete classification of the singular points of the ordinary differential system (1.3) when p and q are linear functions, p(x) = ax1 + bx2 and q(x) = cx1 + dx2 [11], i.e. dx2 dx1 = ax1 + bx2 , = cx1 + dx2 , dt dt with a, b, c, d ∈ R. The equation a − λ b D(λ) = = λ2 − (a + d)λ + ad − bc = 0 c d − λ
(1.4)
is called the characteristic equation of (1.4). Let α = −(a + d), β = ad − bc, then D(λ) = λ2 + αλ + β = 0, it has two roots −α ± α 2 − 4β λ1,2 = . 2 We have the following classification of the singular points: • • • • • • • • •
Saddle Point: β < 0; Stable Node: β > 0, α > 0, α 2 − 4β > 0; Stable Focus: β > 0, α > 0, α 2 − 4β < 0; Unstable Node: β > 0, α < 0, α 2 − 4β > 0; Unstable Focus: β > 0, α < 0, α 2 − 4β < 0; Critical Stable Node/Degenerate Stable Node: β > 0, α > 0, α 2 − 4β = 0; Critical Unstable Node/Degenerate unstable Node: β > 0, α < 0, α 2 − 4β = 0; Center: β > 0, α = 0; Infinite Singularities: β = 0.
In this paper, we propose an effective numerical method, called the “tailored finite point method”, to solve the problem (1.1)–(1.2) numerically. Especially, we shall use our method
J Sci Comput
to simulate the problem with different ‘singular points’. We will also show some numerical results for the case in which the coefficients p(x) and q(x) are nonlinear functions.
2 Finite Point Method In this section we present our new approach to construct a discrete scheme for the (1.1). The new scheme is called a “a tailored finite point method”, because the finite point method has been tailored to some particular properties or solutions of the problem. The finite point method [3, 14, 15, 19] is a development of finite difference method, in which the meshless technique is emphasized. For one dimensional singular perturbation problem, the method is very close to the method of “exponential fitting” discussed in [1, 10, 16] and [20] (Chap. I, Sect. 2). We have presented the tailored finite-point schemes for the Laplace equation and a singular perturbation problem with constant coefficients in [7], we also presented the TFPM for Helmholtz equation in heterogeneous media [6] and interface problems [8]. Here we shall construct a tailored-five-point scheme for the singular perturbation equation (1.1). First, we approximate the coefficients p(x) and q(x) by piecewise constant functions, i.e. we assume that p(x) ≈ p(x 0 ) ≡ p0 and q(x) ≈ q(x 0 ) ≡ q0 in the domain 0 (cf. Fig. 1). Let p0 x1 q0 x2 uε (x1 , x2 ) = v ε (x1 , x2 ) exp + . (2.1) 2ε 2ε Then v ε satisfies: −v ε +
p02 + q02 ε v = 0. 4ε 2
(2.2)
Let α0 =
p0 , 2ε
β0 =
q0 , 2ε
μ0 =
α02 + β02 .
(2.3)
The solution v ε of (2.2) can be expanded at x = x 0 as: v ε (x1 , x2 ) = a0 I0 (μ0 r) +
∞ n=1
Fig. 1 The location mesh of points x 0 , x 1 , . . . , x 9
In (μ0 r) (an cos nθ + bn sin nθ ) ,
(2.4)
J Sci Comput
Fig. 2 The numerical solutions of TFPM for Example 3.2 (saddle point): h = 1/32
where (r, θ ) is the ordered pair of polar coordinates of x = (x1 , x2 ) with the pole at x 0 and In is the n-th order modified Bessel function of first kind. Therefore the solution uε of (1.1) can be expanded at x = x 0 as: uε (x1 , x2 ) = b00 + er(α0 cos θ+β0 sin θ) ∞ 0
0 0 × a0 I0 (μ0 r) + In (μ0 r) an cos nθ + bn sin nθ .
(2.5)
n=1
We take the first few terms in (2.5) to construct an approximation u |0 ≈ ε
b00
+e
r(α0 cos θ+β0 sin θ)
a00 I0 (μ0 r) +
N n=1
In (μ0 r)
an0 cos nθ
+ bn0 sin nθ
.
(2.6)
J Sci Comput
Fig. 3 The numerical solutions of TFPM for Example 3.3 (stable node): h = 1/32
Now, we define the tailored 5-point scheme as follows. Set
W 4 = u(r, θ ) u = b00 + er(α0 cos θ+β0 sin θ) a00 I0 (μ0 r) + I1 (μ0 r) a10 cos θ + b10 sin θ , ∀an0 , bn0 ∈ R, n = 0, 1 . Consider μ0 > 0. First, suppose √ |p0 | ≥ ( 2 + 1)|q0 |,
or
√ |q0 | ≥ ( 2 + 1)|p0 |.
(2.7)
In this case, the tailored finite point equation uses the points x 0 , x 1 , x 2 , x 3 , x 4 . The equation at x 0 is 4 n=0
cn u(x n ) = 0
(2.8)
J Sci Comput
where the numbers cn (n = 0, 1, 2, 3, 4), satisfy 4
cn w(x n ) = 0,
∀w ∈ W 4 .
(2.9)
n=0
Let αh = α0 h, βh = β0 h, and
C1 (μ0 h) = 2I0 (μ0 h) cosh(αh ) − cosh(βh ) .
(2.10)
After some algebra we arrive at the formulas c1 = −
c0 e−αh I0 (μ0 h) − cosh(βh ) , C1 (μ0 h)
c0 e−βh cosh(αh ) − I0 (μ0 h) , C1 (μ0 h) c0 eαh I0 (μ0 h) − cosh(βh ) , c3 = − C1 (μ0 h) c2 = −
c4 = −
c0 eβh cosh(αh ) − I0 (μ0 h) . C1 (μ0 h)
(2.11) (2.12) (2.13) (2.14)
If (2.7) does not hold, the tailored finite point equation uses the points x 0 , x 5 , x 6 , x 7 , x 8 . The equation at x 0 in this case is c0 u(x 0 ) +
8
cn u(x n ) = 0
(2.15)
n=5
where the numbers cn (n = 0, 5, 6, 7, 8), satisfy c0 w(x 0 ) +
8
cn w(x n ) = 0,
∀w ∈ W 4 .
(2.16)
n=5
Let
√
C2 (μ0 h) = 2I0 ( 2μ0 h) cosh(αh + βh ) − cosh(αh − βh ) .
(2.17)
After some algebra we then arrive at the formulas e−αh −βh √ I0 ( 2μ0 h) − cosh(αh − βh ) , C2 (μ0 h) √ eαh −βh cosh(αh + βh ) − I0 ( 2μ0 h) , c6 = −c0 C2 (μ0 h) eαh +βh √ I0 ( 2μ0 h) − cosh(αh − βh ) , c7 = −c0 C2 (μ0 h) √ e−αh +βh cosh(αh + βh ) − I0 ( 2μ0 h) . c8 = −c0 C2 (μ0 h) c5 = −c0
(2.18) (2.19) (2.20) (2.21)
Equations (2.8) and (2.15) give the 5 point tailored finite point scheme at the point x 0 for the problem (1.1).
J Sci Comput
Fig. 4 The numerical solutions of TFPM for Example 3.4 (stable focus): h = 1/32
We can also obtain a tailored 9-point scheme in the similar way. Setting 3 0
0 r(α0 cos θ+β0 sin θ) 0 0 a0 I0 (μ0 r) + W = u(r, θ ) u = b0 + e In (μ0 r) an cos nθ + bn sin nθ ,
8
n=1
∀an0 , bn0 ∈ R, n = 0, 1, 2, 3 , we seek numbers cn (n = 0, . . . , 8), such that 8 n=0
cn w(x n ) = 0,
∀w ∈ W 8 .
(2.22)
J Sci Comput
Fig. 5 The numerical solutions of TFPM for Example 3.5 (center): h = 1/32
Let √
C3 (μ0 h) = 2I0 ( 2μ0 h) cosh(αh ) + cosh(βh ) − 2I0 (μ0 h) cosh(αh − βh ) + cosh(αh + βh ) . (2.23) After some algebra we arrive at the tailored 9-point scheme √ c0 e−αh 2I0 ( 2μ0 h) − cosh(αh + βh ) − cosh(αh − βh ) , 2C3 (μ0 h) √ c0 e−βh 2I0 ( 2μ0 h) − cosh(αh + βh ) − cosh(αh − βh ) , c2 = − 2C3 (μ0 h) √ c0 eαh 2I0 ( 2μ0 h) − cosh(αh + βh ) − cosh(αh − βh ) , c3 = − 2C3 (μ0 h) √ c0 eβh 2I0 ( 2μ0 h) − cosh(αh + βh ) − cosh(αh − βh ) , c4 = − 2C3 (μ0 h)
c1 = −
J Sci Comput
Fig. 6 The numerical solutions of TFPM for Example 3.6 (unstable node): h = 1/32
c5 = −
c0 e−αh −βh −2I0 (μ0 h) + cosh(αh ) + cosh(βh ) , 2C3 (μ0 h)
c6 = −
c0 eαh −βh −2I0 (μ0 h) + cosh(αh ) + cosh(βh ) , 2C3 (μ0 h)
c7 = −
c0 eαh +βh −2I0 (μ0 h) + cosh(αh ) + cosh(βh ) , 2C3 (μ0 h)
c8 = −
c0 e−αh +βh −2I0 (μ0 h) + cosh(αh ) + cosh(βh ) . 2C3 (μ0 h)
For μ0 = 0, i.e. p0 = q0 = 0, (1.1) is reduced to the Laplace equation. Then we should use the standard five-point finite difference scheme or nine-point finite difference scheme for Laplace equation, i.e. u1 + u2 + u3 + u4 − 4u0 = 0, h2
J Sci Comput
Fig. 7 The numerical solutions of TFPM for Example 3.7 (unstable focus): h = 1/32
or 4(u1 + u2 + u3 + u4 ) + (u5 + u6 + u7 + u8 ) − 20u0 = 0. h2
3 Numerical Implementation In this section, we give some numerical examples to show the efficiency of our new method. In the following tables and figures, the ‘exact’ solution u is obtained numerically on a very fine mesh. uFh P 5 is the solution solved by our five point scheme. Furthermore, we let
EhF P ,2 =
uFh P 5 − u L2 () , u L2 ()
EhF P ,ε =
uFh P 5 − u ε, , u ε,
F P 5/9 = the maximum nodal errors of 5-point TFPM or 9-point TFPM, E∞
J Sci Comput
Fig. 8 The numerical solutions of TFPM for Example 3.8 (critical stable node): h = 1/32
where u 2ε, = ε|u|21, + u 2L2 () , and
|u|21,
≡
|∇u(x)|2 dx1 dx2 ,
u 2L2 ()
≡
|u(x)|2 dx1 dx2 .
Example 3.1 (Accuracy test) First, we consider the following problem: −εu + pux1 + qux2 = 0, u|x1 =−1 = x2 + 1, u|x2 =1 = 1 − x1 ,
∀x ∈ ,
u|x2 =−1 = x1 + 1,
(3.1) u|x1 =1 = 1 − x2 ,
(3.2)
J Sci Comput Table 1 Example 3.1: Relative errors of the numerical solutions for TFPM with = 0.01 Mesh size h
1/8
1/16
1/32
1/64
EhF P ,2
3.561E−1
9.299E−2
2.389E−2
5.973E−3
1.9
2.0
2.0
1.214E−1
3.096E−2
7.815E−3
2.0
2.0
2.0
Convergence order EhF P ,ε
4.791E−1
Convergence order
Table 2 Example 3.1: The nodal errors of the numerical solutions by TFPM, h = 1/32
0.5
0.1
0.01
0.001
FP5 E∞
5.741E−2
5.742E−2
5.846E−2
6.024E−2
FP9 E∞
3.473E−2
3.467E−2
3.410E−2
3.439E−2
where
= x = (x1 , x2 ) ∈ R2 −1 < x1 < 1, −1 < x2 < 1 ,
(3.3)
and p(x) = x1 + x2 ,
q(x) = x1 − x2 .
The relative errors of our approximations are shown in Table 1. For our five-point TFPM, it seems that we have the second order convergence rate for L2 norm and energy norm. These convergent rates coincide with the corresponding results in 1D which were proved in [8]. From Table 2, we can find that our TFPM has a uniformly nodal errors as → 0. In what follows, we show numerical results for all kinds of singular points of the dynamic system (1.4). We consider different choices of the coefficients p and q, including cases when p and q are nonlinear functions of x. Example 3.2 (Saddle Point) We consider the problem (3.1)–(3.3) with p(x) = x2 ,
q(x) = x1 .
In this case, we have two real eigenvalues with opposite sign for the linear dynamic system (1.4). The results are shown in Fig. 2. We can see that, although the solution has very sharp boundary layers as → 0, our TFPM can capture them in very coarse mesh without any pseudo-oscillation. Example 3.3 (Stable Node) We consider the problem (3.1)–(3.3) with p(x) = −x1 + 0.5x2 ,
q(x) = 0.5x1 − x2 .
In this case, we have two negative real eigenvalues for the linear dynamic system (1.4). The results are shown in Fig. 3. In this case, the solution has four interior layers, and the origin is a singular point. Our TFPM also can detect these results in coarse mesh.
J Sci Comput
Fig. 9 The numerical solutions of TFPM for Example 3.9 (degenerate stable node): h = 1/32
Example 3.4 (Stable Focus) We consider the problem (3.1)–(3.3) with p(x) = −x1 − x2 ,
q(x) = x1 − x2 .
In this case, we have two conjugated complex eigenvalues with negative real parts for the linear dynamic system (1.4). The results are shown in Fig. 4. Here the solution has four curved interior layers. Our TFPM can get them very clearly. Example 3.5 (Center) We consider the problem (3.1)–(3.3) with p(x) = x2 ,
q(x) = −x1 .
In this case, we have two conjugated imaginary eigenvalues for the linear dynamic system (1.4). The results are shown in Fig. 5.
J Sci Comput
Fig. 10 The numerical solutions of TFPM for Example 3.10 (critical unstable node): h = 1/32
Example 3.6 (Unstable Node) We consider the problem (3.1)–(3.3) with p(x) = x1 + 0.5x2 ,
q(x) = 0.5x1 + x2 .
In this case, we have two positive real eigenvalues for the linear dynamic system (1.4). The results are shown in Fig. 6. Example 3.7 (Unstable Focus) We consider the problem (3.1)–(3.3) with p(x) = x1 − x2 ,
q(x) = x1 + x2 .
In this case, we have two conjugated complex eigenvalues with positive real parts for the linear dynamic system (1.4). The results are shown in Fig. 7. From Figs. 5–7, we can see that our TFPM can capture the sharp boundary layers in very coarse mesh without any artificial oscillation.
J Sci Comput
Fig. 11 The numerical solutions of TFPM for Example 3.11 (degenerate unstable node): h = 1/32
Example 3.8 (Critical Stable Node) We consider the problem (3.1)–(3.3) with p(x) = −x1 ,
q(x) = −x2 .
In this case, we have two double negative real eigenvalues with single elementary divisor for the linear dynamic system (1.4). The results are shown in Fig. 8. In this case and Example 3.9 (cf. Fig. 9), the solutions have four interior layers when → 0. Our TFPM can obtain them very robustly. Especially, the contour curves are very smooth even when h . Example 3.9 (Degenerate Stable Node) We consider the problem (3.1)–(3.3) with p(x) = −x1 ,
q(x) = x1 − x2 .
In this case, we have two double negative real eigenvalues with double elementary divisor for the linear dynamic system (1.4). The results are shown in Fig. 9.
J Sci Comput
Fig. 12 The numerical solutions of TFPM for Example 3.12: h = 1/32
Example 3.10 (Critical Unstable Node) We consider the problem (3.1)–(3.3) with p(x) = x1 ,
q(x) = x2 .
In this case, we have two double positive real eigenvalues with single elementary divisor for the linear dynamic system (1.4). The results are shown in Fig. 10. Example 3.11 (Degenerate Unstable Node) We consider the problem (3.1)–(3.3) with p(x) = −x1 ,
q(x) = x1 − x2 .
In this case, we have two double positive real eigenvalues with double elementary divisor for the linear dynamic system (1.4). The results are shown in Fig. 11. Example 3.12 We consider the problem (3.1)–(3.3) with p(x) = x1 + x2 ,
q(x) = 0.
J Sci Comput
Fig. 13 The numerical solutions of TFPM for Example 3.13: h = 1/32
In this case, we have a line x1 + x2 = 0 of singular points for the linear dynamic system (1.4). The results are shown in Fig. 12. Example 3.13 We consider the problem (3.1)–(3.3) with p(x) = x1 − x2 ,
q(x) = −x1 + x2 .
In this case, we have a line x1 − x2 = 0 of singular points for the linear dynamic system (1.4). The results are shown in Fig. 13. In Examples 3.10–3.13, the solutions have very sharp boundary layers. But our TFPM can get them without any oscillation in very coarse mesh. Example 3.14 We consider the problem (3.1)–(3.3) with nonlinear coefficients p(x) = x1 x2 ,
q(x) = x22 − x14 .
J Sci Comput
Fig. 14 The numerical solutions of TFPM for Example 3.14: h = 1/32
In this case, the origin O is a high order singular point for the nonlinear dynamic system (1.3). The results are shown in Fig. 14. Here the solution has some boundary layers and interior layers. Our TFPM can get all of them with high accuracy in coarse mesh. Example 3.15 We consider the problem (3.1)–(3.3) with nonlinear coefficients p(x) = x1 x2 − 3x13 ,
q(x) = x22 − 6x12 x2 + x14 .
In this case, the origin O is also a high order singular point. The results are shown in Fig. 15.
Example 3.16 We consider the problem (3.1)–(3.3) with nonlinear coefficients p(x) = x13 x2 − x1 x23 ,
q(x) = x12 x22 − x24 + x16 .
J Sci Comput
Fig. 15 The numerical solutions of TFPM for Example 3.15: h = 1/32
In this case, the origin O is also a high order singular point. The results are shown in Fig. 16. From Figs. 15–16, we can see that our TFPM has high accuracy in coarse mesh even the solution has very complex structures. 4 Conclusion In this paper, we present a tailored-finite-point method for a number of linear singular perturbation problems with variable coefficients in two dimensional case. First, we approximate the coefficients with piecewise constants. Then we use the eigenfunction expansion of the local problem with constant coefficients to get the tailored scheme. We consider all kinds of the singular points of the reduced first order linear differential equation, including the saddle point, the stable/unstable node, the stable/unstable focus, etc. Furthermore, we also consider the nonlinear cases. All of our numerical results show that our new method can achieve very high accuracy with very coarse mesh whenever ε is small or large, i.e. we can achieve high accuracy even when our mesh size is much larger than the thickness of the boundary layers and the interior layers. By contrast, as we know, the traditional finite element method
J Sci Comput
Fig. 16 The numerical solutions of TFPM for Example 3.16: h = 1/32
or finite difference method can not get satisfactory numerical results with the same mesh, especially for small ε. We will study the error estimates of our TFPM soon. Acknowledgement
The authors thank Prof. Bruce Kellogg for fruitful discussions on this work.
References 1. Berger, A.E., Han, H.D., Kellogg, R.B.: A priori estimates and analysis of a numerical method for a turning point problem. Math. Comput. 42, 465–492 (1984) 2. Brayanov, I., Dimitrova, I.: Uniformly convergent high-order schemes for a 2D elliptic reaction-diffusion problem with anisotropic coefficients. Lect. Notes Comput. Sci. 2542, 395–402 (2003) 3. Cheng, M., Liu, G.R.: A novel finite point method for flow simulation. Int. J. Numer. Methods Fluids 39, 1161–1178 (2002) 4. Ge, L., Zhang, J.: High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids. J. Comput. Phys. 171, 560–578 (2001) 5. Han, H.D.: The artificial boundary method—numerical solutions of partial differential equations on unbounded domains. In: Li, T., Zhang, P. (eds.) Frontiers and Prospects of Contemporary Applied Mathematics. Higher Education Press/World Scientific, Singapore (2005)
J Sci Comput 6. Han, H., Huang, Z.: A tailored finite point method for the Helmholtz equation with high wave numbers in heterogeneous medium. J. Comput. Math. 26, 728–739 (2008) 7. Han, H., Huang, Z., Kellogg, B.: A Tailored finite point method for a singular perturbation problem on an unbounded domain. J. Sci. Comput. 36, 243–261 (2008) 8. Huang, Z.: Tailored finite point method for the interface problem. Netw. Heterogeneous Media 4, 91–106 (2009) 9. Hemker, P.W.: A singularly perturbed model problem for numerical computation. J. Comput. Appl. Math. 76, 277–285 (1996) 10. Il’in, A.M.: Differencing scheme for a differential equation with a small parameter affecting the highest derivative. Math. Notes 6, 596–602 (1969) 11. Iooss, G., Joseph, D.D.: Elementary Stability and Bifurcation Theory. Springer, New York (1980) 12. Li, J., Chen, Y.: Uniform convergence analysis for singularly perturbed elliptic problems with parabolic layers. Numer. Math. Theor. Methods Appl. 1, 138–149 (2008) 13. Li, J., Navon, I.M.: Uniformly convergent finite element methods for singularly perturbed elliptic boundary value problems: convection-diffusion. Comput. Methods Appl. Mech. Eng. 162, 49–78 (1998) 14. Lin, H., Atluri, S.N.: The Meshless Local Petrov-Galerkin (MLPG) method for solving incompressible Navier-Stokes equations. CMES 2, 117–142 (2001) 15. Mendeza, B., Velazquez, A.: Finite point solver for the simulation of 2-D laminar incompressible unsteady flows. Comput. Methods Appl. Mech. Eng. 193, 825–848 (2004) 16. Miller, J.J.H.: On the convergence, uniformly in ε, of difference schemes for a two-point boundary singular perturbation problem. In: Hernker, P.W., Miller, J.J.H. (eds.) Numerical Analysis of Singular Perturbation Problems, pp. 467–474. Academic Press, San Diego (1979) 17. Morton, K.W.: Numerical Solution of Converction-Diffusion Problems. Applied Mathematics and Mathematical Computation, vol. 12. Chapman and Hall, London (1996) 18. Morton, K.W., Stynes, M., Süli, E.: Analysis of a cell-vertex finite volume method for convectiondiffusion problems. Math. Comput. 66, 1389–1406 (1997) 19. Oñate, E., Idelsohn, S., Zienkiewicz, O.C., Taylor, R.L.: A finite point method in computational mechanics. Applications to convective transport and fluid flow. Int. J. Numer. Methods Eng. 39, 3839–3866 (1996) 20. Roos, H.-G., Stynes, M., Tobiska, L.: Numerical Methods for Singularly Perturbed Differential Equations. Springer, New York (1996) 21. Shishkin, G.I.: A finite difference scheme on a priori adapted meshes for a singularly perturbed parabolic convection-diffusion equation. Numer. Math. Theor. Methods Appl. 1, 214–234 (2008) 22. Stynes, M.: Steady-state convection-diffusion problems. Acta Numer. 14, 445–508 (2005) 23. Wesseling, P.: Uniform convergence of discretization error for a singular perturbation problem. Numer. Methods Partial Differ. Equ. 12, 657–671 (1996)
