NUMERICAL SOLUTION OF A TURNING POINT PROBLEM ^{∗}

WEI PAI TANG^{†}

Abstract.

The turning point problem

\{ \begin{matrix}−εΔu + x\frac{∂u}{∂x} + y \frac{∂u}{∂y} = 0 & (x, y) ∈ [(−1, 1) × (−1, 1)] \\ u(−1, y) = V_{a},  u(1, y) = V_{b}, \\ u(x, −1) = V_{c},  u(x, 1) = V_{d}, \end{matrix}

is known to have some extremely small eigenvalues. No successful numerical solution to this problem has been reported. In this paper, a numerical procedure is proposed. All four boundary layers are well defined and the numerical singularity is successfully removed.

Key words. boundary layer, domain decomposition, overlap, Schwarz Alternating Method(SAM), turning point problem

AMS(MOS) subject classifications. 65F10, 65N20

1. Introduction. The singularly perturbed elliptic partial differential equation

− εΔu + a(x, y)\frac{∂u}{∂x} + b(x, y)\frac{∂u}{∂y} + k(x, y)u = 0 (1)

has been extensively studied by many researchers. When ε → 0, the solution to this problem becomes difficult. A typical property of the solution is the existence of boundary layers. This problem has many important applications including the solution of the Navier-Stokes equations and stochastic differential equations [13].

In (1), the lower order operator represents the deterministic flow field while the second order part represents a slow diffusion of particles. Therefore, the results will depend on the nature of the underlying flow. In [13], Matkowsky classifies the singularly perturbed elliptic boundary value problem into three cases according to how the particles are diffusing (see Fig. 1):

∗ This research was supported by the Natural Sciences and Engineering Research Council of Canada.

† Department of Computer Science, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1.

Fig. 1.

Type (1). with a flow, Type (2). across a flow, Type (3). against a flow. In particular, the problem of type (3) is called the turning point problem, since

the flow changes its direction inside of the region. Asymptotic analysis of the three types of diffusion can be traced back to the early 1950’s [10, 11]. The third type of diffusion is the most difficult one. The first few results were obtained by O’Malley; Ventsel and Freidlin [14, 20]. Later, a stronger result was reported by Ludwig [12]. More results have been published since then [7, 15].

The investigation of the numerical solution of the turning point problem started from stiff-ODE problems. Dorr (1971) first reported that the turning point problem in the one-dimensional case is extremely ill-conditioned [4]. For example, the condition number of the matrix equation for the following problem

\{ \begin{matrix}−εy′′ + xy = 0, & x ∈ (−1, 1) \\ y(−1) = a, & y(1) = b \end{matrix}

is κ(A) = 1.1021 · 10^{14}, where ε = 0.003 and n = 100. Many numerical techniques for the stiff-ODE and turning point problems in the one-dimensional case are found in [1]. However, successful numerical techniques for turning point problems in higher dimensional spaces have not been reported. This is because of the extreme numerical singularity of the resulting matrix equation. B. Zhu has shown [21] that the smallest eigenvalue of the following problem

\{ \begin{matrix}−εΔu + x\frac{∂u}{∂x} + y \frac{∂u}{∂y} = λu & (x, y) ∈ [(−1, 1) × (−1, 1)] \\ u(x, y)|_{Γ} = 0 \end{matrix}

has the estimate

|λ| ≤ O(ε^{−1/2}e^{−α/ε})

where α > 0.

Much work has been done for boundary layer problems in higher dimensional space. For example, Hedstrom and Osterheld first studied the effect of ε on the boundary layer [6]. Segal discussed the different aspects of numerical methods which relate to the computation for singular perturbation problems [18]. Rodrigue and Reiter [16] investigated the application of the Schwarz Alternating Method (SAM) to (1). Brown, Chin, Hedstrom, Manteuffel and Scroggs [2, 3, 9, 17] studied other domain decomposition techniques. Elman and Golub reported their sequence of studies on iterative methods for the convection-diffusion problem [5].

In this paper, a numerical procedure using the domain decomposition approach, or more precisely SAM, for the solution of turning point problems is presented. This procedure can remove the severe numerical singularity in the original form; therefore, a successful numerical solution for this problem becomes feasible. The same SAM can also be used to accurately define the sharp boundary layers.

2. Turning point problem. In this paper, the following turning point problem is considered:

\{ \begin{matrix}−εΔu + x\frac{∂u}{∂x} + y \frac{∂u}{∂y} = 0 & (x, y) ∈ [(−1, 1) × (−1, 1)] \\ u(−1, y) = V_{a},  u(1, y) = V_{b}, \\ u(x, −1) = V_{c},  u(x, 1) = V_{d}, \end{matrix} (2)

This problem was suggested by C. Holland [8]. It is typical of a wide class of turning point problems. The solution to this problem has four sharp boundary layers of width ε at each side of the square and the solution at the origin is

u(0, 0) = \frac{1}{4}(V_{a} + V_{b} + V_{c} + V_{d}).

The technique we use for solving this problem can also be applied to other cases [19].

Applying Zhu’s result to this problem, the linear system of equations discretized from (2) will have an extremely ill-conditioned matrix if ε is small. No meaningful numerical solution can be obtained if special techniques are not used [8]. However, if we carefully examine the relationship between ε and the condition number of the matrix, the following three observations are very important to the construction of our numerical procedure for the solution of (2).

1. The smallness of ε is relative to the size of the solution domain. If we reduce the size of the solution region in (2) the same small ε will result in a different condition number. It is clear that the solution for a smaller region is less difficult. Thus, if a domain decomposition approach is used, the solution of each subproblem will be easier.

2. When the domain decomposition approach is applied to the turning point problem, most of the subproblems are not turning point problems! Therefore, the solution of these subproblems is not an issue. There is only one subdomain which contains a turning point^{1}. Fortunately, it can be made to

1 We only discuss the single turning point case in this paper. There is no conceptual difficulty in generalizing our procedure to the case of more than one turning point.

have a very small size.

3. The asymptotic analysis of this problem [14, 21] showed that the expansion of the solution of (2) in a series of ε has one base term. In particular, this term is independent of ε . Therefore, if ε_{1} and ε_{2} are close, the corresponding solutions of the same boundary value problem (2) are also close. Most of the changes happen around the boundary layers. Numerical computations have verified their analysis.

Based on the first two observations domain decomposition approaches (and the Schwarz alternating method (SAM) in particular) seem to be very helpful in overcoming the numerical singularity which appeared in the original problem. Specifically, if we decompose the solution region into many small overlapping subregions, there is only one subdomain which contains a turning point. The rest of the subproblems are “easier” boundary layer problems, for which many known numerical techniques can be used to solve them. Since the size of that subdomain which contains the turning point is smaller, the numerical difficulty of solving this subproblem will not be as severe as the original problem. However, this naive approach still faces the stability problem if a poor initial guess is given for a very small ε . Fortunately, the third observation can lead us to a successful solution process. We start from a “large” ε for which the turning point problem can be solved without a stability problem. Then we reduce the ε and decompose the solution region into several overlapping subregions. Using the solution from the large ε as the initial guess for the new decomposition, SAM can be applied to obtain the solution for the new smaller ε . Two issues in this process are important:

• The new ε should not change too rapidly so that the solution corresponding to the new problem has no large changes in parts of the region.

• The size of the subdomain which contains the turning point has to be small enough such that the solution of this problem has no stability problems.

3. Convergence analysis. For problem (2), the matrices resulting from many finite difference or finite element discretizations are unsymmetric. The convergence of the SAM for these matrices needs to be justified.

It is known that the matrix from the discretization of (2) is diagonally dominant, if an upwind scheme is used. When the grid size h is small enough, the matrix resulting from a central scheme is also diagonally dominant [5]. With this diagonal dominance condition, a convergence result can be shown.

Let

Ax = b

be the matrix equation from the discretization of (2), and Ω_{i}, i = 1, ⋯, k be the overlapping subdomains such that all grid nodes on the artificial boundary are located at least in the interior of one subdomain. Denote x_{Ω_{i}} as the unknown vector in subdomain Ω_{i} and and A_{Ω_{i}} the corresponding principal submatrix for these unknowns. The description of the Schwarz alternating method applied for this particular decomposition is given in Fig. 3.

Procedure Turning point problem () Choose ε_{0} such that

\{ \begin{matrix}−ε_{0}Δu^{(0)} + x\frac{∂u^{(0)}}{∂x} + y \frac{∂u^{(0)}}{∂y} = 0 & (x, y) ∈ [(−1, 1) × (−1, 1)] \\ u^{(0)}(−1, y) = V_{a},  u^{(0)}(1, y) = V_{b}, \\ u^{(0)}(x, −1) = V_{c},  u^{(0)}(x, 1) = V_{d}, \end{matrix}

can be solved without stability problem.

Let i = 0

While ε_{i} > ε do Pick ε_{i+1} < ε_{i} and a new set of overlapping subdomains. Use u^{(i)} as initial guess start the SAM iteration for

\{ \begin{matrix}−ε_{i+1}Δu^{(i+1)} + x\frac{∂u^{(i+1)}}{∂x} + y \frac{∂u^{(i+1)}}{∂y} = 0 \\ (x, y) ∈ [(−1, 1) × (−1, 1)] \\ u^{(i+1)}(−1, y) = V_{a},  u^{(i+1)}(1, y) = V_{b}, \\ u^{(i+1)}(x, −1) = V_{c},  u^{(i+1)}(x, 1) = V_{d}, \end{matrix}

until it converges.

i = i + 1 end

end

Fig. 2. Numerical procedure for turning point problem.

Given initial guess x^{(0)},

While ∥Ax^{(l)} − b∥ ≥ T ol do x^{(l+1)} = x^{(l)} For i = 1 to k

Solve^{a}
A_{Ω_{i}}x_{Ω_{i}}^{(l+1)} = b_{Ω_{i}}^{(l+1)},
Update x^{(l+1)} by x_{Ω_{i}}^{(l+1)}

end l = l + 1

end

a Note that b_{Ω_{i}}^{(l+1)} contains the information from the artificial boundary as well.

Fig. 3. Description of the SAM algorithm.

Then the following lemma can lead to the convergence of the SAM algorithm in this case.

Lemma 3.1. If the matrix A_{Ω_{i}} is diagonally dominant, then

∥x_{Ω_{i}}∥_{∞} < γ∥x_{Γ_{i}}∥_{∞}

where x_{Γ_{i}} is a vector which contains all boundary node values (including both true and artificial boundary nodes) of Ω_{i} and γ < 1

Proof. We prove this lemma by contradiction.

4. Boundary layers. To gain efficiency and to maintain the diagonal dominance, the upwind scheme is used for defining the overall solution. However, it is known that this scheme is diffusive on the boundary layer. The central scheme is used around the boundary layer to improve accuracy. Many much smaller overlapping subregions (say size of 30ε × 30ε ) along the boundary are allocated. Of course, the grid size of these subdomains is refined to the size of ε . Apply SAM to these new boundary subdomains and a convergence solution is ensured. Our numerical results indicate the improvement of sharpness in the boundary layers.

5. Numerical testing. Numerical testing is carried out for this model turning problem (2). We start with ε_{0} = 0.02. The solution of this problem has no stability problem. Then we set ε_{k} = 0.5∗ε_{k−1}; decompose the square into (k+1)×(k+1) equal sized subdomains and apply SAM to the turning point problem with the smaller ε_{k} . It is interesting that the solution of the turning point subproblem does not change any further after k > 3. In a real application, many of the subproblems which are away from the boundary layers need not be recomputed after reducing ε a few times. This phenomenon can be explained by the asymptotic expansion for this problem [14, 21]. We present the surface plots of the solution for ε = 0.0025, 0.00125, and 0.000625. In particular, we prresent both solutions before and after using centeral scheme to sharpen the boundary layers. Our computation stoped when ε_{k} = 1/1600. If the computation would continue for even smaller ε , the decompostion strategy we used here need to be updated to ensure the stability problem of the turning point subproblem.

Acknowledgment. The author would like to thank the helpful discussion with C. Holland, G. Rodrigue and B. Zhu.

REFERENCES

[1] U. Ascher, R. Mattheij, and R. Russell, Numerical solution of boundary value problems for ordinary differential equations, Prentice Hall, New Jersey, 1988.

[2] D. L. Brown, R. C. Y. Chin, G. Hedstrom, and T. A. Manteuffel, Layer tracking, asymptotics, and domain decomposition, Research Report UCRL-JC-106336, Lawrence Livermore National Lab, 1991.

[3] R. C. Y. Chin and T. A. Manteuffel, An analysis of block successive overrelaxation for a class of matrices with complex spectra, SIAM J. Num. Anal., 25 (1988), pp. 564–585.

[4] F. W. Dorr, An example of ill conditioning in the numerical solution of singular perturbation problems, Mathematics Of Computation, 25 (1971), pp. 271–283.

Fig. 4. Surface plots for ε = 1/50, 1/1600, 1/200. 7 Fig. 5. Surface plots for ε = 1/400, 1/800.

8 [5] H. C. Elman and G. H. Golub, Iterative methods for cyclically reduced non-selfadjoint

linear systems, Mathematics Of Computation, 54 (1990), pp. 671–700.

[6] G. W. Hedstrom and A. Osterheld, The effect of cell Reynolds number on the computation of a boundary layer, Journal of Computational Physics, 37 (1980), pp. 399–421.

[7] C. J. Holland, Singular perturbations in elliptic boundary value problems, Journal of Differential Equations, 20 (1976), pp. 248–265.

[8] , Comments on singularly perturbed elliptic partial equations, (1991). (Personal communication).

[9] F. A. Howes and G. W. Hedstrom, A domain decomposition methods for a convectiondiffusion equation with turning points, in Domain Decomposition Methods, T. Chan, R. Glowinski, J. Periaux, and O. Widlund, eds., SIAM, 1989, pp. 38–46.

[10] R. A. Khasminski1, On diffusion processes with a small parameter, Izv. Akad. Nauk. SSSR Ser. Mat., 27 (1963), pp. 1281–1300.

[11] N. Levinson, The first boundary value problem for εδu + au_{x} + bu_{y} + cu = d for small ε, Ann. Math., 51 (1950), pp. 428–445.

[12] D. Ludwig, Persistence of dynamical systems under random perturbations, SIAM Review, 17 (1975), pp. 605–640.

[13] B. J. Matkowsky, Singular perturbations, stochastic differential equations, and applications, in Singular Perturbations and Asymptotics, Academic Press, 1980, pp. 109–147.

[14] R. E. O’Malley, On boundary value problems for a singular perturbed differential equation with a turning point, Journal of Differential Equations, 20 (1976), pp. 248–265.

[15] , A singular singularly perturbed linear boundary value problem, SIAM Journal On Mathematical Analysis, 10 (1979), pp. 695–708.

[16] G. Rodrigue and E. Reiter, A domain decomposition method for boundary layer problems, in Domain Decomposition Methods, T. Chan, R. Glowinski, J. Periaux, and O. Widlund, eds., SIAM, 1989, pp. 226–234.

[17] J. S. Scroggs, A physically motivated domain decomposition for singularly perturbed equations, SIAM J. Num. Anal., 28 (1991), pp. 168–178.

[18] A. Segal, Aspects of numerical methods for elliptic singular perturbation problems, SIAM J. Sci. Statist. Comp., 3 (1982), pp. 327–348.

[19] W.-P. Tang and B. Zhu, Numerical solutions of a turning point problem and its conjugate problem, research report, Department of Computer Science, University of Waterloo, Waterloo, Ontario, Canada, 1991. (in preparation).

[20] A. D. Ventsel and M. I. Friedlin, On small random perturbation of dynamical systems, Russian Math. Survey, 25 (1970), pp. 1–55.

[21] B. Zhu, Multi-dimensional boundary layer problems with presence of Ackerberg-O’Malley resonance, in Proceedings BAIL 5, Shanghai, 1988, pp. 436–441.