A Mixed Finite Element Method for Singularly Perturbed Fourth Oder Convection - Reaction - Diffusion Problems on Shishkin Mesh

Charuka D. Wickramasinghe Department of Mathematics, University of Kentucky, Lexington, KY 40508, USA [email protected]

Abstract.

This paper introduces an approach to decoupling singularly perturbed boundary value problems for fourth-order ordinary differential equations that feature a small positive parameter $\epsilon$ multiplying the highest derivative. We specifically examine Lidstone boundary conditions and demonstrate how to break down fourth-order differential equations into a system of second-order problems, with one lacking the parameter and the other featuring $\epsilon$ multiplying the highest derivative. To solve this system, we propose a mixed finite element algorithm and incorporate the Shishkin mesh scheme to capture the solution near boundary layers. Our solver is both direct and of high accuracy, with computation time that scales linearly with the number of grid points. We present numerical results to validate the theoretical results and the accuracy of our method.

Keywords: Shishkin mesh; finite element algorithm; boundary layers; convection-diffusion problems.

1. Introduction

In this paper, we consider the following stationary state fourth order singularly perturbed differential equation with Lidstone boundary conditions.

\left.\begin{aligned} L_{\epsilon}&=-\varepsilon u^{iv}(x)-a(x)u^{{}^{\prime\prime\prime}}(x)+b(x)u^{{}^{\prime\prime}}(x)=-f(x)\\ u(0)&=0,\quad u(1)=0,\quad u^{{}^{\prime\prime}}(0)=0,\quad u^{{}^{\prime\prime}}(1)=0\end{aligned}\right\}\quad

(1.1)

The function $f(x)$ , $a(x)$ , and $b(x)$ are all smooth and satisfy $a(x)\geq\alpha>0$ and $b(x)\geq\beta>0$ . The parameter $\epsilon$ is assumed to be a small positive value, such that $0<\epsilon\leq 1$ . While it may be easy to analytically solve this problem in some cases, finding the solution $u$ with analytical techniques can be difficult or even impossible for a general function $f$ . The well-posedness of problem (1.1) has been discussed in more detail in sources such as [1] and [2].

The convection-diffusion-reaction equation is used in three processes [28]: convection, which involves the movement of materials from one region to another; diffusion, which involves the movement of materials from an area of high concentration to an area of low concentration; and reaction, which involves decay, adsorption, and the reaction of substances with other components. Singularly perturbed problems have many applications in engineering and applied mathematics, including chemicals and nuclear engineering [5], linearized Navier-Stokes equation at high Reynolds number ([3], [4]), elasticity, aerodynamics, oceanography, meteorology, modeling of semiconductor devices [15], control theory, and oil extraction from underground reservoirs [13], and in many other fields. However, solving these problems numerically presents major computational difficulties due to boundary layers, where the solution changes rapidly. The study of second-order singularly perturbed differential equations are quite large, as seen in the extensive literature ([6], [7],[8],[9],[10],[11],[12],[13],[14]) and their references. However, very few studies have addressed singularly perturbed fourth-order boundary value problems in the literature.

The following paragraph presents a brief overview of analytical and numerical methods used for solving singularly perturbed fourth-order differential equations. In [16], Vrabel, R studied the existence of solutions for fourth-order nonlinear equations with the Lidstone boundary conditions. In [17], Essam R introduced the Differential Transform Method (DTM) as an alternative to existing methods for solving higher-order singularly perturbed boundary value problems (SPBVPs). In [18], Vrabel, R provided a detailed analysis of the boundary layer phenomenon subject to the Lidstone boundary conditions by analyzing the integral equation associated with the SPBVPs. In [19], Sun, G, and Martin S developed Piecewise polynomial Galerkin finite-element methods on a Shishkin mesh, which achieved almost optimal convergence results in various norms. In 1991, Ross and Stynes studied a fourth-order problem and generated an approximate solution using patched basis functions. The method is uniformly first-order convergent in the $H^{1}$ [0,1] norm [20]. In [21], Shanthi, V., and N. Ramanujam transformed the SVBVP into a system of weakly coupled system of two second-order ODEs and then used the fitted operator method (FOM), fitted mesh method (FMM), and boundary value technique (BVT) to approximate the solution.

The problem being studied involves rapidly changing solutions in very thin regions near the boundary. Traditional numerical methods often fail to accurately capture these changes, which can result in errors across the entire domain. To address this issue, various methods such as Bakhavalov and Gartland meshes have been developed ([22], [23], [24]). In this study, we analyze a standard finite element method combined with the Shishkin mesh, which is a type of local refinement strategy introduced by a Russian mathematician Grigorii Ivanovich Shishkin in 1988 [25]. Finite element methods on Shishkin meshes in 1D were first studied 1995 by Sun and Stynes [2], the analysis for second order problems was also published in the two books ([22], [26]) from 1996. The goal is to propose a new mixed finite element algorithm that is reliable, effective, and easy to implement, and can be used to solve higher order singularly perturbed differential equations. As far as the authors are aware, there is currently no finite element algorithm available to solve (1.1) using their decoupled formulation under the Shishkin mesh. The main advantage of decoupling system is that it reduces both memory space and time requirements.

Remark 1.1.

For simplicity, the current paper focuses on analyzing a 1-dimensional problem where the coefficients $a(x)$ and $b(x)$ are both constant. However, the analysis could be extended to 2 dimensions and variable coefficients, although this may present some challenges. Additionally, the problem could be expanded to include non-homogeneous boundary conditions through a simple linear transformation.

The rest of the article is organized as follows: In section 2 we introduce the decouple formulation of (1.1). In section 3 we present Shishkin mesh method and the finite element algorithm. We also present error estimate results of the decouple formulations. In section 4 we present numerical results to validate out theoretical results and conclusion is the section 5. Throughout the following text, the generic positive constants $C$ , $a$ and $b$ may take different values in different formulas but is always independent of the mesh and the small positive parameter $\epsilon$ .

2. The Decouple Formulation

Hereafter, we will consider the following problem by setting $a(x)=a=$ constant and $b(x)=b=$ constant in equation (1.1).

\left.\begin{aligned} L_{\epsilon}&=-\varepsilon u^{iv}(x)-au^{{}^{\prime\prime\prime}}(x)+bu^{{}^{\prime\prime}}(x)=-f(x)\\ u(0)&=0,\quad u(1)=0,\quad u^{{}^{\prime\prime}}(0)=0,\quad u^{{}^{\prime\prime}}(1)=0\end{aligned}\right\}\quad

(2.1)

The function $f(x)$ , $a(x)$ , and $b(x)$ are assumed to be sufficiently smooth for $0\leq x\leq 1$ where,

\left.\begin{aligned} a&\geq\alpha>0\\ b&\geq\beta>0\\ a+\frac{1}{2}b^{{}^{\prime}}&>c>0\hskip 14.22636pt\mbox{for all }x\in[0,1]\end{aligned}\right\}\quad

(2.2)

Under the conditions in (2.2) the problem 2.1 is well posed [2]. Let $(\cdot,\cdot)$ denote the usual $L^{2}(0,1)$ inner product. We define the bilinear form of the equation (2.1) as follows:

A_{\varepsilon}(u,v)=(-\varepsilon u^{{}^{\prime\prime}},v^{{}^{\prime\prime}})+(au^{{}^{\prime\prime}},v^{{}^{\prime}})-(bu^{{}^{\prime}},v^{{}^{\prime}})=-(f,v)\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptu,v\in H^{2}_{0}(0,1)

(2.3)

The weighted energy norm is given by

|||v|||=\{\varepsilon|v|^{2}_{2}+||v||^{2}_{1}\}^{\frac{1}{2}}\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptv\in H^{2}_{0}(0,1)

Lemma 2.1.

Assume (2.2) holds. Then there exist positive constants $C_{1}$ , and $C_{2}$ such that for all $u,v\in H^{2}_{0}(0,1)$ ,

|A_{\varepsilon}(u,v)|\leq C_{1}\varepsilon^{-{\frac{1}{2}}}|||u|||\cdot|||v|||

(2.4)

and

C_{2}|||u|||^{2}\leq|A_{\varepsilon}(u,u)|

(2.5)

Proof.

: Based on Cauchy-Schwarz inequality it is easy to see that

|A_{\varepsilon}(u,v)|\leq\begin{cases}C_{1}|||u|||\cdot||v||_{2}\\ C_{1}||u||_{2}\cdot|||v|||\end{cases}

(2.6)

Since $0<\varepsilon\leq 1$ (2.4) follows immediately. The proof of (2.5) can be shown by following the proof of the lemma (2.1) in [19]. ∎

The weak formulation of (2.1) is to find $u\in H^{2}_{0}(0,1)$ such that

A_{\varepsilon}(u,v)=(-f,v),\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptv\in H^{2}_{0}(0,1)

(2.7)

By the Lax-Milgram lemma, (2.5) has a unique solution $u$ in $H^{2}_{0}(0,1)$ .

We are now able to present our decoupled formulation. It is widely recognized that numerical solutions of higher order problems, such as (2.1), are significantly more challenging than those of lower order problems. To address this issue, we decouple (2.1) into a system of lower order differential equations, as follows:

\left.\begin{aligned} -w^{{}^{\prime\prime}}(x)&=f(x)\hskip 14.22636pt\mbox{for}\hskip 14.22636ptx\in(0,1)\\ w(0)&=0,\quad w(1)=0\end{aligned}\right\}\quad

(2.8)

\left.\begin{aligned} -\epsilon u^{{}^{\prime\prime}}(x)-au^{{}^{\prime}}(x)+bu(x)&=w(x)\\ u(0)&=0,\quad u(1)=0\end{aligned}\right\}\quad

(2.9)

The equation represented by Equation (2.8) is a standard Poison equation, which has the same source term $f(x)$ as Equation (2.1). Assuming that $f(x)$ belongs to $L^{2}(0,1)$ and the given boundary conditions for Equation (2.8) are met, the problem defined by Equation (2.8) is well-posed, according to [27]. These kind of finite element decouple formulations can be found from ([29]), ([30]) for some particular problems in 2D.

Equation (2.9) is a second-order problem that involves a convection-reaction-diffusion process with a singular perturbation. The source term for this equation is represented by $w(x)$ , which is the solution to the problem defined by Equation (2.8). Under following assumptions,

\left.\begin{aligned} a&\geq\beta>0\\ b&\geq 0\\ b+\frac{a^{{}^{\prime}}}{2}&>0\hskip 14.22636pt\mbox{for all }x\in[0,1]\end{aligned}\right\}\quad

(2.10)

the problem defined by Equation (2.9) is also well-posed, as stated in [20].

In order to establish a connection between the solution $u$ obtained from Equations (2.8) and (2.9) and the fourth-order problem represented by Equation (2.1), we present the following lemma. Let us define $H^{m}(0,1)$ as the Sobolev space comprising functions whose $i^{th}$ derivative, $0\leq i\leq m$ , is square-integrable.

Lemma 2.2.

The solution $u\in H^{4}(0,1)$ obtained through (2.8) and (2.9) satisfies the following fourth order differential equation.

\left.\begin{aligned} -\varepsilon u^{iv}(x)-au^{{}^{\prime\prime\prime}}(x)+bu^{{}^{\prime\prime}}(x)&=-f(x)\\ u(0)=0,\quad u(1)=0,\quad u^{{}^{\prime\prime}}(0)=0,\quad u^{{}^{\prime\prime}}(1)&=0\end{aligned}\right\}\quad

(2.11)

Proof.

: We first apply the differential operator $L=d^{2}/dx^{2}$ to both sides of Equation (2.9), which gives us:

-\varepsilon u^{iv}(x)-au^{{}^{\prime\prime\prime}}(x)+bu^{{}^{\prime\prime}}(x)=w^{{}^{\prime\prime}}(x)

(2.12)

As we know from Equation (2.8), $-w^{{}^{\prime\prime}}(x)=f(x)$ . Therefore, the equation:

-\varepsilon u^{iv}(x)-au^{{}^{\prime\prime\prime}}(x)+bu^{{}^{\prime\prime}}(x)=-f(x)

holds.

To verify the boundary conditions, we apply the differential operator $L$ to the boundary conditions $u(0)=0$ and $u(1)=0$ in Equation (2.9). This yields the boundary conditions $u^{{}^{\prime\prime}}(0)=0$ and $u^{{}^{\prime\prime}}(1)=0$ . Thus, the conclusion of Lemma (2.2) holds. ∎

3. The Finite Element Method on Shishkin Mesh

In this section, we present a linear finite element method to solve the singularly perturbed boundary value problem (2.1) based on the results obtained in the previous section. We utilize the Shishkin mesh to present error estimates for the singularly perturbed convection reaction diffusion problem (2.9), and we adopt the same definitions and notations as in [2].

3.1. Layer Adapted Shishkin Mesh

Given an even positive integer N, the Shishkin mesh $X_{s}^{N}$ is constructed by dividing the interval [0,1] into two subintervals. We shall consider a mesh $X_{s}^{N}:0=x_{0}<x_{1}<x_{2}<\cdot\cdot\cdot<x_{n-1}<x_{n}=1$ that is equidistant in $[\tau,1]$ but graded in $[0,\tau]$ , where we choose the transition point $\tau$ as Shishkin does:

\tau=min\{1/2,(s+1)\alpha^{-1}\varepsilon N\}

which depends on $\varepsilon$ and N where, s is the order of the highest derivative.

Refer to caption — Figure 1. 1D Shishkin mesh with transition point $\tau$ .

Assumption 3.1.

In this paper we will assume that

\epsilon\leq CN^{-1}

(3.1)

as is generally the case for discretizations of convection-dominated problems.

3.2. The Finite Element Algorithm

Let V be a Hilbert space with norm $\|\cdot\|_{v}$ (but we shall often omit the subscript v to simplify the notation) and scalar product . In the discretization of second-order differential equations with domain $I$ , one generally chooses $V$ as a subset of the Sobolev space $H^{1}(I)$ .

Let $I=[0,L]$ be an interval and let the n+1 node points $\{x_{i}\}_{i=0}^{n}$ define a partition. $I:0=x_{0}<x_{1}<x_{2}<\cdot\cdot\cdot<x_{n-1}<x_{n}=L$ of $I$ into n sub intervals of length $h_{i}=x_{i}-x_{i-1}$ for $i=1,\cdot,\cdot,\cdot,n$ , and $H=max_{i}h_{i}$ . On the mesh $I$ we define the space $V_{n}\subseteq V$ of continuous piecewise linear functions by $V_{n}=\{v:v\in C^{0}(I),v|_{I_{i}}\in P_{I}(I_{i})\}$ where $C^{0}(I)$ denotes the space of continuous functions on $I$ , and $P_{I}(I_{i})$ denotes the space of linear functions on $I_{i}$ .

Lemma 3.2.

Let, $A^{n}_{\varepsilon}(u,v)$ be the discrete bilinear form of $A_{\varepsilon}(u,v)$ in the equation (2.3). There exists a positive constant $h_{0}$ (independent of $\varepsilon$ ) such that for $H\leq h_{0}$ , we have

C_{1}|||v|||^{2}\leq A^{n}_{\varepsilon}(v,v)\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptv\in H^{2}_{0}

(3.2)

Proof.

: The proof can be done by using the argument of Lemma 4.1 of [19], with an additional integration by parts. ∎

Then the Galerkin discretization of problem (2.1) is to find $u_{n}\in V_{n}$ such that

A^{n}_{\varepsilon}(u_{n},v)=(f,v)\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptv\in V_{n}

(3.3)

We will use the following notations:

	$\displaystyle a(u,v)=\int_{0}^{1}u^{{}^{\prime}}(x)v^{{}^{\prime}}(x)\,dx$
	$\displaystyle(f,v)=\int_{0}^{1}f(x)v(x)\,dx$
	$\displaystyle b(u,v)=\varepsilon a(u,v)+(u^{{}^{\prime}},v)+(u,v)\hskip 14.22636pt\mbox{for all}\hskip 14.22636ptu,v\in V:=H_{0}^{1}(0,1)$

The variational formulation of (2.8) is to find $w\in V:=H_{0}^{1}(0,1)$ such that

a(w,v)=(f,v)\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptv\in V

(3.4)

The variational formulation of (2.9) is to find $u\in V:=H_{0}^{1}(0,1)$ such that

b(u,v)=(f,v)\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptv\in V

(3.5)

The following algorithm summarizes the basic steps in computing the finite element solution $w_{n}$ for the two point boundary value problem (2.8).

Algorithm 1

Step 1: Create a mesh

X_{s}^{N}:0=x_{0}<x_{1}<x_{2}<\cdot\cdot\cdot<x_{n-1}<x_{n}=1

and define the corresponding space of continuous piecewise linear functions

V_{n,0}=\{v\in V_{n}:v(0)=v(1)=0\}

with hat basis functions.

\varphi_{j}(x_{i})=\begin{cases}1&\text{if $i=j$ },\\ 0&\text{if $i\neq j$ }\hskip 14.22636pt\mbox{for }\hskip 14.22636pti,j=0,1,2,\cdot\cdot\cdot,n.\end{cases}

Step 2: Compute the

(n-1)\times(n-1)

matrix A and the

(n-1)\times 1

vector b,with entries

A_{i,j}=\int_{0}^{1}\varphi^{{}^{\prime}}_{j}\varphi^{{}^{\prime}}_{i}\,dx\hskip 28.45274ptb_{i}=\int_{0}^{1}f\varphi_{i}\,dx

(3.6)

Step 3: Solve the linear system

A\xi=b

(3.7)

Step 4: Set

w_{n}=\sum_{j=1}^{n-1}\xi_{j}\varphi_{j}

(3.8)

In the same way, we can also define the finite element algorithm for $u_{n}$ in the two-point boundary value problem (2.9). With these algorithms established, we can now introduce our main mixed finite element algorithm for solving the fourth-order singularly perturbed convection reaction diffusion equation (2.11). To obtain the finite element solution of the singularly perturbed boundary value problem (2.11), we utilize the decomposition method described in the previous section.

Algorithm 2

Let

k

be the order of the interpolation polynomial. For, any

f\in H^{-1}(0,1)

and

k\geq 1

, we consider the following steps:

Step 1: Find

w_{n}\in V_{n}^{k}

in the weak formulation of the Poisson equation (2.8) on a Shishkin mesh

X_{s}^{N}

a(w_{n},v)=(f,v)\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptv\in V

(3.9)

Step 2: Set,

-w_{n}=f_{\epsilon}

(3.10)

where

f_{\epsilon}

is the source terms of the equation singularly purtubed second order differential equation (2.9).

Step 3: Then find

u_{n}\in V_{n}^{k}

in the weak formulation of the equation (2.9) on the Shishkin mesh

X_{s}^{N}

for sufficiently large N (an even positive integer) independent of

\epsilon

b(u_{n},v)=(f_{\epsilon},v)\hskip 8.5359pt\mbox{for all}\hskip 8.5359ptv\in V

(3.11)

The discretized linear systems corresponding to the stiffness matrix S, convection matrix C, and mass matrix M are obtained as shown below. Equation (3.12) represents the discretized linear system of (3.9), while equation (3.13) corresponds to the discretized linear system of (3.11).

\left(\frac{\epsilon}{h}S\right)W=F_{1}

(3.12)

\left(\frac{\epsilon}{h}S-C+\frac{h}{6}M\right)U=F_{2}

(3.13)

where $U,F\in\mathbb{R}^{n-1}$ and $S,C,M\in\mathbb{R}^{{n-1}\times{n-1}}$ with:

W:=\begin{bmatrix}w_{1}\\ w_{2}\\ \vdots\\ w_{n-1}\end{bmatrix}\;,\hskip 8.5359ptF_{1}:=\begin{bmatrix}(f,\phi_{1})\\ (f,\phi_{2})\\ \vdots\\ (f,\phi_{n-1})\end{bmatrix}\;,\hskip 8.5359ptU:=\begin{bmatrix}u_{1}\\ u_{2}\\ \vdots\\ u_{n-1}\end{bmatrix}\;,\hskip 8.5359ptF_{2}:=\begin{bmatrix}(f_{\epsilon},\phi_{1})\\ (f_{\epsilon},\phi_{2})\\ \vdots\\ (f_{\epsilon},\phi_{n-1})\end{bmatrix}\;

S:=\begin{bmatrix}2&-1\\ -1&2&-1\\ &\ddots&\ddots&\ddots\\ &&-1&2&-1\\ &&&-1&2\\ \end{bmatrix},\hskip 8.5359ptC:=\frac{1}{2}\begin{bmatrix}0&1\\ -1&0&1\\ &\ddots&\ddots&\ddots\\ &&-1&0&1\\ &&&-1&0\\ \end{bmatrix}

M:=\begin{bmatrix}2&1\\ 1&4&1\\ &\ddots&\ddots&\ddots\\ &&1&4&1\\ &&&1&2\\ \end{bmatrix}

3.3. The Error Estimates

In this section, we present maximum norm error estimate result for our model problem using linear finite elements.

L_{\epsilon}=-\varepsilon u^{iv}(x)-au^{{}^{\prime\prime\prime}}(x)+bu^{{}^{\prime\prime}}(x)=-f(x)

u(0)=0,\quad u(1)=0,\quad u^{{}^{\prime\prime}}(0)=0,\quad u^{{}^{\prime\prime}}(1)=0

As a special case to our model problem we present maximum norm error estimate result with $\varepsilon=1$ for linear finite elements.

Theorem 3.3.

Let u be the solution of problem (2.1). Let, $u_{n}$ be the solution of (3.3) on the Shishkin mesh $X_{s}^{N}$ and $k$ be the order of the interpolation polynomial. Then for N sufficiently large (independently of , $\varepsilon$ ) we have

\|u-u_{n}\|_{\infty}\leq C(N^{-1}lnN)^{min(2,k+1)}

(3.14)

Proof.

: The result follows from the argument of Corollary 5.1 in ([2]). ∎

Remark 3.4.

The proposed decouple formulation is also valid for the special case of $\varepsilon=1$ . In this case there is no boundary layer and thus the Shishkin mesh can be replaced by an uniformly refined mesh (equidistant).

Lemma 3.5.

Suppose that we use a sufficiently accurate quadrature rule, namely, that $k+1\geq 2$ . Let u be the exact solution to the equation (2.1) with $\varepsilon=1$ . and $u_{n}$ be the finite element approximation to the weak formulation of the equation (2.9) with $\varepsilon=1$ . on a uniform mesh. Then we have

\|u-u_{n}\|_{\infty}\leq CN^{-2}

(3.15)

Proof.

: On a uniformly refined mesh it is well known that one has

|||u-u_{n}|||\leq CN^{-2}

(3.16)

It is easy to see that

\|u-u_{n}\|_{\infty}\leq\|u-u_{n}\|_{1}\leq|||u-u_{n}|||

(3.17)

Combining the equations (3.16) and (3.17) the conclusion holds. ∎

4. Numerical Results

In this section we present a few numerical experiments to illustrate the computational method discussed in this paper. The numerical experiments are performed on a laptop computer with MATLAB R2022a in MacBook Air with M1 chip. The linear finite elements are used to solve our model problem. We use the following numerical convergence rate to validate the theoretical convergence rates:

R=\frac{ln\|u-u_{n}\|_{\infty}-ln\|u-u_{n-1}\|_{\infty}}{ln2}

where $u_{n}$ is the finite element approximation after n mesh refinements and $u$ is the exact solution at the same mesh level as the finite element approximation $u_{n}$ is calculated. We use the following model problem to validate theoretical results over the following three examples.

-\epsilon u^{iv}(x)-u^{{}^{\prime\prime\prime}}(x)+u^{{}^{\prime\prime}}(x)=-f(x)\hskip 14.22636ptfor\hskip 14.22636ptx\in(0,1)

(4.1)

u(0)=0,u(1)=0,u^{{}^{\prime\prime}}(0)=0,u^{{}^{\prime\prime}}(1)=0

(4.2)

where $f(x)$ is chosen so that the exact solution is

u(x)=c_{1}e^{r_{1}x}+c_{2}e^{r_{2}x}-\frac{x^{2}+x+1}{2}-\epsilon

(4.3)

where,

c_{1}=-c_{2}+(\frac{1}{2}+\epsilon)

c_{2}=\frac{e^{r_{1}}(\frac{1}{2}+\epsilon)-(\frac{3}{2}+\epsilon)}{e^{r_{1}}-e^{r_{2}}}

r_{1}=\frac{2}{1+\sqrt{(1+4\epsilon)}}

r_{2}=\frac{2}{1-\sqrt{(1+4\epsilon)}}

Example 4.1.

The purpose of this example is to demonstrate that the solution obtained through our decoupled formulation converges to the exact solution. We compare the finite element solution to equation (4.1) with its exact solution (4.3) using both a uniformly refined mesh and a Shishkin mesh. We present the maximum norm error $||u-u_{n}||_{\infty}$ for different values of $\varepsilon$ and different mesh sizes $N$ . Table 1 summarizes the maximum norm error for $\varepsilon=10^{-10},10^{-8},$ and $10^{-6}$ for both uniformly refined meshes and Shishkin meshes. It is well known that singularly perturbed problems have boundary layers, and uniformly refined meshes do not capture the solution near boundary layers well. This phenomenon is evident from the Table 1 values for uniformly refined meshes. To address this issue, a Shishkin mesh is introduced to capture the solution. The values in Table 1 demonstrate that our decoupled formulation solution converges to the true solution under Shishkin meshes.

Table 2 presents the maximum norm error $||u-u_{n}||_{\infty}$ for $\varepsilon=10^{-4},10^{-2},$ and $1$ for both uniformly refined meshes and Shishkin meshes. Although a Shishkin mesh is not required for $\varepsilon=1,$ we included the error results for completeness. By comparing the results in Table 1 and Table 2, we can observe that as we increase the value of $\varepsilon$ from $\varepsilon=10^{-10}$ to $\varepsilon=1,$ the solution under uniformly refined meshes gradually converges to the true solution. This observation is strongly in agreement with existing theoretical results.

Remark 4.2.

All of the $N$ and $\varepsilon$ values in Table 1 satisfy Assumption 1. However, for some values of $N$ and $\varepsilon$ in Table 2, we have $\varepsilon\geq CN^{-1}$ .

Figures 3 and 5 illustrate the numerical solutions and the exact solution under uniformly refined meshes for different values of $N$ and $\varepsilon$ . Figures 4 and 6 show the numerical solutions and the exact solution under corresponding values of $N$ and $\varepsilon$ as in Figures 3 and 5. When the singularly perturbed effect is high, it can be observed from Figures 3 and 5 that the numerical solution does not converge to the true solution under a uniform mesh due to the boundary layer occurring at the point $x=0$ . However, from Figures 4 and 6, it can be seen that our algorithm converges to the true solution under a Shishkin mesh.

Table 1.

\|u-u_{n}\|_{\infty}

for

\varepsilon=10^{-10},10^{-8}\mbox{ and }10^{-6}\mbox{ for uniform and Shishkin meshes}

	$\epsilon=10^{-10}$		$\epsilon=10^{-8}$		$\epsilon=10^{-6}$
$N$	Uniform	Shishkin	Uniform	Shishkin	Uniform	Shishkin
4	0.0354	0.0558	0.0354	0.0559	0.0354	0.0559
8	0.0297	0.0152	0.0297	0.0153	0.0297	0.0153
16	0.0293	0.0039	0.0293	0.0040	0.0293	0.0040
32	0.0296	9.8165e-04	0.0296	0.0010	0.0296	0.0010
64	0.0299	2.3380e-04	0.0299	2.5559e-04	0.0301	2.5535e-04
128	0.0300	5.2171e-05	0.0300	6.4006e-05	0.0311	6.3812e-05
256	0.0301	1.5676e-05	0.0301	1.5991e-05	0.0342	1.5830e-05
512	0.0302	5.3830e-06	0.0303	4.8273e-06	0.0429	4.7611e-06
1024	0.0302	1.8953e-06	0.0309	1.4789e-06	0.0512	1.4328e-06
2048	0.0302	7.0754e-07	0.0329	4.4459e-07	0.0514	4.2862e-07
4096	0.0303	2.8955e-07	0.0395	1.3130e-07	0.0510	1.2852e-07
8192	0.0307	4.0176e-08	0.0500	3.8094e-08	0.0501	3.8150e-08

Table 2.

\|u-u_{n}\|_{\infty}

for

\varepsilon=10^{-4},10^{-2}

and 1 for uniform and Shishkin meshes.

	$\epsilon=10^{-4}$		$\epsilon=10^{-2}$		$\epsilon=1$
$N$	Uniform	Shishkin	Uniform	Shishkin	Uniform	Shishkin
4	0.0355	0.0558	0.0371	0.0429	0.0011	6.9689e-04
8	0.0299	0.0152	0.0366	0.0087	2.2395e-04	1.7362e-04
16	0.0306	0.0039	0.0277	0.0020	4.9304e-05	4.3367e-05
32	0.0350	9.8165e-04	0.0141	0.0010	1.1558e-05	1.0840e-05
64	0.0450	2.3380e-04	0.0046	5.2399e-04	2.7982e-06	2.7117e-06
128	0.0491	5.2171e-05	0.0010	2.8929e-04	6.8864e-07	6.7792e-07
256	0.0467	1.5676e-05	2.4386e-04	1.5611e-04	1.7081e-07	1.6948e-07
512	0.0422	5.3830e-06	6.0942e-05	8.2720e-05	4.2536e-08	4.2370e-08
1024	0.0342	1.8953e-06	1.5157e-05	4.3196e-05	1.0613e-08	1.0593e-08
2048	0.0221	7.0754e-07	3.7826e-06	2.2282e-05	2.6510e-09	2.6486e-09
4096	0.0097	2.8955e-07	9.4508e-07	1.1371e-05	6.6242e-10	6.6407e-10
8192	0.0027	1.3217e-07	2.3621e-07	5.7469e-06	1.6553e-10	1.7331e-10

Example 4.3.

The purpose of this example is to demonstrate the convergent rates and validate the results in Theorem 5 and Lemma 5. Tables 3 and 4 present the convergent rates of the numerical solutions for different $\varepsilon$ values and mesh sizes $N$ . In Table 3, for $\varepsilon=10^{-10},10^{-8},10^{-6}$ , there is no convergent rate under a uniformly refined mesh due to the high singularly perturbed features of the problem. However, as we decrease the singular perturbed properties of the problem, the uniform mesh starts showing convergent rates. For example, in Table 4, when $\varepsilon=1$ , the convergent rate $R=2$ is observed under a uniformly refined mesh, which is in strong agreement with Lemma 5. The convergent rates under a Shishkin mesh for $\varepsilon=10^{-10},10^{-8},10^{-6},10^{-4},10^{-2}$ are also in strong agreement with Theorem 5. This concludes that our proposed algorithm works well on fourth-order singularly perturbed problems like (2.1).

Figure 7 supports the same conclusion as the tables. It shows that there is no convergent rate for uniformly refined meshes with high singularly perturbed properties of the problem. However, with a Shishkin mesh, we can observe the expected convergent rates from Theorem 5.

Table 3. Convergent rate

R

for

\varepsilon=10^{-10},10^{-8}\mbox{ and }10^{-6}\mbox{ for uniform and Shishkin meshes}

	$\epsilon=10^{-10}$		$\epsilon=10^{-8}$		$\epsilon=10^{-6}$
$N$	Uniform	Shishkin	Uniform	Shishkin	Uniform	Shishkin
16	0.0210	1.9404	0.0210	1.9404	0.0204	1.9406
32	0.0156	1.9737	0.0156	1.9737	0.0179	1.9740
64	0.0134	1.9899	0.0135	1.9899	0.0230	1.9909
128	0.0080	1.9975	0.0084	1.9976	0.0457	2.0006
256	0.0043	2.0008	0.0059	2.0009	0.1372	2.0111
512	0.0023	1.7275	0.0086	1.7280	0.3273	1.7333
1024	0.0014	1.7051	0.0260	1.7067	0.2550	1.7325
2048	0.0016	1.7288	0.0919	1.7340	0.0059	1.7410
4096	0.0044	1.7425	0.2631	1.7596	0.0117	1.7377
8192	0.0163	1.7330	0.3392	1.7853	0.0236	1.7522

Table 4. Convergent rate

R

for

\varepsilon=10^{-4},10^{-2}\mbox{ and }1\mbox{ for uniform and Shishkin meshes}

	$\epsilon=10^{-4}$		$\epsilon=10^{-2}$		$\epsilon=1$
$N$	Uniform	Shishkin	Uniform	Shishkin	Uniform	Shishkin
16	0.0297	1.9514	0.4017	2.1152	2.1834	2.0012
32	0.1943	2.0032	0.9771	0.9764	2.0927	2.0003
64	0.3657	2.0699	1.6094	0.9658	2.0464	1.9990
128	0.1240	2.1640	2.1549	0.8570	2.0227	2.0000
256	0.0711	1.7347	2.0872	0.8899	2.0113	2.0000
512	0.1478	1.5421	2.0005	0.9163	2.0057	2.0000
1024	0.3023	1.5060	2.0074	0.9373	2.0028	2.0000
2048	0.6301	1.4215	2.0026	0.9551	2.0012	1.9998
4096	1.1936	1.2890	2.0009	0.9705	2.0007	1.9958
8192	1.8130	1.1314	2.0004	0.9845	2.0006	1.9380

Example 4.4.

In this example, we are comparing the CPU time required by our finite element algorithm to solve problem (2.1) on uniform and Shishkin meshes. The results are presented in Table 5 and Table 6. From the tables, we can observe that the algorithm converges to the true solution faster under the Shishkin mesh compared to the uniform mesh.

Table 5. The

L^{\infty}

error

\|u_{R}^{j}-\phi_{j}\|_{L^{\infty}(\Omega)}

in the L-shaped domain on quasi-uniform meshes.

	$\epsilon=10^{-10}$		$\epsilon=10^{-8}$		$\epsilon=10^{-6}$
$N$	Uniform	Shishkin	Uniform	Shishkin	Uniform	Shishkin
512	0.05	0.09	0.08	0.08	0.05	0.07
1024	0.13	0.13	0.15	0.11	0.11	0.11
2048	0.43	0.26	0.41	0.26	0.40	0.27
4096	1.43	0.76	1.42	0.78	1.24	0.81
8192	5.21	2.64	4.90	2.88	4.12	2.92
16384	20.06	13.27	18.69	14.05	17.69	13.52

Table 6. The

L^{\infty}

error

\|u_{R}^{j}-\phi_{j}\|_{L^{\infty}(\Omega)}

in the L-shaped domain on quasi-uniform meshes.

	$\epsilon=10^{-4}$		$\epsilon=10^{-2}$		$\epsilon=1$
$N$	Uniform	Shishkin	Uniform	Shishkin	Uniform	Shishkin
512	0.06	0.07	0.05	0.09	0.01	0.07
1024	0.14	0.11	0.11	0.12	0.05	0.11
2048	0.37	0.24	0.34	0.27	0.22	0.27
4096	1.15	0.72	1.13	0.78	0.91	0.85
8192	4.00	2.75	3.86	2.73	3.78	2.88
16384	17.34	14.11	17.27	13.04	17.00	14.23

5. Conclusion

It is worth noting that the proposed method can be applied to a wide range of singularly perturbed problems, not only the specific type of problem studied in this work. The decoupling technique presented in this paper can be adapted to other types of singularly perturbed problems, and the numerical scheme can be modified accordingly. Additionally, this work serves as a foundation for future research in the field of singularly perturbed problems, and can be used as a benchmark for comparing the performance of other numerical methods. Ultimately, we expect that it may be possible apply the proposed method for 2D with a completely new mesh algorithm as an alternative to a standard Shishkin mesh as shown in the figure 6. This is the subject of an ongoing work. Overall, the proposed method provides a promising approach for solving singularly perturbed problems efficiently and accurately.

Acknowledgments

The author thank the anonymous reviewers for their valuable time.

Funding

This research is supported by postdoctoral scholar program in the Department of Mathematics at the University of Kentucky in Lexington, Kentucky.

References

[1] Ehme, Jeffrey, Paul W. Eloe, and Johnny Henderson. ”Upper and lower solution methods for fully nonlinear boundary value problems.” Journal of Differential Equations 180.1 (2002): 51-64.
[2] Sun, Guangfu, and Martin Styne. ”Finite-element methods for singularly perturbed high-order elliptic two-point boundary value problems. II: convection-diffusion-type problems.” IMA journal of numerical analysis 15.2 (1995): 197-219.
[3] Amando, R. S. ”Numerical computation of internal and external flows, volume 1: Fundamentals of numerical discretization: C. Hirsch John Wiley & Sons, Ltd., 1988.” (1989): 371.
[4] Kreiss, Heinz-Otto, and Jens Lorenz. Initial-boundary value problems and the Navier-Stokes equations. Society for Industrial and Applied Mathematics, 2004.
[5] Alhumaizi, Khalid. ”Flux-limiting solution techniques for simulation of reaction–diffusion–convection system.” Communications in Nonlinear Science and Numerical Simulation 12.6 (2007): 953-965.
[6] Doolan, Edward P., John JH Miller, and Willy HA Schilders. Uniform numerical methods for problems with initial and boundary layers. Boole Press, 1980.
[7] Farrell, AJ. ”Sufficient conditions for uniform convergence of a class of difference schemes for a singularly perturbed problem.” IMA journal of numerical analysis 7.4 (1987): 459-472.
[8] Kadalbajoo, Mohan K., and Y. N. Reddy. ”Numerical treatment of singularly perturbed two point boundary value problems.” Applied mathematics and computation 21.2 (1987): 93-110.
[9] Kadalbajoo, Mohan K., and Y. N. Reddy. ”Approximate method for the numerical solution of singular perturbation problems.” Applied Mathematics and Computation 21.3 (1987): 185-199.
[10] Natesan, Srinivasan, and N. Ramanujam. ”A computational method for solving singularly perturbed turning point problems exhibiting twin boundary layers.” Applied Mathematics and computation 93.2-3 (1998): 259-275.
[11] OMalley Jr, Robert E. Introduction to Singular Perturbations. Volume 14. Applied Mathematics and Mechanics. New York Univ Ny Courant Inst of Mathematical Sciences, 1974.
[12] Hirsch, Charles. Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics. Elsevier, 2007.
[13] Ewing, Richard E., ed. The mathematics of reservoir simulation. Society for Industrial and Applied Mathematics, 1983.
[14] Kopteva, Natalia, and Martin Stynes. ”Approximation of derivatives in a convection–diffusion two-point boundary value problem.” Applied numerical mathematics 39.1 (2001): 47-60.
[15] Polak, S. J., et al. ”Semiconductor device modelling from the numerical point of view.” International Journal for Numerical Methods in Engineering 24.4 (1987): 763-838.
[16] Vrabel, R. ”On the lower and upper solutions method for the problem of elastic beam with hinged ends.” Journal of mathematical analysis and applications 421.2 (2015): 1455-1468.
[17] El-Zahar, Essam R. ”Approximate analytical solutions of singularly perturbed fourth order boundary value problems using differential transform method.” Journal of king Saud university-science 25.3 (2013): 257-265.
[18] Vrabel, R. ”Formation of boundary layers for singularly perturbed fourth-order ordinary differential equations with the Lidstone boundary conditions.” Journal of Mathematical Analysis and Applications 440.1 (2016): 65-73.
[19] Sun, Guangfu, and Martin Stynes. ”Finite-element methods for singularly perturbed high-order elliptic two-point boundary value problems. I: reaction-diffusion-type problems.” IMA journal of numerical analysis 15.1 (1995): 117-139.
[20] Roos, H. G., and M. Stynes. ”A uniformly convergent discretization method for a fourth order singular perturbation problem.” Bonner Math. Schriften 228 (1991): 30-40.
[21] Shanthi, V., and N. Ramanujam. ”Asymptotic numerical methods for singularly perturbed fourth order ordinary differential equations of convection–diffusion type.” Applied Mathematics and computation 133.2-3 (2002): 559-579.
[22] Miller, John JH, Eugene O’riordan, and Grigorii I. Shishkin. Fitted numerical methods for singular perturbation problems: error estimates in the maximum norm for linear problems in one and two dimensions. World scientific, 1996.
[23] Roos, Hans-Gorg, Martin Stynes, and Lutz Tobiska. Robust numerical methods for singularly perturbed differential equations: convection-diffusion-reaction and flow problems. Vol. 24. Springer Science & Business Media, 2008.
[24] Bakhvalov, Nikolai Sergeevich. ”On the optimization of the methods for solving boundary value problems in the presence of a boundary layer.” Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki 9, no. 4 (1969): 841-859.
[25] Kopteva, Natalia, and Eugene O’Riordan. ”Shishkin meshes in the numerical solution of singularly perturbed differential equations.” Int. J. Numer. Anal. Model 7.3 (2010): 393-415.
[26] Roos, Hans-Görg, Martin Stynes, and Lutz Tobiska. Robust numerical methods for singularly perturbed differential equations: convection-diffusion-reaction and flow problems. Vol. 24. Springer Science & Business Media, 2008.
[27] Ciarlet, Philippe G. The finite element method for elliptic problems. Society for Industrial and Applied Mathematics, 2002.
[28] Makungu, James, Heikki Haario, and William Charles Mahera. ”A generalized 1-dimensional particle transport method for convection diffusion reaction model.” Afrika Matematika 23 (2012): 21-39.
[29] Li, Hengguang, Charuka D. Wickramasinghe, and Peimeng Yin. ”A $C^{0}$ finite element method for the biharmonic problem with Dirichlet boundary conditions in a polygonal domain.” arXiv preprint arXiv:2207.03838, 2022.
[30] Wickramasinghe, Charuka Dilhara, ”A $C^{0}$ Finite Element Method For The Biharmonic Problem In A Polygonal Domain” (2022). Wayne State University Dissertations. 3704. https://digitalcommons.wayne.edu/oa_dissertations/3704.