Singularly perturbed differential equations involve the highest-order derivative term being multiplied by a small perturbation parameter, ε. In the context of singularly perturbed problems, differential equations with two small parameters affecting the diffusion and convection terms are intriguing areas of research. The main purpose of this study is to obtain a robust higher-order numerical method to solve a two-parameter singularly perturbed boundary value problem. Find Υ such that

subject to the boundary conditions

where ε(0 < ε ≪ 1) and μ(0 < μ ≪ 1) are the small parameters that make the differential equation singularly perturbed. Let us assume the functions a(θ), b(θ) and f(θ) are sufficiently smooth and bounded to ensure the existence of a unique solution. Assume that there exist constants α, β and ζ independent of ε and μ such that for any θ = [0, 1] the conditions,

hold for some constant α, β and ζ. Assume δ≈minx∈Ω̄b(θ)a(θ). When μ = 1, Equation 1 is reduced to the singularly perturbed convection-diffusion boundary value problem. When μ = 0, Equation 1 reduces to a well-known singularly perturbed reaction-diffusion boundary value problem. The mathematical model for the kind of Equation 1 often arises in chemical reactor theory [1], transport phenomena in chemistry and biology [2], and lubrication theory [3]. The nature of the two-parameter problem was asymptotically examined by O'Malley [4], where the ratio of μ to ε has significant role in the solution. For this problem, two boundary layers occur at θ = 0 and θ = 1. Because of the presence of these layers, some standard numerical methods in Lodhi et al. [5], Kambampati et al. [6], Pandit and Kumar [7], Khandelwal and Khan [8] are applied to a uniform mesh, which provides an oscillatory numerical solution. Consequently, considerable attention has been devoted to using non-uniform meshes to solve two-parameter singularly perturbed boundary value problems, as discussed in works by O'Riordan and Pickett [9], Roos and Uzelac [10], Kadalbajoo and Yadaw [11], Kadalbajoo and Jha [12], Brdar and Zarin [13], Luo et al. [14], Brdar and Zarin [15, 16], Padmaja et al. [17], Andisso and Duressa [18], Linß and Roos [19], Cheng [20], Zhang and Lv [21], Valarmathi and Ramanujam [22], Patidar [23] to solve two-parameter singularly perturbed boundary value problems. The majority of the previously developed methods to solve the problem at hand are less accurate and of lower order. Inspired by this, the goal of this study is to offer a high order and more accurate fitted operator method with the help of post-processing technique to solve the considered problem. Because of the presence of boundary layer in the solution of Equation 1, devising a higher-order convergent numerical method is a big challenge. In this study, we apply the well-known post-processing technique of Richardson extrapolation method to obtain fourth order uniformly convergent numerical solution of Equation 1. The present approach yields a more accurate solution in terms of maximum absolute errors than previous approaches found in the literature.

Some robust numerical methods for one-parameter and two-parameter singularly perturbed problems in studies by Hassen and Duressa [24], Mohye et al. [25], Tesfaye et al. [26], Cheru et al. [27], Daba et al. [28], Gupta et al. [29], singularly perturbed turning point problem in Gupta et al. [30], time-fractional singularly perturbed convection-diffusion problem [31], singularly perturbed problems with spatio-time delays the study by Ejere et al. [32] and singularly perturbed problems Burger-Huxley problem in the study by Daba and Duressa [33] and Derzie et al. [34]. Some methods are employed studies by Li et al. [35], El Ahmadi et al. [36], Nachaoui [37], Mardanov et al. [38] to solve different types of differential equations.

The remaining part of the article is arranged as follows: In Section 2, we have provided a brief description of the present method for the numerical solution of Equations 1, 2. The convergence analysis of the method is presented in Section 3. Section 4 presents the numerical results, and comparisons are made with other existing methods. Finally, the conclusion is provided at the end of the article in Section 5.

In this section, we provide a priori bounds for the solution and its corresponding derivatives. The governing problem in Equations 1, 2 exhibits twin boundary layers with different layer widths depending on the relation between the values of ε and μ. If αμ² ≤ δε, then the reduced problem corresponding to Equations 1, 2 is given by

Thus, the boundary layers are expected near θ = 0 and θ = 1, with a width of O(ε) if Υ₀(0) ≠ γ₀ and Υ₀(1) ≠ γ₁. If αμ² ≥ δε, then the reduced problem corresponding to Equations 1, 2 is given by

Thus, the boundary layer of width O(ε/μ) is expected in the right neighborhood of θ = 0 if Υ_μ(0) ≠ γ₀ and the boundary layer of width O(μ) is expected in the left neighborhood of θ = 1 if Υ_μ(1) ≠ γ₁. The assumptions given for Equations 1, 2 ensure that the differential operator Π_{ε, μ} satisfies the following maximum principle.

Lemma 1. Let ϖ(θ) be a smooth function such that ϖ(0) ≥ 0 and ϖ(1) ≥ 0. If Π_ε,μ ϖ(θ) ≥ 0, for θ ∈ Ω, then ϖ(θ) ≥ 0, for θ∈Ω̄.

Proof. Assume θ^* be such that ϖ(θ*)=minθ∈Ω̄ϖ(θ) and ϖ(θ^*) < 0. Then, it is obvious that θ^* ∉ {0, 1}. Hence, ϖ′(θ^*) = 0 and ϖ^′′(θ^*) ≥ 0. Now, Πε,μϖ(θ*)≡-εϖ′′(θ*)+μa(θ)ϖ′(θ*)+b(θ)ϖ(θ*)<0, which is contradiction. Hence, we conclude that ϖ(θ) ≥ 0 for all θ ∈ [0, 1].     □

The following lemma proves the stability estimate to obtain a unique solution.

Lemma 2. On θ∈Ω̄, the solution Υ(θ) to the problem in Equations 1, 2 satisfy the bound

Proof. Defining two functions ζ^±(θ) such that

It is straightforward that ζ^±(0) ≥ 0 and ζ^±(1) ≥ 0. Now, for θ ∈ Ω,

Therefore, the desired result follows by applying Lemma (1).     □

The solution to the reduced problem μa(θ)Υ′(θ)+b(θ)Υ(θ) = f(θ), in general, does not satisfy both the boundary conditions; therefore, there exist boundary layers at both the boundaries, θ = 0 and θ = 1 [10]. To describe these boundary layers, the characteristic equation for the homogeneous part of Equation 1 with constant coefficients that are the minimum values of the corresponding variable coefficients is considered as

The characteristic equation in Equation 3 has two real solutions

The solution Ψ₀(θ) < 0 describes the boundary layer at θ = 0, whereas Ψ₁(θ) > 0 describes the boundary layer at θ = 1. To bound the solution and its derivatives, we define

The two real solutions Ψ₀(θ) < 0 and Ψ₁(θ) < 0 describing the boundary layers, respectively, at θ = 0 and θ = 1 are based on the two cases below.

Case 1. If μ2ε→0 as ε → 0, then

The governing Equation 1 has two boundary layers that behave like the reaction-diffusion case (μ≈0) with each of width O(ε) at θ = 1 and θ = 0. The complementary function of Equation 1 may be expressed as

where A₁ and B₁ are real constant numbers.

Case 2. If εμ2→0 as μ → 0, then

In this case, the governing Equation 1 has two boundary layers near θ = 0 and θ = 1 with different layer widths O(ε/μ) and O(μ), respectively. Now, the complementary function of Equation 1 can be given as

where A₂ and B₂ are real constant numbers. Note that most numerical methods give an accurate numerical solution for case 1, since μ≈0 behaves like a reaction-diffusion problem. For case 2, it is challenging to produce an accurate numerical solution. Therefore, in this study, we focus on case 2.

Theorem 1. For any 0 < p < 1, we have up to a certain order q that depends on the smoothness of the data. If a,b,f∈Ck(Ω̄), then the solution u(θ) satisfies

Proof. The details of the proof is well established in studies by Roos and Uzelac [10].     □

Let N be a positive integer and [0,1] be the closed domain, where N is the subinterval such that 0 = θ₀ < θ₁ < ⋯ < θ_N = 1 and θi=iℓ, ℓ=1N i = 0, 1, ⋯ , N. Using the notation Υ_i as a numerical approximation to the analytical solution Υ(θ_i) and second-order central finite difference approximations for the second- and first derivatives, we have the discrete problem as

with the discrete boundary conditions

where TE=ℓ2(112Υi(4)+16Υi(3))≈O(h2). In order to regulate the solution behavior of the singular perturbation parameter ε, we have introduced the fitting factor η on the homogeneous part of Equation 4.

The discrete problem in Equation 6 can be written as the three-term recurrence relation of the form

with the discrete boundary conditions in Equation 5 and where the coefficients are given by

The developed method is considered as an exponentially fitted operator finite difference method to solve the problem in Equations 1, 2. The coefficients L_i, M_i, and R_i are given to satisfy the conditions |L_i| > 0, |M_i| > 0, |R_i| > 0 and |M_i| ≥ |L_i|+|R_i|. These conditions guarantee that the linear system is diagonally dominant and can be solved by a tri-diagonal solver, which is the Thomas algorithm.

In this section, we prove the stability and convergence analysis of the discrete problem. First, we want to prove the discrete comparison principle for the discrete scheme in Equation 12.

Theorem 2. Assume Πε,μN be discrete operator and Θ_i be comparison function such that Πε,μNΥi≤Πε,μNΘi,∀i=1,2,⋯ ,N-1. If Υ₀ ≤ Θ₀ and Υ_N ≤ Θ_N, then Υ_i ≤ Θ_i ∀i = 0, ⋯ , N.

Proof. The matrix associated with operator Πε,μN is of size (N − 1) × (N − 1) and satisfies the property of M-matrix. That is, the inverse matrix exists, and it is nonnegative. See the detailed proof in Kellogg and Tsan [40]. This guarantees the existence and uniqueness of the discrete solution.     □

Lemma 3. Let Υ_i be the discrete solution. Then, we have the following bound

Proof. Let Ξ=β-1max∀i∈[0,1]|Πε,μNΥi|+max{|γ0|,|γ1|} and define the two barrier functions Ψi± by Ψi±=Ξ±Υi. At the boundary points, we have Ψ0±=Ξ±Υ0=Ξ±γ0≥0, and ΨN±=Ξ±ΥN=Ξ±γN≥0. On the discretized domain 1 ≤ i ≤ N − 1, we have

where b_i ≥ β > 0 and from Theorem (2), we get Ψi±≥0, for θi∈Ω̄.     □

We use the truncation error given in Equation 5 to show the convergence analysis of the present method as follows:

Theorem 3. Let Υ(θ_i) be the continuous solution and Υ_i be the discrete solution. Then, the error bound satisfies

where C is a constant independent of ε, μ and the mesh lengths h.

The above theorem shows that the present method is second-order convergent, independent of the parameters ε and μ. Next, we develop the post-processing technique to improve the accuracy of the present method and order of convergence.

In this section, we undertake computerized calculations to validate the efficacy of the proposed method against the theoretical results described in previous sections.

Example 1. Consider variable coefficient two parameter singularly perturbed problem

Example 2. Consider variable coefficient two parameter singularly perturbed problem

Since the exact solutions for each example are not available, the double mesh principle was employed to compute the maximum absolute errors for each (ε, μ):

where ΥiN is the numerical solution with N mesh points and Υi2N is the numerical solution at the finer mesh with 2N mesh points.

Example 3. Consider constant coefficient two parameter singularly perturbed problem

for which the analytical solution is given by

where

The maximum absolute errors for each (ε, μ) may be determined using the following formula, as the exact solution for Example (3) is known.

where ΥiN is the numerical solution with N mesh points and ΥN(θi) is the analytical solution. The (ε, μ)-maximum errors for all the Examples are calculated using the following formula

Furthermore, we compute the numerical rate of convergence before and after post-processing technique with the following formulas, respectively:

The (ε, μ)−maximum rates of convergence before and after post-processing techniques were calculated using the following formulas, respectively

Tables 1–3 show the calculated maximum errors eε,μN and the parameter-uniform errors e^N for Examples (1)–(3), respectively. These findings demonstrate that the current approach provides parameter-uniform convergence for both the before and after post-processing technique. Figures 1–3 display the plots of the numerical simulations for Examples (1)–(3). Figures 1–3 illustrate the plots of the numerical solution profile for Examples (1)–(3) for fixed μ and varying ε. From these figures, we observe that for fixed μ as ε → 0, strong layers are formed. Figures 4–6, respectively, show the log–log scale plots of the maximum errors for Examples (1), (2), and (3). The numerical findings show that the current higher-order fitted operator finite difference technique provides a numerical solution with more accuracy. The application of post-processing approach improves the accuracy of the numerical solution and speeds up the rate of convergence, as demonstrated by the numerical findings in all of the Tables.

A higher-order exponentially fitted operator finite difference method for two parameter singularly perturbed boundary value problems is presented in this study. The stability and uniform convergence of the current method are well established, ensuring second-order convergence. The post-processing technique is then applied to enhance the convergence order of the method and improve accuracy in terms of maximum errors. Theoretically, we have proven that the post-processing technique provides fourth-order parameter-uniform convergence. Three numerical examples are computed for various perturbation parameter values in order to verify the applicability of the current method. The present method can be applied to singularly perturbed parabolic problem with or without delay.