Richardson extrapolation technique for singularly perturbed convection-diffusion problem with non-local boundary conditions

This study presents the Richardson extrapolation techniques for solving singularly perturbed convection-diffusion problems (SPCDP) with non-local boundary conditions (NLBC). A numerical approach is presented using an upwind finite difference scheme a piecewise-uniform (Shishkin) mesh. To handle the non-local boundary conditions, the trapezoidal rule is applied. The study establishes an error bound for numerical solutions and determines the numerical approximation for scaled derivatives. To enhance convergence and accuracy, we utilize Richardson extrapolation. This elevates accuracy from O (cid:0) N − 1 ln N (cid:1) to O (cid:0) N − 2 ln 2 N (cid:1) using this technique, where N is the number of mesh intervals. Numerical results are presented to validate the theoretical findings, demonstrating the effectiveness and accuracy of the proposed technique.


Introduction
This article investigates the convection-diffusion type's singularly perturbed differential equations (SPDE) with non-local boundary conditions.Non-local boundary conditions are indeed necessary and commonly employed when studying SPCDP [1,2,4,20].Such boundary value problems (BVP) with integral boundary conditions (IBC) are commonly encountered in various domains, including electrochemistry, thermoelasticity, heat conduction, etc.Only a few authors [1,9,17] have addressed SPDE with non-local boundary conditions.Singular perturbation problems (SPPs) arising in the context of convection-diffusion problems with non-local boundary conditions are prevalent across various fields of applied mathematics and engineering [3].BVPs with integral boundary conditions multiplying the leading derivative term by a small parameter ε are called singularly perturbed problems with integral boundary conditions.SPCDP exhibits a characteristic where a small parameter ε(0 < ε ≪ 1) multiplies some or all of the highest-order terms in the differential equation.This parameter represents the degree of perturbation in the system and significantly influences the solution's behavior.The Navier-Stokes equation, governing fluid flow dynamics, is an example of a SPCDP [30].Introducing the small parameter ε, the equation becomes: where u and v represent velocity components along x and y directions, respectively, and p denotes pressure.The Reynolds number Re is a dimensionless parameter relating to the fluid's length scale, velocity scale, and kinematic viscosity.At large Reynolds numbers (Re), the equation transforms into a SPCDP with non-local boundary conditions.The diffusion terms ∂ 2 u ∂x 2 and ∂ 2 u ∂y 2 are scaled by ε 2 , indicating their relatively smaller influence compared to the convective terms ∂(u 2 +p) ∂y and ∂(uv) ∂y .The small parameter ε signifies the existence of a boundary layer where the solution exhibits rapid variations.SPCDPs involve significant contributions from both convection and diffusion, with ε amplifying either the convection or diffusion term in the equations.Finite Difference Methods (FDMs) are commonly used to approximate such solutions, though research on approximating their derivatives has been relatively limited.These approximations are valuable in certain applications like flux or flag calculations.SPPs, characterized by a small parameter ε(0 < ε ≪ 1) multiplying the highest derivative term, have been extensively studied in the field of differential equations [14,28].These problems exhibit rapid changes in the solution within specific domain regions.To obtain accurate numerical solutions for such problems, it is essential to develop appropriate approaches that provide error estimates independent of the small parameter.One of the most straightforward and practical approaches for developing such methods involves employing a category of piecewise uniform (Shishkin) mesh.Numerical methods for equations with non-local boundary conditions have also been widely investigated [2].
Richardson extrapolation has been employed by various scholars as a technique to solve SPCDPs with non-local boundary conditions.Relevant works include what is mentioned in the reference (for example [10,13,19,22,31]).The primary goal of this work is to present Richardson's extrapolation to improve numerical solution accuracy and efficiency, as well as to investigate and analyze and post-processing method.This improvement in convergence rate is particularly targeted for problems that are discretized using a Shishkin mesh [24].The development of SPCDPs with non-local boundary conditions arises from the need to accurately model physical phenomena exhibiting both convective and diffusive behavior while considering non-local effects at the boundaries.Studying SPCDPs with non-local boundary conditions aims to accurately represent specific systems' behavior, enhancing the methods' accuracy, stability, and efficiency.The use of these techniques improves the effectiveness of investigating real-world events and allows for consistent results in a wide range of applications [7].
M. Cakir and G. M. Amiraliyev [4,5] developed a second-order numerical method for SPPs with non-local boundary conditions, and investigated FDM for the same problem with non-local boundary conditions.Amiraliyev and Raja [1], focus on the well-posedness of SPDEs with non-local boundary conditions.These works will investigate solutions' existence, uniqueness, and stability to ensure that the issue is well-posed.On the other hand, Kopteva and Stynes [16,24,31] focus on obtaining derivative approximations in SPCDPs that consider scale variations between the boundary layer region and the outer region.Their research addresses SPCDPs where the convection and diffusion terms exhibit significantly different magnitudes.They aim to accurately capture the solution behavior in different problem regions by appropriately scaling the derivatives.As the ε(0 < ε ≪ 1) approaches zero, the derivative solution of SPDE becomes unbounded.To approximate these derivatives accurately, scaling techniques are necessary.
The research work of R. Mythili Priyadharshini and N. Ramanujam primarily revolves around the development of approximation techniques for computing derivatives, particularly scaled first and second derivatives.Their contributions are highlighted in the papers cited as [23,25].One of the main objectives of their research is to address the numerical approximation of derivatives, especially in the context of problems with disparate scales or SPP.The research conducted by Debela and Duressa [11] presents a significant advancement in the field of numerical methods for solving SPCDPs with non-local boundary conditions.Later, Debela and Duressa proposed a computational method for the class of SPCDE with IBC using the Richardson extrapolation technique in [9].Motivated by these considerations, we proposed using the Richardson extrapolation technique to solve SPCDPs with non-local boundary conditions (2.1).To the best of our knowledge, the authors' approach began with proposing an approximation method for scaled solution derivatives.We worked on problems with upwind schemes on a Shishkin mesh using scaled derivatives, obtaining approximately a first-order convergence rate.The authors most likely used Richardson extrapolation on the upwind FDM within the Shishkin mesh to improve accuracy from almost first to almost second order.The rest of the article is organized as follows: Section 2, defines problems with non-local boundary conditions and presents analytical results.Section 3, introduces a numerical method based on an upwind FDM.Section 4, discusses numerical scaled derivatives.Section 5, includes numerical examples to validate the theoretical results.In Section 6, the paper provides final conclusions.Throughout the paper, C refers to a generic constant independent of ε and discretization parameters N .

Problems with non-local boundary conditions and some analytical results
Find ξ such that: where 0 < ε ≪ 1 is a singular perturbation parameter and the coefficient functions µ(x), η(x) are smooth, bounded and satisfy µ(x) and it satisfies The above problem satisfies the maximum principle and stability result.The detail proof is given in [1,4,27].

Analytical results
To develop sharp bounds we write the analytical solution in the form ξ(x) = p(x) + q(x), where p(x) is the smooth component and q(x) is the singular component.The smooth component p(x) can be expressed as an asymptotic expansion p(x) = p 0 (x) + εp 1 (x) + ε 2 p 2 (x) satisfies the following equations Hence, the smooth component of the solution satisfies (2.5) Now, we must extend the scaling to the point where x = 1.Consequently, the differential operator Therefore, the boundary component of the solution satisfies the following conditions: we choose q(z) = m+1 i=0 ε i q i (z).
Theorem 2.1.Let ξ(x) be the solution of (2.1) and p 0 (x) be the solution of (2.5).Then, there exists a constant C > 0 such that for all x ∈ Ω, we have Proof.The detailed proof is given in [21].
Lemma 2.2.The solution ξ can be decomposed into the sum ξ = p + q, where p(x) is the smooth component and q(x) is the boundary component respectively.Furthermore, these components and their derivatives satisfy the following bounds: ) (2.9) Proof.The detail proof is given in [4,12].The detailed proof of the lemma can be established by utilizing appropriate barrier functions, making use of Theorem 1, and employing the proof technique described in the reference [12] (p.46).

Numerical Methods
This section discusses the mesh selection strategy for solving the problem (2.1).We concentrate on the piecewise uniform mesh.The numerical computations use an upwind FDM, which accounts for non-local boundary conditions.

Construction of piecewise-uniform (Shishkin) mesh
This mesh is extensively discussed in the references [8,13,18,24].To effectively handle the boundary layer at x = 1 in the SPCDPs with non-local boundary conditions (2.1), we utilize a piecewise-uniform mesh.This mesh includes a transition point at 1 − τ , where To ensure numerical solvability, we choose the parameter τ based on a specific condition: where N is significantly larger than 1 ε .A piecewise-uniform mesh is constructed by dividing the domain Ω = [0, 1] into two sub-intervals: [0, 1 − τ ] and [1 − τ, 1].Each sub-interval is uniformly subdivided into N/2 intervals to create the mesh.The mesh widths h i = x i − x i−1 .Now, according to the definition of x i 's, the spatial mesh sizes can be expressed as follows: The piecewise-uniform mesh, denoted as Ω N 1 , is entirely defined by user-specified parameters N and τ .The interior mesh points are given by: The step between consecutive interior mesh points is given by ĥ = h i+1 + h i .

Numerical Scheme
The discrete problem corresponding to (2.1) is as follows: Find Z N (x i ) such that where the first and second-order finite differences are defined as

Note:
The above numerical scheme satisfies the discrete maximum principle and discrete stability result.Hence, the matrix associated with this scheme is an M-Matrix.
Equation (3.4) can be expressed as the following system of algebraic equations: where the coefficients in the upwind finite difference scheme are given by . (3.6) Theorem 3.1.Let ξ be the solution of (2.1) and Z N be the numerical solution defined by (3.4).
Then, the following inequality holds: where C is a constant independent of ε and N .
Proof.The detailed proof can be found in [26].
Prior to the extrapolation analysis, we introduce a crucial lemma for the subsequent section.We define the piecewise (0, 1)-Padé approximation of exp −αx i ε on the mesh Ω N 1 , where i = 0, 1, . . ., N , as the following mesh functions: ,where by convection S 0 = S ′ 0 = 1.The numerical scheme for an upwind scheme (3.4) on Ω N can be expressed as follows: For i = 1, . . ., N − 1, Lemma 3.2.The mesh functions S i satisfy the following property for i = 1, . . ., N −1: there exists a positive constant C such that Furthermore, for i = N/2 + 1, . . ., N − 1, there exists a constant C 1 such that Proof.The detailed proof can be found in [24] 3.

Richardson extrapolation Technique
This article aims to enhance the accuracy of the upwind scheme (3.4) using Richardson extrapolation, which has proven to improve numerical solutions for differential equations [29].The ξ N (x) is computed on the mesh Ω N 1 and Ω 2N 1 .This mesh has 2N sub-intervals and the transition point 1 − τ as Ω N 1 .The two meshes are interconnected, that is, Clearly Ω 2N  1 be the mesh obtained by bisecting the mesh intervals in Ω N 1 , and let ξ N represent the approximation of the solution on Ω 2N  1 .Thus, where C is a constant independent of mesh size h i and ε.The remainder term R N (x i ) ≤ Ch 2 i .Now, from the definition of transition parameters τ = 2ε α ln N and ln N = ατ 2ε .Substituting ln N into (3.9),we find that Ω2N 1 is solved using the same 1 − τ transition point: where Z 2N denotes the solution of discrete problem (3.4) and where the remainders Subtracting equation (3.11) from equation (3.9) or eliminating the initial term O(N −1 ) from both equations, we acquire: The truncation error of the remainder R N originates from combining the first-order derivative difference scheme and the trapezoidal rule's truncation error for a non-uniform mesh.Thus, the approximation for the truncation error of the remainder R N from equation (3.12) is as follows: This suggests that the Richardson extrapolation technique enhances the convergence rate from almost first-order to almost second-order.

Discrete solution decomposition
Analogous to the continuous solution, we can decompose the discrete solution Z N into the following sum , where P N represents the smooth component and satisfies the following discrete problems: IJFMR23069377 Volume 5, Issue 6, November-December 2023 and Q N represents the boundary component, and satisfies the following discrete problem: Following a similar approach, we define Z 2N (x i ) = P 2N (x i ) + Q 2N (x i ) to continue the decomposition.Thus, the error in the discrete solution is presented by decomposing the solution Z N over Ω N 1 as described in equations (3.13) and (3.14):

Extrapolated solution of P N
Let's first consider the error in the smooth part of ξ, denoted by P N (x i ) − p(x i ).The lemma provides the bound for the truncation error p. Lemma 3.3.Assuming ε ≤ N −1 and for all x i ∈ Ω, we can calculate the local truncation error for the smooth component p in the following way: Proof.Based on the bounds established in [21], where

and all
x ∈ [0, 1], we can use Lemma 2.2 for derivative bounds on p, along with Taylor's series expansion, to derive the following equation: where For every x in Ω, Ẽ is defined as a non-local solution to the boundary value problem using Keller's classical approach [15].The function Ẽ is defined as the solution to the following BVP: where Φ(x) is given by Φ where ψ and λ represent the smooth and boundary parts of Ẽ, respectively.Now by using (2.8), we have the following bounds: (3.17) Therefore, we have shown that Lemma 3.4.Under assumptions ε ≤ N −1 .We have Proof Given x i ∈ (0, 1) fixed, Taylor's expansion yields: Furthermore, using (3.16), we find: From the truncation error, we have: Thus, Now, let us define discrete mesh functions Then, apply the difference operators on M i by using Lemma 3.2 and ε ≤ N −1 for 0 ≤ i ≤ N/2, it follows that and for By choosing an adequately large value for C 3 , it effectively operates as a barrier function denoted by M i for ± P N (x i ) − p(x i ) − Hψ(x i ) .Employing the discrete maximum principle to the barrier function M i , we subsequently ascertain: Thus, for i = 1, . . ., N/2, we have: where we used Ẽ = ψ + λ.Now we can show that extrapolation improves the accuracy of P N on (0, 1 − τ ] interval.
IJFMR23069377 Volume 5, Issue 6, November-December 2023 Lemma 3.5.For all x i ∈ [0, 1 − τ ], we have achieved second-order convergence for the smooth component |p( where C is a constant. Since the subinterval mesh width of Ω 2N 1 is half of those of Ω N 1 , and the function Ẽ(.) depends on τ , we utilize Lemma 3.4 to deduce: Similarly, by maintaining a constant value for τ on the mesh Ω 2N  1 , we obtain: Using the extrapolation formulas (3.12), 3.4, and (3.23), we arrive to the following conclusion: The subsequent lemma demonstrates the error of P N (x i ) after extrapolating over (1 for some constant C. Proof.We define the function Ĝ(x) on [1 − τ, 1] by Within the domain Ω N 1 , we introduce a discrete approximation ĜN of Ĝ as outlined below: For the convergence of the upwind scheme, we have Ĝ We define φ(x Using a barrier function of the form C(1 + x i )N −1 (ln N ) 2 , we obtain the following estimate for φ(x i ): Furthermore, noticing that G(1 − τ ) ≤ C, we can deduce the following relationship: Similarly, on the mesh Ω 2N , one has IJFMR23069377 Volume 5, Issue 6, November-December 2023 E-ISSN: 2582-2160 Website:www.ijfmr.comEmail: editor@ijfmr.com Proof.Referring to (Theorem 1 of Cen [6]), we obtain: To study the impact of extrapolation (1 − τ, 1), we introduce the function F over the interval [1 − τ, 1] through a BVP with nonlocal boundary conditions.Suppose F is the solution of the following BVPs with IBC: (3.28) Then, F depends upon τ and independent of N .Now using the fact that ||F ′ (0)|| ≤ Cε −1 , we have Lemma 3.8.For all x i ∈ [1 − τ, 1], we have Proof.The detailed proof is given [24,31].
We will now illustrate how extrapolation improves the accuracy of Q N (x i ) for x i within the range of Lemma 3.9.For some constant C and for all Proof.Assuming that x i ∈ [1−τ, 1], we can reconfigure equation (3.30) to more explicitly showcase its reliance on the selection of both N and τ .
Given that Ω2N 1 is solved utilizing the identical transition point 1 − τ , it can be concluded that: Replacing the first term (3.31) and (3.32), we have Hence, the desired result is required.
Theorem 3.10.(Error after extrapolation) Assume that ε ≤ N −1 , then there exists a positive constant C such that: where C are positive constants.Proof.
For each We obtain the desired result by combining the results of Lemmas 3.5, 3.6 for the smooth component and 3.7, 3.8, 3.9 for the layer component of the above equations.

Numerical Scaled Derivatives
In this section, we establish the approximation of the scaled derivative εξ ′ (x) of the continuous solution (2.1), as well as the scaled upwind discrete derivative εD − Z(x) of the numerical solution (3.4), at internal points x i (i = 1, . . ., N − 1) within both the fine and coarse regions.The errors, represented as e(x i ) = Z(x i ) − ξ(x i ), satisfy the following equation: where, according to Theorem 3.1, [η(x i )e(x i )] = O(N −1 ln N ).We utilize the equation above in the proofs of the subsequent lemmas and theorems.
Lemma 4.1.For each mesh point x i ∈ Ω N and all x ∈ Ωi = [x i−1 , x i ], we have: where ξ is the solution of any problem given in (2.1).Proof.To bound D − ξ(x i ) − ξ ′ (x) , we can express it as follows: we can precisely bound the expression by considering each term independently on the right-hand side (RHS).Using the triangle inequality, we can write: For any mesh point x i ∈ Ω N ∪ {0}, and any function φ ∈ C 2 ( Ωi ) that satisfies the given differential equation, we have the following identity: from which it follows that We apply to this p and, using Lemma (2.2) and the solution's decomposition, we observe that this yields the required bound in the first term.We have plans for the second term To bound each term separately, consider the case τ = 1/2 where When τ = 1/2, we see from the argument leading to (4.4) that For all x i ∈ [0, 1 − τ ), we have that Hence, let's focus on the remaining case τ = 2ε α ln N and x i ∈ [1 − τ, 1).Then, we can use the fact that Lq = 0 and integration by parts, we obtain: From this, we can conclude that for all s ∈ Ωi : Combining these inequalities completes the proof.
Lemma 4.2.Let p and P N be the exact and discrete regular components of the solution to the problem (2.5) and (3.13) respectively.Then for every mesh point x i ∈ Ω N , we have Proof.We denote the notation of error and truncation error at each mesh point by First, we proved that for all i, 0 To establish the outcome for N/2 + 1 ≤ i ≤ N − 2, we express the relationship τ (x i ) = L N e(x i ) in the following form: By summing and rearranging these equations from j = i o N − 1, we obtain Proceeding to bound each term separately, we have already established that the first term is bound as For the second term, since |τ To bound the last term, let's consider the telescoping effect.We observe as follows: and .
By summing and rearranging the expressions, we see the first RHS bracket term simplify through telescoping, and the last RHS bracket is non-zero solely for j = N/2 and j = N/2 + 1. Therefore,for 1 < i ≤ N/2 we have: Similarly for i ≥ N/2, all terms on the RHS of the expression are zero.By utilizing the given facts: Using these facts, we deduce that for all i, Finally, considering the boundedness of µ(x j ) and e(x), and the fact that x i+1 − x i = O(N −1 ), we conclude the expression is bounded 0 ≤ i ≤ N − 2 as follows: Lemma 4.3.Let q be the exact solution of (2.6) and let Q N be the discrete singular component of the solution of Equations (3.14).Then, when τ = 1 2 , we have for all where C is a constant.When τ = 2ε α ln N , we have ) Proof.For detailed proof, one can refer to ([21], Lemma 3.15).
Lemma 4.4.Let q and Q N be the exact and discrete singular components, respectively, of the solution defined by (2.6) and (3.14).Then, for all x i ∈ Ω N , we have where C is a constant independent of N and ε.
Proof.We represent the error and truncation error at each mesh point by We will prove the result separately for two cases: τ = 1 2 and τ = 2ε α ln N .The local truncation error is bound in the standard way and using Lemma 2.2, we obtain First of all, we consider the case τ = 2ε α ln N , the mesh is non-uniform.For all , we can use the triangle inequality to obtain: We bound each term on the right-hand side individually.According to lemma 4.1, this yields: Now, we bound each term separately.In the first case x i ∈ ΩN ∩ (0, 1 − τ ].Therefore, to bound the first term on the right-hand side, IJFMR23069377 Volume 5, Issue 6, November-December 2023 By utilizing Lemma 4.3 and (2.9), along with the triangle inequality, we can express this as: For For 1 ≤ i < N/2 − 1, we write the equation τ ( (4.13) Summing and rearranging we get As a result, using the result at x N/2 + 1 and (4.12), we obtain 2 , using the above arguments along with (4.12) and (4.13), we arrive at the result: Thus, Lemma 4.4 has been successfully proven.
IJFMR23069377 Volume 5, Issue 6, November-December 2023 Theorem 4.5.Let ξ be the solution of the continuous problem (2.1), and let Z be the solution of the discrete problem (3.4).For any i within the range 0 ≤ i ≤ N − 1, the following ε-uniform error estimate satisfies: where C is independent of ε and N .

Numerical Examples, Results and Discussion
We perform numerical tests to confirm theoretical findings, employing model problems from equations (2.1) and utilizing the numerical scheme in equation (3.4).This section presents two examples, and since exact solutions are unknown, we assess the maximum point-wise error using the double mesh principle [2,24,26].The calculation of the maximum pointwise double-mesh differences is given as follows: where Z N and Z 2N represent the numerical solutions acquired with N and 2N mesh intervals, respectively.We compute the ε-uniform maximum pointwise double-mesh differences, denoted as We define the computed corresponding ε-uniform numerical convergence rates for all N as: The numerical results are shown for ε values from ε ∈ R ε = 2 −20 , . . ., 2 −2 , 2 −1 , where R ε define the ranges for the singular perturbation parameter.
Example 5.1.Consider the following SPP: with boundary conditions Example 5.2.Consider the following SPP:      Using the provided data, Tables 1 and 2 show calculated values of E N and S N for the scaled derivative of the solution εξ ′ .Computation utilizes scaled discrete upwind method, shown in Examples 5.1 and 5.2.In Tables 3 and 4, we can find precise maximum pointwise errors and convergence rates for Examples 5.1 and 5.2.These tables indicate nearly first-order convergence.Additionally, the tables summarize the maximum pointwise errors and convergence orders for the examples.As we review the results in Tables 3 and 4, we will notice a consistent decrease in the computed ε-uniform errors E N for Examples 5.1 and 5.2 as N increases.This confirms the ε-uniform convergence of the upwind scheme (3.4) both before and after extrapolation.We can see from the numerical solution plots in Figures 1 and 3, as well as the log-log plots of maximum pointwise errors in Figures 2 and 4, that Richardson extrapolation effectively increases the order of convergence of the upwind scheme.The upwind scheme's order of convergence improves from O N −1 ln N to O N −2 ln 2 N , which is consistent with the theoretical bounds established in Theorems 3.1, 3.10, and 4.5.These experimental results validate the effectiveness of Richardson extrapolation.

Conclusion
This article explores using Richardson extrapolation to solve SPCDPs with non-local boundary conditions (2.1).We discretized the domain using a piecewise-uniform mesh and utilized the upwind finite difference scheme.To handle the non-local boundary conditions, we employed the trapezoidal rule for numerical integration.Our initial approximation of the scaled derivatives (εξ ′ ) and scaled discrete derivatives (εD − Z) resulted in an almost first-order convergence rate.In addition, we utilized the Richardson extrapolation method to greatly enhance accuracy, resulting in nearly a second-order convergence rate.Our analysis revealed an improvement in convergence rate from about O N −1 ln N to O N −2 ln 2 N with respect to ε, leading to more dependable and precise solutions with fewer errors at the nodes.We presented two instances that illustrated the highest pointwise errors and convergence rates for different values of ε and N .Finally, a comparison is made that demonstrates how post-processing techniques produce better, more accurate results.Future work could extend this method to handle problems with two parameters, PDE, and equations with a discontinuous source term with non-local boundary conditions.

Figure 4 :
Figure 4: The log-log plot of the maximum pointwise errors for Example 5.2

Table 1 :
Values of E N and S N for the solution of ξ and the scaled first derivatives εξ

Table 2 :
Values of E N and S N for the solution of ξ and the scaled first derivatives εξ