Numerical Solution of Singularly Perturbed Problems via both Galerkin and Subdomain Galerkin methods

TL;DR

This paper presents numerical solutions for singularly perturbed boundary value problems using Galerkin and subdomain Galerkin finite element methods with quadratic B-splines, demonstrating their effectiveness through comparison with analytical solutions.

Contribution

It introduces the application of quadratic B-spline based Galerkin and subdomain Galerkin methods to singularly perturbed problems on graded meshes, providing new insights into their accuracy.

Findings

Methods accurately approximate analytical solutions

Subdomain Galerkin method improves convergence

Effective for problems with boundary layer behavior

Abstract

In this paper, numerical solutions of singularly perturbed boundary value problems are given by using variants of finite element method. Both Galerkin and subdomain Galerkin method based on quadratic B-spline functions are applied over the geometrically graded. Results of some text problems are compated with analytical solutions of the singularly perturbed problem

Equations

$$-\epsilon u^{\prime \prime }+p(x)u^{\prime }+q(x)u=f(x),x\in [0,1]$$

$$u\left(0\right)=\lambda \text{ and }u\left(1\right)=\beta ,\lambda ,\beta \in \mathbb{R}$$

$$0=x_0<x_1<\cdot \cdot \cdot <x_N=1$$

$$h_0+h_1+\cdot \cdot \cdot +h_{N-1}=1$$

$$h_0=\frac{1}{1+\sigma +{\sigma }^2+\cdot \cdot \cdot +{\sigma }^{N-1}}.$$

$$\begin{tabular}[]{l}$Q_{m-1}$\\ $Q_{m}$\\ $Q_{m+1}$\end{tabular}=\dfrac{1}{h_{m}^{2}}\left\{\ \begin{array}[]{l}(h_{m}-\xi)^{2}\sigma,\\ h_{m}^{2}+2h_{m}\sigma\xi-\left(1+\sigma\right)\xi^{2},\\ \xi^{2}\end{array}\ \right.$$

$$u_N=\sum _{m=-1}^N{\delta }_m{Q}_m$$

$$\begin{array}[]{ll}u_{m}=u(x_{m})&=\sigma\delta_{m-1}+\delta_{m},\vspace{0.2cm}\\ u_{m}^{\prime}=u^{\prime}(x_{m})&=\dfrac{2\sigma}{h_{m}}(\delta_{m}-\delta_{m-1}).\end{array}$$

$$-\epsilon v{u}^{\prime \prime }(x)+p(x)v{u}^{\prime }(x)+q(x)vu(x)=f(x)$$

$$\underset{x_m}{\overset{x_{m+1}}{\int }}(-\epsilon v^{\prime }{u}^{\prime }(x)+vp(x)u^{\prime }+vq(x)u(x))dx-\epsilon v{u}^{\prime }(x)\stackrel{x_{m+1}}{\underset{x_m}{|}}-\underset{x_m}{\overset{x_{m+1}}{\int }}vf(x)dx=0$$

$$\sum _{j=m-1}^{m+1}[-\epsilon \underset{0}{\overset{h_m}{\int }}({\varphi }_i^{\prime }{\varphi }_j^{\prime }+p{\varphi }_i{\varphi }_j^{\prime }+q{\varphi }_i{\varphi }_j)d\xi ]{\delta }_j-\epsilon {\varphi }_i{\varphi }_j^{\prime }\stackrel{h_m}{\underset{0}{|}}{\delta }_j-\underset{0}{\overset{h_m}{\int }}{\varphi }_{i}f(x_m+\xi )d\xi =0$$

$$\begin{array}[]{lll}a_{ij}=\mathop{\displaystyle\int}\limits_{0}^{h_{m}}\phi_{i}^{\prime}\phi_{j}^{\prime}d\xi,&&r_{ij}=\phi_{i}\phi_{j}^{\prime}\overset{h_{m}}{\underset{0}{\mid}},\\ &&\\ b_{ij}=\mathop{\displaystyle\int}\limits_{0}^{h_{m}}\phi_{i}\phi_{j}^{\prime}d\xi,&&f_{i}=\mathop{\displaystyle\int}\limits_{0}^{h_{m}}\phi_{i}f(x_{m}+\xi)d\xi,\\ \lx@intercol\hfil\hfil\lx@intercol&\lx@intercol\hfil\hfil\lx@intercol&\lx@intercol\hfil\hfil\lx@intercol\\ \lx@intercol\hfil c_{ij}=\mathop{\displaystyle\int}\limits_{0}^{h_{m}}\phi_{i}\phi_{j}d\xi,\hfil\lx@intercol&\lx@intercol\hfil\hfil\lx@intercol&\lx@intercol\hfil\hfil\lx@intercol\end{array}$$

$$A^{(m)}$$

$$B^{(m)}$$

$$C^{(m)}$$

$$R^{(m)}$$

$$F^{(m)}=\frac{1}{h_m^2}\left[\begin{array}{ccc}{\vartheta }_1& {\vartheta }_1& {\vartheta }_1\\ {\vartheta }_2& {\vartheta }_2& {\vartheta }_2\\ {\vartheta }_3& {\vartheta }_3& {\vartheta }_3\end{array}\right]$$

$${\vartheta }_1$$

$${\vartheta }_2$$

$${\vartheta }_3$$

$$(-\epsilon A^{(i)}+p{B}^{(i)}+q{C}^{(i)}-\epsilon R^{(i)}){\delta }^{(i)}=F^{(i)}$$

$${\delta }^{(i)}=({\delta }_{m-1}^{(i)},{\delta }_m^{(i)},{\delta }_{m+1}^{(i)},{\delta }_{mi+2}^{()}),F^{(i)}=(f_{m-1}^{(i)},f_m^{(i)},f_{m+1}^{(i)},f_{m+2}^{(i)})$$

$$(-\epsilon A+pB+qC-\epsilon R)\delta =F$$

$$A=\left[\begin{array}[]{cccccccccc}\sigma_{0,0}^{(0)}&\sigma_{0,1}^{(0)}&\sigma_{0,2}^{(0)}&&&&&&&\\ \sigma_{1,0}^{(0)}&\sigma_{1,1}^{\ast(1)}&\sigma_{1,2}^{\ast(1)}&\sigma_{1,3}^{\ast(1)}&&&&&&\\ \sigma_{2,0}^{(0)}&\sigma_{2,1}^{\ast(1)}&\sigma_{2,2}^{\ast(2)}&\sigma_{2,3}^{\ast(2)}&\sigma_{2,4}^{\ast(2)}&&&&&\\ &\sigma_{3,0}^{\ast(1)}&\sigma_{3,1}^{\ast(2)}&\sigma_{3,2}^{\ast(3)}&\sigma_{3,3}^{\ast(3)}&\sigma_{3,4}^{\ast(3)}&&&&\\ &&\ddots&\ddots&\ddots&\ddots&\ddots&&&\\ &&&\sigma_{i,i-2}^{\ast(i-1)}&\sigma_{i,i-1}^{\ast(i)}&\sigma_{i,i}^{\ast(i+1)}&\sigma_{i,i+1}^{\ast(i+1)}&\sigma_{i,i+2}^{\ast(i+1)}&&\\ &&&&\ddots&\ddots&\ddots&\ddots&\ddots&\\ &&&&&\sigma_{n-1,n-3}^{\ast(n-2)}&\sigma_{n-1,n-2}^{\ast(n-1)}&\sigma_{n-1,n-1}^{\ast(n)}&\sigma_{n-1,n}^{\ast(n)}&\sigma_{n-1,n+1}^{(n)}\\ &&&&&&\sigma_{n,n-2}^{\ast(n-1)}&\sigma_{n,n-1}^{\ast(n-1)}&\sigma_{n,n}^{\ast(n-1)}&\sigma_{n,n+1}^{(n)}\\ &&&&&&&\sigma_{n+1,n-1}^{(n)}&\sigma_{n+1,n}^{(n)}&a_{n+1,n+1}^{(n)}\end{array}\right]$$

$$\displaystyle\begin{array}[]{lll}\sigma_{1,1}^{\ast(1)}=\sigma_{1,1}^{(0)}+\sigma_{1,1}^{(1)},&&\sigma_{1,2}^{\ast(1)}=\sigma_{1,2}^{(0)}+\sigma_{1,2}^{(1)},\\ &&\\ \sigma_{2,1}^{\ast(1)}=\sigma_{2,1}^{(0)}+\sigma_{2,1}^{(1)},&&\sigma_{2,2}^{\ast(2)}=\sigma_{2,2}^{(0)}+\sigma_{2,2}^{(1)}+\sigma_{2,2}^{(2)},\\ &&\\ \sigma_{2,3}^{\ast(2)}=\sigma_{2,3}^{(2)}+\sigma_{2,3}^{(1)},&&\sigma_{i,i-1}^{\ast(i)}=\sigma_{i,i-1}^{(i-1)}+\sigma_{i,i-1}^{(i)},\end{array}$$

$$\displaystyle\begin{array}[]{lll}\sigma_{i,i}^{\ast(i)}=\sigma_{i,i}^{(i-1)}+\sigma_{i,i}^{(i)}+\sigma_{i,i}^{(i+1)},&&\sigma_{i,i+1}^{\ast(i)}=\sigma_{i,i+1}^{(i)}+\sigma_{i,i+1}^{(i+1)},\\ &&\\ \sigma_{n-1,n-2}^{\ast(n-1)}=\sigma_{n-1,n-2}^{(n-2)}+\sigma_{n-1,n-2}^{(n-1)},&&\sigma_{n-1,n-1}^{\ast(n)}=\sigma_{n-1,n-1}^{(n-2)}+\sigma_{n-1,n-1}^{(n-1)}+\sigma_{n-1,n-1}^{(n)},\\ &&\\ \sigma_{n-1,n}^{\ast(n)}=\sigma_{n-1,n}^{(n-1)}+\sigma_{n-1,n}^{(n)},&&\sigma_{n,n-1}^{\ast(n-1)}=\sigma_{n,n-1}^{(n-1)}+\sigma_{n,n-1}^{(n)},\\ &&\\ \sigma_{n,n}^{\ast(n-1)}=\sigma_{n,n}^{(n-1)}+\sigma_{n,n}^{(n)}.&&\end{array}$$

$$F=\left[\begin{array}[]{l}f_{0}^{0}\\ f_{1}^{0}+f_{1}^{1}\\ f_{2}^{0}+f_{2}^{1}+f_{2}^{2}\\ \vdots\\ f_{i}^{i-1}+f_{i}^{i}+f_{i}^{i+1}\\ \vdots\\ f_{n-2}^{n-1}+f_{n-2}^{n-2}+f_{n-2}^{n-3}\\ f_{n-1}^{n-1}+f_{n-1}^{n-2}\\ f_{n}^{n-1}\end{array}\right]$$

$${\delta }_{-1}=\frac{\lambda -{\delta }_0}{\alpha },{\delta }_N=\beta -\alpha {\delta }_{N-1}.$$

$$\underset{x_0}{\overset{x_n}{\int }}[-\epsilon u^{\prime \prime }(x)+p{u}^{\prime }(x)+qu(x)-f(x)]dx=0$$

$$V_{n}=\left\{\begin{array}[]{cc}1,&x_{m}\leq x<x_{m+1}\\ 0,&\text{other case}\end{array}\right.$$

$$-\epsilon u^{\prime }(x)\stackrel{x_{m+1}}{\underset{x_m}{|}}+pu(x)\stackrel{x_{m+1}}{\underset{x_m}{|}}+q\underset{x_m}{\overset{x_{m+1}}{\int }}u(x)dx=\underset{x_m}{\overset{x_{m+1}}{\int }}f(x)dx$$

$$\begin{array}[]{l}[-\varepsilon(\mathop{\displaystyle\sum}\limits_{J=m-1}^{m+1}\phi_{i}^{\prime}\overset{hm}{\underset{0}{\mid}})\delta_{j})+p(x_{m}+\xi)(\mathop{\displaystyle\sum}\limits_{J=m-1}^{m+1}\phi_{i}\overset{hm}{\underset{0}{\mid}}\delta_{j})+q(x_{m}+\xi)\mathop{\displaystyle\sum}\limits_{J=m-1}^{m+1}\mathop{\displaystyle\int}\limits_{x_{m}}^{x_{m+1}}\phi_{i}\delta_{j}dx\\ =\mathop{\displaystyle\int}\limits_{x_{m}}^{x_{m+1}}f(x_{m}+\xi)dx\end{array}$$

$$\underset{x_m}{\overset{x_{m+1}}{\int }}{u}^{\prime \prime }(x)dx$$

Full text

**Keywords: **Subdomain Galerkin, graded mesh, spline, singularly perturbed.

1 Introduction

This paper contains numerical solutions of one dimensional singularly perturbation problems

[TABLE]

with boundary conditions

[TABLE]

where $\varepsilon$ is a small positive parameter, $p(x),$ $q(x),$ $f(x)$ are sufficiently smooth functions with $p(x)\geq p^{\ast}>0,$ $q(x)\geq q^{\ast}>0.$ These problems depend on $\varepsilon$ in such a way that the solution varies rapidly in some parts and varies slowly in some other parts. So, typically there are thin transition layers where the solutions can jump abruptly, while away from the layers the solution behaves regularly and vary slowly. The numerical treatment of the singular perturbation problems is far from the trivial because of the boundary layer behavior of the solution. There are a wide variety of asymptotic techniques for solving singular perturbation problems.

These problems occur in many areas of engineering and applied mathematics such as chemical reactor theory, optimal control, quantum mechanics, fluid mechanics, reaction-diffusion process, aerodynamics, heat transport problems with large Peclet numbers and Navier–Strokes flows withlarge Reynolds numbers etc.

Many authors have studied on this problem and tried to overcome the above-mentioned difficulties. M.K. Kadalbajoo and Vikas Gupta [6] proposed B-spline collocation method on a non-uniform mesh of Shishkin type.to solve singularly perturbed two-point boundary value problems with a turning point exhibiting twin boundary layers. J.Vigo-Aguiar and S.Natesan [1] consider a class of singularly perturbed two-point boundary-value problems for second-order ordinary differential equations. They suggested an iterative non-overlapping domain decomposition method in order to obtain numerical solution to these problems. Tirmizi et al. [11] have proposed a generalized scheme based on quartic non-polynomial spline functions in order to designed for numerical solution of singularly perturbed two-point boundary-value problems. D.J.Fyfe [5] used cubic splines on equal and unequal intervals and compared the results. He observed that very little advantage is gained by using unequal intervals. M.K.Kadalbajoo and K.C.Patidar [7] gave some difference schemes using spline in tension. They showed that these methods are second-order accurate. Employing coordinate stretching a Galerkin-spectral method is applied to the singularly perturbed boundary value problems by W.Liu and T.Tang [8]. G.Beckett and J.A.Mackenzie [2] gave a $p$th order Galerkin finite element method on a non-uniform grid. In their study the grid is constructed by equidistributing a strictly positive monitor function. After the appropriate selection of the monitor function parameters they obtained insensitive numerical solution.

The definitions of B-splines over the geomertically graded mesh was given in reference [3]. Dag and Sahin [4] have set up the finite element method employing the quadratic and the cubic B-splines to form the trial function. In this article, we used the finite element method with the quadratic B-splines. After giving the expressions of the mentioned B-splines over the geometrically graded mesh we applied the quadratik Galerkin and quadratik subdomain Galerkin method to Eq.(1).

Briefly, outline is as follows. In Section 2, numerical methods are given. Numerical experiments are carried out for one test problem and errors of those methods are compared in Section 3. Finally conclusion is given in last section.

2 B-spline Galerkin Methods

For numerical purpose, let us divide the solution domain $\left[0,1\right]$ into subintervals by the knots $x_{m}$ such that

[TABLE]

where $x_{m+1}=x_{m}+h_{m}$ and $h_{m}$ is the size of interval $\left[x_{m},x_{m+1}\right]$ with relation $h_{m}=\sigma h_{m-1}.$ Here $\sigma$ is mesh ratio constant.

To construct the geometrically graded mesh, determination of the first element size $h_{0}$ is necessary. Since

[TABLE]

it is easy to write

[TABLE]

This partition will be uniform if the mesh ratio $\sigma$ is taken as unity. To obtain finer mesh at the left boundary, $\sigma$ must be chosen as $\sigma>1.$ On the other hand, to make the mesh size smaller at the right boundary, $\sigma$ must be chosen as $\sigma<1$ . Mentioned selection of $\sigma$ will be done by experimentally.

2.1 Quadratic B-spline Galerkin method (QM)

The expression of the quadratic B-splines over the geometrically graded mesh may be given in the following form [3]:

[TABLE]

where $\xi=x-x_{m}$ and $0\leq\xi\leq h_{m}.$ A quadratic B-spline covers $3$ elements. Any quadratic B-spline $Q_{m}$ and its derivatives vanish outside of the interval $\left[x_{m-1},x_{m+2}\right]$ and therefore an element is covered by $3$ successive quadratic B-splines. The set of the quadratic B-splines $\left\{Q_{-1},Q_{0},...,Q_{N}\right\}$ forms a basis for the functions defined on the solution domain [10]. Thence, an approximation $u_{N}$ to the analytical solution $u$ can be written as

[TABLE]

where $\delta_{m}$ are unknown parameters. By the substitution of the value of $Q_{m}$ at the knots $x_{m}$ in Eq.(4), the nodal value $u$ and its derivative $u^{\prime}$ are expressed in terms of $\delta_{m}$ by

[TABLE]

Both sides of the weight function by multiplying the differential equation and is integrated over the range $[x_{m},x_{m+1}]$ following equation is obtained.

[TABLE]

partial integration is applied to the first term in the above integral is obtained as follows.

[TABLE]

weight functions selected as $Q_{j},$ $j=m-1,m,m+1$ and used (4), following integral is obtained.

[TABLE]

So, following values can be computed.

[TABLE]

Where $i,j=m-1,m,m+1$ and

[TABLE]

and

[TABLE]

Where

[TABLE]

Defined in terms of local matrices $A^{(i)},$ $B^{(i)},$ $C^{(i)},$ $R^{(i)}$ and $F^{(i)}$, equation (6) can be represented following form.

[TABLE]

Where

[TABLE]

By combining the local matrices which is defined on $[x_{i},x_{i+1}],$ $i=0,\ldots,N-1$, the global system in the range of $[x_{0},x_{N}]$ can be defined as follows.

[TABLE]

The matrix of $A$ is

[TABLE]

Where

[TABLE]

$B,$ $C,$ $R$ matrices are obtained similarly. Also the matrix of $F$ is computed as follows.

[TABLE]

The matrix system (7) has $N+1$ equations and $N+3$ unknowns. In order to solve this system, the numbers of equations and unknowns must be equal. From the boundary conditions (2) and Eq.(5) it is easy to write

[TABLE]

Using these equalities, $\delta_{-1}$ and $\delta_{N}$ can be eliminated from the system and then matrix equation (15) can be solved with Thomas algorithm. Substituting the obtained parameters $\delta_{m}$ in Eq.(5), the numerical solution is found at the knots $x_{m}.$

2.2 Quadratic B-spline Subdomain Galerkin method (QM)

If on each side of the equation (1), multiplied by weight function $V_{n}$ and is integrated over the range $[x_{m},x_{m+1}]$ then following integral obtained.

[TABLE]

Here the wight function is

[TABLE]

partial integration is applied to the first term in the above integral is obtained as follows.

[TABLE]

By the substitution of the nodal value $u$ and its derivative $u^{\prime}$ in last equation is expressed following integral.

[TABLE]

With the help of the division points values, Qauadratik B-spline shape functions defined on geometrically increasing intervals $[x_{m},x_{m+1}]$

[TABLE]

The value of $Q_{m-1,}Q_{m},$

$Q_{m+1}$ substitute in (13) is computed as follows.

[TABLE]

If we replace integrals calculated by (13) the following results are obtained.

[TABLE]

When we apply the Galerkin method to the system (10), (11), (12) and (13) derivatives are replaced with their equals obtained from Eq.(10). This substitution yields the following system:

[TABLE]

With necessary operations, this system can be written in matrix form as

[TABLE]

where

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

The matrix system (14) has $N+1$ equations and $N+3$ unknowns. In order to solve this system, the numbers of equations and unknowns must be equal. From the boundary conditions (2) and Eq.(5) it is easy to write

[TABLE]

Using these equalities, $\delta_{-1}$ and $\delta_{N}$ can be eliminated from the system and then matrix equation (15) can be solved with Thomas algorithm. Substituting the obtained parameters $\delta_{m}$ in Eq.(5), the numerical solution is found at the knots $x_{m}.$

3 Numerical Experiments

We have tested the accuracy of the numerical methods on two examples. Errors are measured with the norm

[TABLE]

Since the boundary layers are at the right boundary in both examples, in order to minimize the error, we have searched the interval $\left(0,1\right)$ for the best choice of  the mesh ratio $\sigma.$ Solution profiles are illustrated in Figs. 1-4 for the first example and in Figs. 5-8 for the second example. These figures are graphed for $N=20$ and two different $\varepsilon.$ In order to see the success of the numerical methods more clear, exact solutions and obtained results are illustrated together in all figures. In all figures, continuous line is used for the exact solutions and the lines$\cdots$o$\cdots$o$\cdots$, $\cdots$+$\cdots$+$\cdots$ are used for QM and CM respectively. Using uniform mesh leads to oscillations, seen in Figs. 1, 3, 5 and 7, in solution profiles because of the boundary layer. As observed from Figs. 2, 4, 6 and 8, after the best choice of mesh ratio $\sigma,$ these oscillations disappear. Using various $\varepsilon$ and $N,$ calculated numerical errors are tabulated and compared in Table 1 and Table 2 for the first and the second examples respectively.

**Example **Our example is

[TABLE]

with the exact solution

[TABLE]

taken from [9].

4 Conclusion

Quadratic and cubic B-spline algorithms are applied to singularly perturbed problems. Difficulties arised from the modelling of the boundary layers in numerical methods are tried to overcome by using B-splines over the geometrically graded mesh. Simplicity of the adaptation of B-splines and obtaining acceptable solutions can be noted as advantages of given numerical methods. Consequently, in getting the numerical solution of the differential equations having boundary layers, B-spline collocation methods over the geometrically graded mesh are advisable.