Asymptotic Expansion with Boundary Layer Analysis for Strongly Anisotropic Elliptic Equations

TL;DR

This paper develops an asymptotic expansion method with boundary layer analysis for strongly anisotropic elliptic equations, addressing ill-posedness and boundary conditions to improve solution approximation accuracy.

Contribution

It introduces a systematic approach combining mean-fluctuation decomposition and boundary layer analysis for higher-order asymptotic expansions in anisotropic elliptic problems.

Findings

Derived asymptotic expansions up to arbitrary order.

Handled ill-posedness via solution decomposition.

Provided composite approximation for the leading order.

Abstract

In this article, we derive the asymptotic expansion, up to an arbitrary order in theory, for the solution of a two-dimensional elliptic equation with strongly anisotropic diffusion coefficients along different directions, subject to the Neumann boundary condition and the Dirichlet boundary condition on specific parts of the domain boundary, respectively. The ill-posedness arising from the Neumann boundary condition in the strongly anisotropic diffusion limit is handled by the decomposition of the solution into a mean part and a fluctuation part. The boundary layer analysis due to the Dirichlet boundary condition is conducted for each order in the expansion for the fluctuation part. Our results suggest that the leading order is the combination of the mean part and the composite approximation of the fluctuation part for the general Dirichlet boundary condition.

Tables

Table 1. Table 1. The discrete L ∞ subscript 𝐿 L_{\infty} norms of the errors r 0 subscript 𝑟 0 r_{0} and r 2 subscript 𝑟 2 r_{2} for different values of ε 𝜀 \varepsilon .

${\epsilon }^2$	${\Vert r_0\Vert }_{\infty }$	${\Vert r_2\Vert }_{\infty }$
0.001	$1.0533E-04$	$5.2834E-07$
0.005	$5.2222E-04$	$7.3587E-06$
0.01	$1.0335E-03$	$2.8475E-05$
0.05	$4.7441E-03$	$5.7746E-04$
0.1	$8.6241E-03$	$1.9240E-03$

Equations

$$\{\begin{array}{ll}-{\epsilon }^{-2}{\partial }_x^2{u}_{\epsilon }(x,y)-{\partial }_y^2{u}_{\epsilon }(x,y)=f(x,y)& \text{in }D,\\ {\partial }_x{u}_{\epsilon }(0,y)={\partial }_x{u}_{\epsilon }(1,y)=0,& 0⩽y⩽1,\\ u_{\epsilon }(x,0)={\varphi }_0(x),u_{\epsilon }(x,1)={\varphi }_1(x),& 0⩽x⩽1.\end{array}$$

$$\{\begin{array}{ll}& \mathrm{d}{X}_t={\epsilon }^{-1}\mathrm{d}{W}_t,\\ & \mathrm{d}{Y}_t=\mathrm{d}{B}_t,\end{array}$$

$$u_{\varepsilon}(x,y)=\mathbb{E}\biggl{[}\phi_{Y_{\tau}}(X_{\tau})+\frac{1}{2}\int_{0}^{\tau}f(X_{t},Y_{t})\,\mathrm{d}t~{}\bigg{|}~{}(X_{0},Y_{0})=(x,y)\biggr{]},$$

$$\bar{u}(y):={\int }_0^1{u}_{\epsilon }(x,y)\mathrm{d}x,$$

$${\tilde{u}}_{\epsilon }:=u_{\epsilon }-\bar{u}.$$

$$\{\begin{array}{ll}& -{\bar{u}}^{\prime \prime }(y)=\bar{f}(y),\text{in }(0,1),\\ & \bar{u}(0)={\bar{\varphi }}_0,\bar{u}(1)={\bar{\varphi }}_1,\end{array}$$

$$\bar{f}(y)={\int }_0^{1}f(x,y)\mathrm{d}x,{\bar{\varphi }}_0={\int }_0^1{\varphi }_0(x)\mathrm{d}x,{\bar{\varphi }}_1={\int }_0^1{\varphi }_1(x)\mathrm{d}x.$$

$$\bar{u}(y)=y\biggl{(}\int_{0}^{1}\int_{0}^{z}\bar{f}(t)\,\mathrm{d}t\mathrm{d}z-\bar{\phi}_{0}+\bar{\phi}_{1}\biggr{)}-\int_{0}^{y}\int_{0}^{z}\bar{f}(t)\,\mathrm{d}t\mathrm{d}z+\bar{\phi}_{0}.$$

$$\{\begin{array}{ll}-{\epsilon }^{-2}{\partial }_x^2{\tilde{u}}_{\epsilon }-{\partial }_y^2{\tilde{u}}_{\epsilon }=\tilde{f},& \text{in }D,\\ {\partial }_x{\tilde{u}}_{\epsilon }(0,y)={\partial }_x{\tilde{u}}_{\epsilon }(1,y)=0,& 0⩽y⩽1,\\ {\tilde{u}}_{\epsilon }(x,0)={\tilde{\varphi }}_0(x),{\tilde{u}}_{\epsilon }(x,1)={\tilde{\varphi }}_1(x),& 0⩽x⩽1,\end{array}$$

$$\tilde{f}(x,y)=f(x,y)-\bar{f}(y),{\tilde{\varphi }}_0(x)={\varphi }_0(x)-{\bar{\varphi }}_0,{\tilde{\varphi }}_1(x)={\varphi }_1(x)-{\bar{\varphi }}_1.$$

$${\int }_0^1{\tilde{u}}_{\epsilon }(x,y)\mathrm{d}x=0,\text{for }0⩽y⩽1,$$

$${\int }_0^1\tilde{f}(x,y)\mathrm{d}x=0,\text{for }0⩽y⩽1,$$

$${\int }_0^1{\tilde{\varphi }}_0(x)\mathrm{d}x={\int }_0^1{\tilde{\varphi }}_1(x)\mathrm{d}x=0.$$

$$\{\begin{array}{ll}{\partial }_x^2{\tilde{u}}_0=0,& \text{in }D,\\ {\partial }_x{\tilde{u}}_0(0,y)={\partial }_x{\tilde{u}}_0(1,y)=0,& 0⩽y⩽1,\\ {\int }_0^1{\tilde{u}}_0(x,y)\mathrm{d}x=0,& 0⩽y⩽1,\\ {\tilde{u}}_0(x,0)={\tilde{\varphi }}_0(x),{\tilde{u}}_0(x,1)={\tilde{\varphi }}_1(x),& 0⩽x⩽1.\end{array}$$

$${\tilde{u}}_{\epsilon }^{\mathbf{ot}}(x,y)=\sum _{n=0}^{\infty }{\epsilon }^n{\tilde{u}}_n^{\mathbf{ot}}(x,y).$$

$$-{\partial }_x^2{\tilde{u}}_0^{\mathbf{ot}}=0,$$

$$-{\partial }_x^2{\tilde{u}}_1^{\mathbf{ot}}=0,$$

$$-{\partial }_x^2{\tilde{u}}_2^{\mathbf{ot}}={\partial }_y^2{\tilde{u}}_0^{\mathbf{ot}}+\tilde{f},$$

$$-{\partial }_x^2{\tilde{u}}_n^{\mathbf{ot}}={\partial }_y^2{\tilde{u}}_{n-2}^{\mathbf{ot}},n⩾3.$$

$$\{\begin{array}{ll}{\partial }_x{\tilde{u}}_n^{\mathbf{ot}}(0,y)={\partial }_x{\tilde{u}}_n^{\mathbf{ot}}(1,y)=0,& 0⩽y⩽1,\\ {\int }_0^1{\tilde{u}}_n^{\mathbf{ot}}(x,y)\mathrm{d}x=0,& 0⩽y⩽1.\end{array}$$

$${\tilde{u}}_0^{\mathbf{ot}}(x,y)$$

$${\tilde{u}}_n^{\mathbf{ot}}(x,y)$$

$${\tilde{u}}_2^{\mathbf{ot}}(x,y)$$

$${\tilde{F}}_n(x,y)=\{\begin{array}{ll}{\int }_0^x{\tilde{F}}_{n-1}(z,y)\mathrm{d}z,& n⩾1,\\ \tilde{f}(x,y),& n=0.\end{array}$$

$${\tilde{u}}_{\epsilon }^{\mathbf{in},0}(x,\xi )=\sum _{n=0}^{\infty }{\epsilon }^n{\tilde{u}}_n^{\mathbf{in},0}(x,\xi ).$$

$$-{\epsilon }^{-2}{\partial }_x^2{\tilde{u}}_{\epsilon }(x,\xi )-{\epsilon }^{-2}{\partial }_{\xi }^2{\tilde{u}}_{\epsilon }(x,\xi )=\tilde{f}(x,\epsilon \xi ).$$

$$\{\begin{array}{ll}& {\partial }_x{\tilde{u}}_{\epsilon }^{\mathbf{in},0}(x,\xi )=0,x=0\text{ or }1,\\ & {\tilde{u}}_{\epsilon }^{\mathbf{in},0}(x,0)={\tilde{\varphi }}_0(x),\end{array}$$

$${\int }_0^1{\tilde{u}}_{\epsilon }^{\mathbf{in},0}(x,\xi )\mathrm{d}x=0,\text{for }\xi ⩾0.$$

$$\tilde{f}(x,\epsilon \xi )=\sum _{n=0}^{\infty }\frac{{\epsilon }^n{\xi }^n}{n!}{\partial }_y^n\tilde{f}(x,0),$$

$$\{\begin{array}{ll}& -{\partial }_x^2{\tilde{u}}_0^{\mathbf{in},0}-{\partial }_{\xi }^2{\tilde{u}}_0^{\mathbf{in},0}=0,0<x<1,\xi >0,\\ & {\partial }_x{\tilde{u}}_0^{\mathbf{in},0}(x,\xi )=0,x=0\text{ or }1,\\ & {\tilde{u}}_0^{\mathbf{in},0}(x,0)={\tilde{\varphi }}_0(x),\\ & {\int }_0^1{\tilde{u}}_0^{\mathbf{in},0}(x,\xi )\mathrm{d}x=0,\text{for }\xi ⩾0,\end{array}$$

$$\{\begin{array}{ll}& -{\partial }_x^2{\tilde{u}}_1^{\mathbf{in},0}-{\partial }_{\xi }^2{\tilde{u}}_1^{\mathbf{in},0}=0,0<x<1,\xi >0,\\ & {\partial }_x{\tilde{u}}_1^{\mathbf{in},0}(x,\xi )=0,x=0\text{ or }1,\\ & {\tilde{u}}_1^{\mathbf{in},0}(x,0)=0,\\ & {\int }_0^1{\tilde{u}}_1^{\mathbf{in},0}(x,\xi )\mathrm{d}x=0,\text{for }\xi ⩾0,\end{array}$$

Full text

Asymptotic Expansion with Boundary Layer Analysis for Strongly Anisotropic Elliptic Equations

Ling Lin 111email: linglin@cityu.edu.hk. Corresponding author. and Xiang Zhou 222email: xiang.zhou@cityu.edu.hk. The research of XZ was supported by the grants from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. CityU 11304314, 11304715 and 11337216).

Department of Mathematics

City University of Hong Kong

Tat Chee Ave, Kowloon

Hong Kong SAR

abstract

In this article, we derive the asymptotic expansion, up to an arbitrary order in theory, for the solution of a two-dimensional elliptic equation with strongly anisotropic diffusion coefficients along different directions, subject to the Neumann boundary condition and the Dirichlet boundary condition on specific parts of the domain boundary, respectively. The ill-posedness arising from the Neumann boundary condition in the strongly anisotropic diffusion limit is handled by the decomposition of the solution into a mean part and a fluctuation part. The boundary layer analysis due to the Dirichlet boundary condition is conducted for each order in the expansion for the fluctuation part. Our results suggest that the leading order is the combination of the mean part and the composite approximation of the fluctuation part for the general Dirichlet boundary condition.

**Keywords: **strongly anisotropic elliptic equation, boundary layer analysis, matched asymptotic analysis.

1. Introduction

The strongly anisotropic elliptic problem we consider in this article is the following equation imposed in the domain ${D}=(0,1)\times(0,1)$ with mixed Dirichlet–Neumann boundary conditions:

[TABLE]

The special feature for this equation is that $\varepsilon$ is a small positive number, that is, the diffusion coefficient along the $x$-direction is very large. Here the Neumann boundary conditions are imposed on the left and right boundaries and the Dirichlet boundary conditions are imposed on the top and bottom boundaries of the rectangular domain.

The equation (1) belongs to a large class of diffusion models with the strongly anisotropic diffusion coefficients from many applications, e.g., image processing [13], flows in porous media [1, 8], semiconductor modeling [9], heat conduction in fusion plasmas [12], and so on. Note that here we use $\varepsilon^{2}$ rather than $\varepsilon$ in the diffusion coefficients for the reason we will mention later. Furthermore, in this simplified model (1), the line field parallels to the $x$-axis. In more realistic models, the line field may not be so simple and could be a closed loop. Our main interest is to examine the asymptotic behaviors of the solution $u_{\varepsilon}$ such as $u_{\varepsilon}\sim u_{0}+\varepsilon u_{1}+\varepsilon^{2}u_{2}+\cdots$. To illustrate the main ideas, the equation (1) serves a good model which can simply lots of technical calculations.

We first briefly review the existing works on the limit of the solution as $\varepsilon\downarrow 0$. Mathematically, the Neumann boundary conditions yield an ill-posed limiting problem as taking the formal limit $\varepsilon\downarrow 0$ in the original problem (1) [2]. The traditional numerical methods for the elliptic equations, such as the standard five-point scheme, suffer from the large condition numbers for tiny values of $\varepsilon$. There have been a lot of efforts focusing on the numerical methods for the strongly anisotropic elliptic problems. In particular, the class of asymptotic preserving method was developed recently by P. Degond et al. in a series of papers, e.g., [4, 3, 5, 2]. Their main idea is to decompose the solution into two parts, a mean part along the strongly diffusive direction and a fluctuation part, and then they reformulate the original equation into a coupled system of the equations for these two parts. More recently, [11] proposed a new innovative approach to replace one of the Neumann boundary condition by the integration of the original equation along the field line. Consequently, the singular terms can be replaced by some regular terms, which yields a well-posed limiting problem.

Heuristically, the vanishing $\varepsilon$ means that the variation of the solution $u_{\varepsilon}(x,y)$ along the $x$-direction is very slow, which then suggests that the limiting solution $u_{0}=\lim_{\varepsilon\downarrow 0}u_{\varepsilon}$ (defined in certain sense) is constant along the $x$-direction, i.e., a function of the $y$ variable only. However, this would be inconsistent with the Dirichlet boundary conditions in (1) if either of $\phi_{0}$ and $\phi_{1}$ were nonconstant. In other words, a function of only the $y$ variable (such as the so-called mean part) is, in general, not capable of describing the limiting solution in the whole domain because there exist boundary layers near each nonconstant Dirichlet boundary. Inside these boundary layers, the function of only the $y$ variable has to be corrected to match the nonconstant Dirichlet boundary conditions. It is noteworthy that for existing numerical examples presented in previous works such as [5, 11], as far as the authors know, the Dirichlet boundary conditions are always homogeneous, and so there are no boundary layers. Actually, those numerical examples are intentionally constructed by choosing a true solution without boundary layers first and then defining the force term $f$ accordingly. However, as we argued in the above, the emergence of the boundary layer is generic. This phenomena certainly makes our asymptotic analysis more complicated.

The existence of the boundary layers for the strongly anisotropic elliptic problem (1) can also be easily seen from its probabilistic interpretation which is connected to a special type of the stochastic fast-slow dynamics. The random walk model corresponding to the elliptic problem (1) is very simple. Let $(X_{t},Y_{t})$ be the position of a particle in ${D}$ satisfying the following stochastic differential equation

[TABLE]

subject to the reflection boundary condition on the left and right boundaries ($x=0,1$) and the absorbing boundary condition on the top and bottom boundaries ($y=0,1$). Here $W_{t}$ and $B_{t}$ are two independent (standard) Brownian motions. By the Feymann-Kac formula [10], the solution to (1) is represented by

[TABLE]

where $\tau=\inf\left\{t>0:Y_{t}=0\text{ or }1\right\}$ is the absorption time of the $Y$ process to the Dirichlet boundaries. The solution to (2) is straightforward: $X_{t}=X_{0}+\varepsilon^{-1}W_{t}\overset{\mathrm{d}}{=}X_{0}+W_{t/\varepsilon^{2}}$ and $Y_{t}=Y_{0}+B_{t}$. (“$\overset{\mathrm{d}}{=}$” means the equality in the sense of distribution.) So $X_{t}$ is a fast process and $Y_{t}$ is a slow process. The particle randomly moves drastically fast with the speed at the order $\mathcal{O}(\varepsilon^{-1})$ along the $x$-direction while at a normal speed at the order $\mathcal{O}(1)$ in the $y$-direction. By the averaging principle, the leading order dynamics as the limit of $\varepsilon\downarrow 0$ is the expectation of the slow dynamics for $Y_{t}$ with respect to the invariant measure of the fast variable $X_{t}$, which is a uniform distribution here. Thus the expectation of the integral part in (3) should have a limit independent of the $x$ variable, which is exactly the so called mean part in [4]. But for the expectation of the first term in (3), it depends on the distribution of the absorbing point $(X_{\tau},Y_{\tau})$. If the starting position $(x,y)$ is away from the absorbing boundary, then the absorbing time $\tau$ is sufficient large compared to the $\mathcal{O}(\varepsilon)$ relaxation time to the equilibrium in the $x$-direction, so that the averaging principle still holds, and thus the limit of (3) is a function of the variable $y$ only. However, the averaging principle breaks down if the initial position $(x,y)$ is too close to the absorbing boundary so that $\tau$ would be too short to allow the fast dynamics to relax to the equilibrium. It is easy to see that this occurs if the distance to the boundary is $\mathcal{O}(\varepsilon)$, thus the thickness of the boundary layers around $y=0,1$ is $\mathcal{O}(\varepsilon)$.

Our main motivation is to give a more detailed understanding of the above probabilistic picture by the tool of asymptotic analysis. The goal is to derive a series of approximate functions to the solution $u_{\varepsilon}$ up to an arbitrary order as $\varepsilon\downarrow 0$. In this note, we shall consider the general Dirichlet boundary conditions in (1). This means that we should include the boundary layer analysis in our asymptotic expansion. We shall show that each asymptotic term $u_{0},u_{1},u_{2},\ldots$ exhibits the boundary layer effect. In particular, the leading order $u_{0}(x,y)$ is not simply the mean part $\bar{u}(y)=\int_{0}^{1}u_{\varepsilon}(x,y)\,\mathrm{d}x$. To attack the ill-posedness arising from the Neumann boundary conditions, we utilize the strategy of decomposing the solution into a mean part and a fluctuation part [4, 3, 5, 2]; to deal with the boundary layers originating from the nonconstant Dirichlet boundary conditions, we adopt the Van Dyke’s method of matched asymptotic expansions [6]. Thus the outer expansion and the inner expansion are both conducted. The series in the outer expansion are described by the $y$-parametrized one-dimensional Neumann boundary value problem in the $x$ variable, while the series in the inner expansion are in the form of the two-dimensional elliptic equations which are solved with the aid of the Fourier series.

The rest of the paper is organized as follows. Section 2 presents our main result of the asymptotic expansion. Section 3 gives a rigorous proof of our formal expansion. Section 4 shows the numerical results to validate the convergence order in $\varepsilon$ and demonstrate the computational efficiency. The last section contains our concluding discussion.

2. Asymptotic Result

2.1. Decomposing the solution into the mean value and the fluctuation

For the solution $u_{\varepsilon}$ to (1), we introduce the mean part $\bar{u}$ along the line field, i.e., the $x$-coordinate

[TABLE]

and denote the residual as the fluctuation part $\tilde{u}_{\varepsilon}$,

[TABLE]

Then by integrating both sides of the equation (1) with respect to $x$ over $[0,~{}1]$, we obtain

[TABLE]

where

[TABLE]

Clearly, (4) is a well-posed linear two-point boundary value problem, and $\bar{u}$ can be solved uniquely. Formally,

[TABLE]

Subtracting (4) from (1) yields the PDE for the fluctuating part:

[TABLE]

where

[TABLE]

Note that by construction, we have

[TABLE]

[TABLE]

and

[TABLE]

2.2. Asymptotic expansions of the fluctuation $\tilde{u}_{\varepsilon}$

Our main task is to seek an asymptotic expansion of the fluctuation $\tilde{u}_{\varepsilon}$. Formally, as $\varepsilon\downarrow 0$ in (6) and (7), the formal limit $\tilde{u}_{0}=\lim_{\varepsilon\downarrow 0}\tilde{u}_{\varepsilon}$ would satisfy

[TABLE]

Clearly, this is an ill-posed problem unless $\tilde{\phi}_{0}(x)\equiv 0$ and $\tilde{\phi}_{1}(x)\equiv 0$, since the first three equations in (8) yield $\tilde{u}_{0}\equiv 0$. This inconsistency implies that we have a singular perturbation problem and anticipate the emergence of two boundary layer regions near the Dirichlet boundaries $y=0$ and $y=1$ respectively. We apply the Van Dyke’s method of matched asymptotic expansions [6] to tackle this problem, i.e., first separately solve the problem in the inner regions within the boundary layers and in the outer region away from the boundary layers, and then match them at the edges of the boundary layers.

2.2.1. Outer expansion

Assume the following outer expansion away from the Dirichlet boundaries $y=0$ and $y=1$:

[TABLE]

Substituting this into the equation in (6) and equating coefficients, we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

These equations are a set of parametric one dimensional differential equations in the $x$ variable and $y$ is in the role of parameters. The outer expansion solutions must also satisfy the Neumann boundary condition as in (6) and the integral condition (7), which gives for any $n$

[TABLE]

Then each $\tilde{u}^{\mathrm{\bf{ot}}}_{n}$ is the unique solution to these Neumann problems due to the second condition in (13). We can solve $\tilde{u}^{\mathrm{\bf{ot}}}_{n}$ recursively from (9)$\sim$(12) together with (13). In particular, we have

[TABLE]

Here $\tilde{F}_{n}$ is defined recursively as

[TABLE]

2.2.2. Inner expansion near $y=0$

Next we explore the inner solution near $y=0$ in terms of the stretched variable $\xi=y/{\varepsilon}$ by assuming

[TABLE]

In terms of $\xi$, the equation in (6) becomes

[TABLE]

Thus the inner expansion $\tilde{u}_{\varepsilon}^{\mathrm{\bf{in}},0}(x,\xi)$ near $y=0$ asymptotically satisfies the equation (17) with the boundary conditions

[TABLE]

and the integral condition

[TABLE]

We can write the Taylor expansion of the fluctuation part of the external force:

[TABLE]

then equate coefficients of the same powers of $\varepsilon$ to obtain the following two-dimensional elliptic equations on the domain $(x,\xi)\in(0,1)\times(0,\infty)$:

[TABLE]

[TABLE]

and for $n\geqslant 2$,

[TABLE]

These problems (18),(19) and (20) do not have the uniqueness of the solutions even the integral conditions $\int_{0}^{1}\tilde{u}_{n}^{\mathrm{\bf{in}},0}(x,\xi)\,\mathrm{d}x=0$ are imposed. The uniqueness comes from the matching to the outer solutions as we will show below.

We next solve $\tilde{u}^{\mathrm{\bf{in}},0}_{n}$ by the method of separation of variables because of the simple geometry of the domain. We first work on the lowest order at $n=0$. We look at the solutions that can be expanded into the form $\tilde{u}^{\mathrm{\bf{in}},0}_{0}=\sum_{k}A_{k}(\xi)B_{k}(x)$. The equation and boundary conditions in (18) show that

[TABLE]

which form a complete orthogonal basis for the space $L^{2}([0,~{}1])$. Thus we can expand $\tilde{u}^{\mathrm{\bf{in}},0}_{0}(\cdot,\xi)$ and $\tilde{\phi}_{0}$ respectively in terms of these Fourier cosine series:

[TABLE]

[TABLE]

where the coefficients for $k\geqslant 1$ are

[TABLE]

[TABLE]

Note that the terms for $k=0$ in these Fourier cosine series disappear since

[TABLE]

[TABLE]

Substituting these Fourier cosine expansions into (20), we deduce that for each $k\geqslant 1$, $A_{k}(\xi)$ satisfies

[TABLE]

Hence we have that for $k\geqslant 1$,

[TABLE]

with the constants $c_{k}$ to be determined later by matching the outer and inner solutions.

Remark 2.1**.**

We comment a bit on the possibility of generalizing the above calculations to a general line field. One can work in the curvilinear coordinate of the field line and obtain the equations for the outer expansions straightforwardly. But since in general the line field may not match the Dirichlet boundary like in our model (1), then the boundary layer may not be a rectangular band with a uniform width $\varepsilon$, and consequently, the stretching variable $\xi$ would not be simply equal to $y/\varepsilon$; the geometric property of the field line and the Dirichlet boundary should be incorporated to derive the equations for the inner expansion near the Dirichlet boundary.

2.2.3. Matching

To determine the constants $c_{k}$’s in the first-term approximation of the boundary layer solution near $y=0$, we make use of the essential point that the inner solution $\tilde{u}^{\mathrm{\bf{in}},0}_{0}(x,\xi)$ and the outer solution $\tilde{u}^{\mathrm{\bf{ot}}}_{0}(x,y)$ should match on the boundary of the layer near $y=0$, that is,

[TABLE]

This gives

[TABLE]

and so

[TABLE]

Consequently,

[TABLE]

2.2.4. Inner expansion near $y=1$

For the other Dirichlet boundary at $y=1$, we proceed in the exactly same way to derive the inner solution. Assume the inner expansion near $y=1$ in terms of the stretched variable $\eta=(1-y)/\varepsilon$,

[TABLE]

In terms of $\eta$, the equation in (6) becomes

[TABLE]

Thus the inner expansion $\tilde{u}_{\varepsilon}^{\mathrm{\bf{in}},1}(x,\eta)$ near $y=1$ must asymptotically satisfy this equation and the boundary and integral conditions

[TABLE]

Again, using Taylor expansion and then equating coefficients of like powers, we obtain

[TABLE]

[TABLE]

and for $n\geqslant 2$,

[TABLE]

By the same token, we solve the above equations by Fourier cosine series and use the matching procedure to determine the constants, then we obtain the lowest order

[TABLE]

where

[TABLE]

2.2.5. Composite Expansion

Now we can get the leading order term of $u_{\varepsilon}$ which is valid on the whole domain. Expressing all the three pieces of expansions in terms of $x$ and $y$, and combining them by adding them together and then subtracting their common parts, eventually we obtain the following composite approximation by noting (5), (14), (21) and (22),

[TABLE]

where $\bar{u}(y)$ is the mean solution given in (5).

2.2.6. Higher order approximations

Higher order approximations can be obtained similarly by the Van Dyke’s method of matched asymptotic expansions [6]. Let us compute the second order expansion for demonstration for illustration. To this end, we need to solve the next two orders in the inner expansion near the boundaries $y=0,1$, i.e., $\tilde{u}^{\mathrm{\bf{in}},0}_{n}$ and $\tilde{u}^{\mathrm{\bf{in}},1}_{n}$ for $n=1,2$.

Using the method of separation of variables again, by expanding in the Fourier cosine series, we obtain

[TABLE]

[TABLE]

where

[TABLE]

$a_{k}$ and $b_{k}$ are undetermined constants. Clearly, the $\mathrm{e}^{k\pi\xi}$ terms should disappear since they asymptotically blow up, this implies that $a_{k}=0$ and $b_{k}=0$. Hence, by guessing $a_{k}\equiv 0$ and $b_{k}\equiv 0$, we have the first three terms of the inner solution $\tilde{u}_{\varepsilon}^{\mathrm{\bf{in}},0}(x,\xi)$:

[TABLE]

Note that the outer expansion is

[TABLE]

To see that these two expansions do match up to the given order, we write the outer expansion in terms of the inner variables and vice versa. By dropping the asymptotically negligible $\mathrm{e}^{-k\pi\xi}$ terms as $\xi\to\infty$, we have the $\mathcal{O}(\varepsilon^{3})$ approximations:

[TABLE]

and

[TABLE]

To show that matching (to this order) has been accomplished, we only need to check that the right hand sides of (25) and (24) are equal. In fact, we have for every $y$, by (13) and (15),

[TABLE]

from (16) and the integration by parts twice, for $k\geqslant 1$,

[TABLE]

Note that in the second equality, we also used the simple fact of the integral condition

[TABLE]

Analogously, we can solve

[TABLE]

[TABLE]

and from

[TABLE]

[TABLE]

we also see that matching (to this order) has been accomplished.

The last step is to combine the three expansions into a composite expansion

[TABLE]

where $u^{[0]}(x,y)$ is the leading order given in (2.2.5).

Clearly, using the above method, we may proceed to derive the asymptotic expansions of $u_{\varepsilon}$ to any order, and in general, the form of the asymptotic expansion up to $\mathcal{O}(\varepsilon^{2n})$, $n=0,1,\dots$, is

[TABLE]

The justification of the approximation orders of these asymptotic expansions will be given in the next section.

2.3. Discussion

For our approximations (26), it is observed that the boundary layer terms would disappear if

[TABLE]

i.e.,

[TABLE]

all happen to be constant functions independent of $x$. Note that this condition is stronger than the homogeneous Dirichlet boundary conditions since it also involves the even order normal derivatives up to $2n-2$ of the external force on the Dirichlet boundaries. In this case free of the boundary layers, we only need to compute the outer expansions by solving the equations (11) and (12) up to $2n$ to get the approximation $u^{[2n]}$. Note that each equation in (11) and (12) is actually a system parametrized by the $y$ variable of independent ordinary differential equations in the $x$ variable rather than a two-dimensional partial differential equation. In computation, each equation can be solved by a numerical integrator in parallel at all grid points of the $y$ variable which is now viewed as a parameter. This can reduce the computational cost to linear scaling. However, in general, the appearance of the boundary layers seems inevitable; then many existing algorithms may need further improvements to approximate the solution on the whole domain.

3. Theoretical Justification

Define the remainders in the asymptotic expansions of $u_{\varepsilon}$:

[TABLE]

where the approximations $u^{[2n]}$’s are given in (26). We shall prove the following estimates of the errors $r_{2n}$ in this section.

Theorem 3.1**.**

[TABLE]

To justify these estimates, we first establish a modified version of the maximum principle for the elliptic equation with mixed boundary value conditions.

Lemma 3.2**.**

Let

[TABLE]

be a uniformly elliptic operator on a connected, bounded, open domain $\Omega$. Suppose that $w\in C^{2}(\Omega)\cap C(\bar{\Omega})$ satisfies $\mathcal{L}w\leqslant 0$ in $\Omega$ and $\partial_{\boldsymbol{n}}w\leqslant 0$ on $\Gamma_{N}\subset\partial\Omega$, where $\boldsymbol{n}$ is the outer unit normal to $\Omega$ on the boundary $\partial\Omega$. Also assume that $\Omega$ satisfies the interior ball condition at every $\boldsymbol{x}\in\Gamma_{N}$. Let $\Gamma_{D}=\partial\Omega\setminus\Gamma_{N}$. Then

[TABLE]

Proof.

It suffices to show that $w$ attains its maximum over $\bar{\Omega}$ on $\Gamma_{D}$. Let us assume $w$ is not constant within $\Omega$, as otherwise the proof is trivial. Then the strong maximum principle [7] states that $w$ cannot attain its maximum over $\bar{\Omega}$ at any interior point. Furthermore, $w$ cannot attain its maximum over $\bar{\Omega}$ at any boundary point $\boldsymbol{x}_{0}\in\Gamma_{N}$ either, since otherwise Hopf’s lemma [7] would imply $\partial_{\boldsymbol{n}}w(\boldsymbol{x}_{0})>0$, which contradicts the assumption $\partial_{\boldsymbol{n}}w\leqslant 0$ on $\Gamma_{N}$. Hence $w$ has to attain its maximum over $\bar{\Omega}$ on $\Gamma_{D}$. ∎

The proof of Theorem 3.1 relies on the following lemma, which is a consequence of the above modified version of the maximum principle.

Lemma 3.3**.**

For any $\varepsilon>0$, let $\mathcal{L}=-\partial_{x}^{2}-\varepsilon^{2}\partial_{y}^{2}$. Suppose that $u\in C^{2}({D})\cap C(\bar{{D}})$ satisfies

[TABLE]

Then

[TABLE]

where

[TABLE]

Proof.

Let $w(x,y)=u(x,y)-v(x)$, where

[TABLE]

Then it is easy to check that

[TABLE]

From the modified version of the maximum principle (Lemma 3.2), we conclude that $w(x,y)\leqslant 0$ in ${D}$, thus

[TABLE]

Applying the above argument to $-u$ yields

[TABLE]

Then we obtain the desired inequality. ∎

Now we give the proof of Theorem 3.1.

Proof of Theorem 3.1.

Direct calculation shows that $r_{2n}$ satisfies

[TABLE]

where $\mathcal{L}$ is the same as in Lemma 3.3. Then the asserted error estimates follow from Lemma 3.3. ∎

4. Numerical Example

We choose the source term

[TABLE]

and the Dirichlet boundary conditions

[TABLE]

Note that the Dirichlet boundary conditions should be consistent with the Neumann boundary conditions in (1), which dictates the compatibility conditions

[TABLE]

Clearly, these compatibility conditions are satisfied in this example.

We choose a grid for which the meshpoints are half-integered in the $x$-direction and integered in the $y$-direction, that is,

[TABLE]

where $\Delta x=1/N$, $\Delta y=1/M$, and $i=1,2,\cdots,N$, $j=1,2,\cdots,M+1$. Then the solution $u_{\varepsilon}$ to (1) is solved numerically by the standard five-point finite difference method. Figure 1 shows respectively the contour plots of the numerical solutions with a very fine mesh size $M\times N=512\times 2048$ to the exact solution $u_{\varepsilon}$, the mean part $\bar{u}$, the first order asymptotic approximation $u^{[0]}$, and the higher order asymptotic approximation $u^{[2]}$ for $\varepsilon^{2}=0.05$. It is observed that there exist two boundary layers near the boundaries $y=0$ and $y=1$ respectively, and the thickness of each boundary layer is roughly $\varepsilon=\sqrt{0.05}\approx 0.22$. This is consistent with the result of the asymptotic analysis in last section. Clearly, as shown in the subfigure (b), the mean solution $\bar{u}$ fails to capture the leading order solution inside these two boundary layers. The first order approximation $u^{[0]}$ is very close to the true solution while the next order approximation $u^{[2]}$ is almost identical to the true solution, as indicated from the subfigures (c) and (d).

To further validate the approximation order, we compute the errors in $L_{\infty}$ norm. Table 1 shows the discrete $L_{\infty}$ norms of the errors $r_{0}$ and $r_{2}$ for different values of $\varepsilon$, and Figure 2 is the convergence plot of these errors versus $\varepsilon^{2}$ on logarithmic scales. The best fitting straight lines through each set of data points are also plotted in Figure 2. From the equations of these fitting lines, we observe that the orders of the approximation errors $r_{0}$ and $r_{2}$ are roughly $\mathcal{O}(\varepsilon^{2})$ and $\mathcal{O}(\varepsilon^{4})$ respectively. These numerical results confirm the theoretical assertion in Theorem 3.1.

5. Conclusion

We have presented a formal expansion of the solution to the strongly anisotropic diffusion equation (1) in the simple rectangular domain. Our expansions take into account the boundary layers near the Dirichlet boundaries. We proved the rigorous convergence of the formal expansion by the maximum principle. In theory, our method can be generalized to complicated cases where the field line is not simply the $x$-axis direction, but in the form of curves. For such cases, the analysis of the boundary layers is more difficult. Nevertheless, we expect that our analysis here would provide some ideas of constructing new efficient numerical methods which incorporate both the boundary layer effect due to the Dirichlet boundary condition and the ill-posedness due to the Neumann or periodic boundary condition in the strongly anisotropic diffusion limit.