Preprint
Article

This version is not peer-reviewed.

Efficient Layer-Resolving Fitted Mesh Finite Difference Approach for Solving a System of n Two-Parameter Singularly Perturbed Convection-Diffusion Delay Differential Equations

A peer-reviewed article of this preprint also exists.

Submitted:

14 February 2025

Posted:

20 February 2025

You are already at the latest version

Abstract

This paper presents a robust fitted mesh finite difference method for solving a system of n singularly perturbed two parameter convection-reaction-diffusion delay differential equations defined on the interval [0,2]. Leveraging a piecewise uniform Shishkin mesh, the method adeptly captures the solution’s behavior arising from delay term and small perturbation parameters. The proposed numerical scheme is rigorously analyzed and proven to be parameter-robust, exhibiting nearly first-order convergence. Numerical Illustration is included to validate the method’s efficiency and to confirm the theoretical predictions.

Keywords: 
;  ;  ;  ;  ;  

1. Introduction

Singularly perturbed differential equations (SPDEs) are integral to many fields, including fluid dynamics, chemical reactor theory, population dynamics and control systems [1,2]. Within this class, singularly perturbed delay differential equations (SPDDEs) present additional complexities due to boundary and interior layers that arise from small perturbation parameters and delay terms, making numerical solutions challenging. Established techniques like fitted mesh [3] and fitted operator methods [4] have provided accurate solutions for SPDEs, with notable work by Cen [5], who used a hybrid scheme with a Shishkin mesh to achieve near-second-order convergence. Similarly, Gracia et al. [6] developed a monotone method for SPDEs with two parameters influencing convection and diffusion. SPDDEs typically involve boundary value problems influenced by two small parameters, μ and ϵ i ( i = 1 , 2 , , n ) , whose interactions generate complex layer behavior governed by the ratio μ 2 / ϵ i . This work aims to construct a parameter-robust numerical method for system of n SPDDEs as both parameters approach zero, addressing boundary and interior layers accurately regardless of perturbation values. Stability is established, and bounds for the solution’s derivatives are derived, supporting the convergence of the fitted mesh finite difference approach, which attains nearly first-order accuracy. Several studies have addressed singular perturbation problems, emphasizing their asymptotic behaviour, the development of parameter uniform methods and delay differential equations to achieve robust and accurate solutions across varying the perturbation parameters [7,8,9]. The novelty of this research lies in addressing a system of n SPDDEs involving two small parameters in a convection-diffusion context, contrasting with prior studies that either considered single delay equations [10], system of two equations without delay terms [11] and system of two equations with delay terms [12]. This paper considers interacting variables affected by both delay and two distinct small perturbation parameters, ϵ i ( i = 1 , 2 , , n ) and μ introducing challenges in boundary and interior layer formation. To overcome these challenges, this paper presents advanced mesh techniques and robust numerical schemes that effectively resolve layers and achieve parameter uniform convergence, significantly enhancing the numerical analysis of SPDDEs.

2. Formulation of the Problem

The system of singularly perturbed delay differential equations involving two small parameters is under consideration,
E u ( ϰ ) + μ A ( ϰ ) u ( ϰ ) B ( ϰ ) u ( ϰ ) + D ( ϰ ) u ( ϰ 1 ) = f ( ϰ ) for all ϰ Ω = ( 0 , 2 ) ,
u = φ on [ 1 , 0 ] , u ( 2 ) = ι .
Here, u = u 1 u 2 u n , f = f 1 f 2 f n ,   E = ϵ 1 0 0 0 ϵ 2 0 0 0 ϵ n , A ( ϰ ) = a 1 ( ϰ ) 0 0 0 a 2 ( ϰ ) 0 0 0 a n ( ϰ ) ,   B ( ϰ ) = b 11 ( ϰ ) b 12 ( ϰ ) b 1 n ( ϰ ) b 21 ( ϰ ) b 22 ( ϰ ) b 2 n ( ϰ ) b n 1 ( ϰ ) b n 2 ( ϰ ) b n n ( ϰ ) ,   D ( ϰ ) = d 1 ( ϰ ) 0 0 0 d 2 ( ϰ ) 0 0 0 d n ( ϰ ) , where ϵ i , i = 1 , 2 , , n are small parameters satisfying 0 < ϵ 1 < ϵ 2 < < ϵ n 1 and μ is another small parameter with 0 < μ 1 . The coefficient functions a i ( ϰ ) , b i j ( ϰ ) , d i ( ϰ ) and f i ( ϰ ) are all sufficiently smooth throughout the domain Ω ¯ = [ 0 , 2 ] and a i ( ϰ ) α > 0 , b i i ( ϰ ) + j = 1 n b i j ( ϰ ) β > 0 , d i ( ϰ ) > 0 ,   b i i ( ϰ ) + j = 1 n b i j ( ϰ ) d i ( ϰ ) κ > 0 , for i , j = 1 , 2 , , n and i j . The φ ( x ) is sufficiently smooth over the interval [ 1 , 0 ] . The value of γ is determined as γ = min ϰ Ω ¯ j = 1 n b i j ( ϰ ) d i ( ϰ ) a i ( ϰ ) for i = 1 , 2 , , n and i j . When μ = 0 , the above problem without delay, has been considered in [13]. The problem demonstrates boundary layers influenced by both ϵ i and μ , in particular, the layers are influenced by the ratio of μ 2 ϵ i . If μ 2 ϵ i γ α , 1 i n , the reduced problem can be expressed as
B ( ϰ ) u ( ϰ ) + D ( ϰ ) u ( ϰ 1 ) = f ( ϰ ) , ϰ ( 0 , 2 ] ,
with boundary conditions u ( ϰ ) = φ ( ϰ ) on [ 1 , 0 ] . This predicts that there will be a boundary layer of width O ( ϵ i ) near ϰ = 0 and an interior layer at ϰ = 1 , arising from the unit delay component, under the assumption that u ( 0 ) φ ( 0 ) . Additionally, a similar boundary layer of width O ( ϵ i ) is anticipated near ϰ = 2 , along with an interior layer at ϰ = 1 , if u ( 2 ) ι . If μ 2 ϵ j γ α , 1 < i < j n , the reduced problem is
μ A ( ϰ ) u ( ϰ ) B ( ϰ ) u ( ϰ ) + D ( ϰ ) u ( ϰ 1 ) = f ( ϰ ) , ϰ ( 0 , 2 ) ,
with boundary conditions u ( ϰ ) = φ ( ϰ ) on [ 1 , 0 ] , u ( 2 ) = ι . This problem remains to exhibit singular perturbation behavior due to the parameter μ . It is expected that a boundary layer with a width O ( μ ) is anticipated near ϰ = 2 , and an interior layer at ϰ = 1 (due to the delay), unless the solution at ϰ = 2 differs from ι . Additionally, a boundary layer of width O ( ϵ i μ ) is anticipated near ϰ = 0 , with an interior layer at ϰ = 1 , if u ( 0 ) φ ( 0 ) . Numerical experiments indicate that the interior right layer weakens considerably when ϵ i μ 2 . Consider the problem
L 1 z ( ϰ ) = E z ( ϰ ) + μ A ( ϰ ) z ( ϰ ) B ( ϰ ) z ( ϰ ) = g ( ϰ ) , ϰ Ω 1 ,
L 2 z ( ϰ ) = E z ( ϰ ) + μ A ( ϰ ) z ( ϰ ) B ( ϰ ) z ( ϰ ) + D ( ϰ ) z ( ϰ 1 ) = f ( ϰ ) , ϰ Ω 2 ,
with boundary conditions z ( 0 ) = φ ( 0 ) , z ( 1 ) = z ( 1 + ) , z ( 1 ) = z ( 1 + ) , z ( 2 ) = ι , where g ( ϰ ) = f ( ϰ ) D ( ϰ ) φ ( ϰ 1 ) , for ϰ Ω 1 , Ω 1 = ( 0 , 1 ) and Ω 2 = ( 1 , 2 ) .

3. Analytical Results

This section presents a minimum principle, establish a stability result for the problem described by Equation (1) and provide its estimates for the derivatives of the solution.
Lemma 3.1.
Let ψ = ( ψ 1 , ψ 2 , , ψ n ) T be such that ψ ( 0 ) 0 , ψ ( 2 ) 0 , L 1 ψ 0 on ( 0 , 1 ) and L 2 ψ 0 on ( 1 , 2 ) , ψ ( 1 ) = 0 , ψ ( 1 ) 0 , then ψ 0 on [ 0 , 2 ] .
Proof.
Assume ϰ * and s * are such that ψ s * ( ϰ * ) = min ϰ Ω , s = 1 , 2 , , n ψ s ( ϰ ) . Suppose ψ s * ( ϰ * ) < 0 . Then, ϰ * cannot be at the boundaries 0 or 2. At ϰ * , the first derivative of ψ s * , denoted as ψ s * ( ϰ * ) = 0 and the second derivative ψ s * ( ϰ * ) 0 .
Claim (i): ϰ * ( 0 , 1 ) . If ϰ * ( 0 , 1 ) , then
( L 1 ψ ) s * ( ϰ * ) = ϵ s * ψ s * ( ϰ * ) + μ a s * ( ϰ * ) ψ s * ( ϰ * ) j = 1 n b s * j ( ϰ * ) ψ j ( ϰ * ) > 0 ,
which contradicts the assumption that L 1 ψ 0 on ( 0 , 1 ) . Thus, ϰ * ( 0 , 1 ) .
Claim (ii): ϰ * [ 1 , 2 ) . If ϰ * [ 1 , 2 ) , then
( L 2 ψ ) s * ( ϰ * ) = ϵ s * ψ s * ( ϰ * ) + μ a s * ( ϰ * ) ψ s * ( ϰ * ) j = 1 n [ b s * j ( ϰ * ) ] ψ j ( ϰ * ) d s * ( ϰ * ) ψ s * ( ϰ * 1 ) > 0 ,
which contradicts the assumption that L 2 ψ 0 on [ 1 , 2 ) . Hence, ϰ * [ 1 , 2 ) . When ϰ * = 1 , the differentiability of ψ s * at ϰ = 1 . If ψ s * ( 1 ) does not exist, then ψ s * ( 1 ) = ψ s * ( 1 + ) ψ s * ( 1 ) > 0 , which contradicts the condition ψ s * ( 1 ) 0 . If ψ s * is differentiable at ϰ = 1 , then μ a s * ( 1 ) ψ s * ( 1 ) j = 1 n b s * j ( 1 ) ψ j ( 1 ) > 0 . Since all entries of A ( ϰ ) , B ( ϰ ) and ψ j ( ϰ ) are continuous on [ 0 , 2 ] , there exists an interval [ 1 h , 1 ) where
μ a s * ( ϰ ) ψ s * ( ϰ ) j = 1 n b s * j ( ϰ ) ψ j ( ϰ ) > 0 .
Next, the second derivative ψ s * ( ϰ ) is examined, if ψ s * ( p ) 0 for any p [ 1 h , 1 ) , then ( L 1 ψ ) s * ( p ) 0 , which leads to a contradiction. Therefore, the assumption ψ s * ( ϰ ) < 0 on the interval [ 1 h , 1 ) . This indicates that ψ s * ( ϰ ) is strictly decreasing over [ 1 h , 1 ) . Given that ψ s * ( 1 ) = 0 and ψ s * is continuous on ( 0 , 2 ) , it follows that ψ s * ( ϰ ) > 0 on [ 1 h , 1 ) . As a result, the continuous function ψ s * ( ϰ ) cannot attain a minimum at ϰ = 1 , which contradicts the assumption that ϰ * = 1 . Thus, ψ 0 on [ 0 , 2 ] . The proof of the lemma is complete.□
Lemma 3.2.(Stability Result)
Let ψ C 2 ( Ω ¯ ) . For ϰ Ω ¯ ,
| ψ i ( ϰ ) | max | ψ ( 0 ) | , | ψ ( 2 ) | , 1 κ L 1 ψ Ω 1 , 1 κ L 2 ψ Ω 2 .
Proof.
Define
M = max | ψ ( 0 ) | , | ψ ( 2 ) | , 1 κ L 1 ψ Ω 1 , 1 κ L 2 ψ Ω 2 .
Consider the functions θ ± ( ϰ ) = M e ± ψ ( ϰ ) , where e = ( 1 , 1 , , 1 ) T . Clearly, θ ± ( 0 ) 0 , θ ± ( 2 ) 0 , L 1 θ ± ( ϰ ) 0 and L 2 θ ± ( ϰ ) 0 for all ϰ Ω . Hence, by Lemma 3.1, proves that | ψ i ( ϰ ) | M , which yields the required result. The proof of the lemma is complete.□
Theorem 3.1.
Let u be the solution of (1) and then, its derivatives satisfy the following bounds on Ω,
| u i ( k ) ( ϰ ) | C ( ϵ i ) k 1 + μ ϵ i k max { u , f } ,
u i ( 3 ) ( ϰ ) C ( ϵ i ) 3 1 + μ ϵ i 3 max u , f , f ,
| u i ( 4 ) ( ϰ ) | C ( ϵ i ) 4 1 + μ ϵ i 4 max u , f , f , f
where the constant C is independent of ϵ i and μ , i = 1 , 2 , , n and k = 1 , 2 .
Proof.
The proof follows the methodology outlined in Lemma 2.2 of [8]. For any ϰ ( 0 , 1 ) , a neighborhood N p = ( p , p + ϵ i ) such that ϰ N p and N p ( 0 , 1 ) . According to the mean value theorem, there exists a y N p satisfying
u i ( y ) = u i ( p + ϵ i ) u i ( p ) ϵ i | u i ( y ) | 2 u i ϵ i .
Now, it follows that u i ( ϰ ) = u i ( y ) + y ϰ u i ( ) d . Thus,
| u i ( ϰ ) | C ϵ i 1 + μ ϵ i max { u , f } .
The bounds for u i is obtained from Equation (1). Similarly, the bounds of u i ( 3 ) and u i ( 4 ) can be established for higher-order derivatives through analogous manipulations. The proof of the theorem is complete. □

4. Shishkin Decomposition of the Solution

For each of the cases α μ 2 γ ϵ i and α μ 2 γ ϵ j , u is expressed by u = v + w L + w R , where
v ( ϰ ) = r ( ϰ ) , for ϰ [ 0 , 1 ) s ( ϰ ) , for ϰ [ 1 , 2 ] ,
w L ( ϰ ) = w L 1 ( ϰ ) = w 1 L 1 ( ϰ ) w 2 L 1 ( ϰ ) w n L 1 ( ϰ ) , for ϰ [ 0 , 1 ) w L 2 ( ϰ ) = w 1 L 2 ( ϰ ) w 2 L 2 ( ϰ ) w n L 2 ( ϰ ) , for ϰ [ 1 , 2 ] , w R ( ϰ ) = w R 1 ( ϰ ) = w 1 R 1 ( ϰ ) w 2 R 1 ( ϰ ) w n R 1 ( ϰ ) , for ϰ [ 0 , 1 ) w R 2 ( ϰ ) = w 1 R 2 ( ϰ ) w 2 R 2 ( ϰ ) w n R 2 ( ϰ ) , for ϰ [ 1 , 2 ]
Case (i): α μ 2 γ ϵ i
In this case, for 1 i n ,
L 1 r ( ϰ ) = g ( ϰ ) , for ϰ ( 0 , 1 ) , r ( 0 ) and r ( 1 ) are selected ,
L 2 s ( ϰ ) = f ( ϰ ) , for ϰ ( 1 , 2 ) , s ( 1 ) and s ( 2 ) are selected ,
s ( ϰ ) = r ( ϰ ) on [ 0 , 1 ) ,
L 1 w L 1 ( ϰ ) = 0 , for ϰ ( 0 , 1 ) , w L 1 ( 0 ) = u ( 0 ) v ( 0 ) c 1 ( ϵ i , μ ) , w L 1 ( 1 ) = 0 ,
L 2 w L 2 ( ϰ ) = 0 , for ϰ ( 1 , 2 ) , w L 2 ( 1 ) = k 1 ( ϵ i , μ ) c 2 ( ϵ i , μ ) , w L 2 ( 2 ) = 0 ,
w L 2 ( ϰ ) = w L 1 ( ϰ ) , for ϰ [ 0 , 1 ) ,
L 1 w R 1 ( ϰ ) = 0 , for ϰ ( 0 , 1 ) , w R 1 ( 0 ) = c 1 ( ϵ i , μ ) , w R 1 ( 1 ) = k 2 ( ϵ i , μ ) ,
L 2 w R 2 ( ϰ ) = 0 , for ϰ ( 1 , 2 ) , w R 2 ( 1 ) = c 2 ( ϵ i , μ ) , w R 2 ( 2 ) = u ( 2 ) v ( 2 ) ,
w R 2 ( ϰ ) = w R 1 ( ϰ ) , for ϰ [ 0 , 1 ) .
Case (ii): α μ 2 γ ϵ j
In this case, for 1 < i < j n ,
L 1 r ( ϰ ) = g ( ϰ ) , for ϰ ( 0 , 1 ) , r ( 0 ) and r ( 1 ) are selected ,
L 2 s ( ϰ ) = f ( ϰ ) , for ϰ ( 1 , 2 ) , s ( 1 ) and s ( 2 ) are selected ,
s ( ϰ ) = r ( ϰ ) on [ 0 , 1 ) ,
L 1 w L 1 = 0 , for ϰ ( 0 , 1 ) , w L 1 ( 0 ) = u ( 0 ) v ( 0 ) c 1 ( ϵ i , μ ) , w L 1 ( 1 ) = 0 ,
L 2 w L 2 = 0 , for ϰ ( 1 , 2 ) , w L 2 ( 1 ) = k 1 ( ϵ i , μ ) c 2 ( ϵ i , μ ) , w L 2 ( 2 ) = 0 ,
w L 2 ( ϰ ) = w L 1 ( ϰ ) , for ϰ [ 0 , 1 ) ,
L 1 w R 1 = 0 , for ϰ ( 0 , 1 ) , w R 1 ( 0 ) = c 1 ( ϵ i , μ ) , w R 1 ( 1 ) = k 2 ( ϵ i , μ ) ,
L 2 w R 2 = 0 , for ϰ ( 1 , 2 ) , w R 2 ( 1 ) = c 2 ( ϵ i , μ ) , w R 2 ( 2 ) = u ( 2 ) v ( 2 ) ,
w R 2 ( ϰ ) = w R 1 ( ϰ ) , for ϰ [ 0 , 1 ) .
To ensure the jump conditions at ϰ = 1 in Equations (10) and (19) are satisfied, the constants k 1 ( ϵ i , μ ) and k 2 ( ϵ i , μ ) must be selected appropriately. Additionally, the constants c 1 ( ϵ i , μ ) and c 2 ( ϵ i , μ ) should be determined independently for the cases α μ 2 γ ϵ i and α μ 2 γ ϵ j , ensuring they meet the bounds required for the singular component. Given that u ( 0 ) and u ( 1 ) are bounded by constants that do not depend on ϵ i and μ , even though c 1 , c 2 , k 1 and k 2 are functions of ϵ i and μ , the magnitudes | c 1 | , | c 2 | , | k 1 | and | k 2 | are constants independent of ϵ i and μ .

5. Bounds on the Regular Component and Its Derivatives

To establish the result, by estimating bounds for the smooth components and their derivatives on the interval [ 0 , 1 ] and then use these bounds to extend the estimates to the interval [ 1 , 2 ] . Specifically, by decomposing each component with respect to ϵ n , then apply ϵ n 1 to the first n 1 components, followed by ϵ n 2 for the first n 2 components, and so on. This step-by-step decomposition approach is as follows for both cases.
Case (i): α μ 2 γ ϵ i
Establishing the bounds of the regular components r and s , it is broken down as
r = y n + ϵ n z n + ϵ n 2 q n + ϵ n 3 p n , s = n + ϵ n g n + ϵ n 2 n + ϵ n 3 t n
where y n = ( y n 1 , y n 2 , . . . , y n n ) T , n = ( n 1 , n 2 , . . . , n n ) T is the solution of
B n ( ϰ ) y n ( ϰ ) = g ( ϰ ) , for ϰ [ 0 , 1 ] ,
B n ( ϰ ) n ( ϰ ) + D n ( ϰ ) y n ( ϰ 1 ) = f ( ϰ ) , for ϰ [ 1 , 2 ] ,
where z n = ( z n 1 , z n 2 , . . . , z n n ) T , g n = ( g n 1 , g n 2 , . . . , g n n ) T is the solution of
B n ( ϰ ) z n ( ϰ ) = ϵ n 1 E y n + μ ϵ n 1 A n y n , for ϰ [ 0 , 1 ] ,
B n ( ϰ ) g n ( ϰ ) D n ( ϰ ) z n ( ϰ 1 ) = ϵ n 1 E n + μ ϵ n 1 A n n , for ϰ [ 1 , 2 ] ,
where q n = ( q n 1 , q n 2 , . . . , q n n ) T , n = ( n 1 , n 2 , . . . , n n ) T is the solution of
B n ( ϰ ) q n ( ϰ ) = ϵ n 1 E z n + μ ϵ n 1 A n z n , for ϰ [ 0 , 1 ] ,
B n ( ϰ ) n ( ϰ ) D n ( ϰ ) q n ( ϰ 1 ) = ϵ n 1 E g n + μ ϵ n 1 A n g n , for ϰ [ 1 , 2 ] ,
where p n = ( p n 1 , p n 2 , . . . , p n n ) T , t n = ( t n 1 , t n 2 , . . . , t n n ) T is the solution of
L 1 p n ( ϰ ) = ϵ n 1 E q n + μ ϵ n 1 A n q n on ( 0 , 1 )
L 2 t n ( ϰ ) = ϵ n 1 E n + μ ϵ n 1 A n n on ( 1 , 2 ) , z n ( 2 ) = 0 and z n ( 0 ) remains to be chosen .
Since the matrix ϵ n 1 E has entries that are bounded, it follows that, for 0 k 3 ,
y n ( k ) C , n ( k ) C , z n ( k ) C , g n ( k ) C , q n ( k ) C , n ( k ) C .
Now using Theorem 3.1 and (34),(35), for the choice of p n n ( 0 ) = 0 , t n n ( 0 ) = 0 , then
| p n n ( k ) | C ϵ n k , | t n n ( k ) | C ϵ n k .
Then from (36) and (37) yields
| r n ( k ) | C ( 1 + ( ϵ n ) 3 k ) , | s n ( k ) | C ( 1 + ( ϵ n ) 3 k ) .
To establish the estimation the bounds r i ( k ) and s i ( k ) for 1 i n 1 and notations are defined for 1 l n , E l = ϵ 1 0 0 0 ϵ 2 0 0 0 ϵ l ,   A l = a 1 0 0 0 a 2 0 0 0 a l ,   B l = b 11 b 12 b 1 l b 21 b 22 b 2 l b l 1 b l 2 b l l ,   D l = d 1 0 0 0 d 2 0 0 0 d l ,   p ˜ l = ( p l 1 , p l 2 , , p l ( l 1 ) ) T , g ( l 1 ) = ( g ( l 1 ) 1 , g ( l 1 ) 2 , , g ( l 1 ) ( l 1 ) ) T , with g ( l 1 ) j = ϵ j ϵ l q l j + μ ϵ l A l q l j + b j l p l l . To proceed with the analysis, let us focus on the first n 1 equations of the system described by equations (34) and (35). It follows that
L ˜ 1 , ( n ) p ˜ n E n 1 p ˜ n ( ϰ ) + A n 1 ( ϰ ) p ˜ n ( ϰ ) B n 1 ( ϰ ) p ˜ n ( ϰ ) = g n 1 ( ϰ ) ,
L ˜ 2 , ( n ) p ˜ n E n 1 t ˜ n ( x ) + A n 1 ( x ) t ˜ n ( x ) B n 1 ( x ) t ˜ n ( x ) + D n 1 ( ϰ ) t ˜ n ( ϰ 1 ) = f n 1 ( x ) .
The next step involves decomposing p ˜ n similarly to equation above
p ˜ n = y n 1 + ϵ n 1 z n 1 + ϵ n 1 2 q n 1 + ϵ n 1 3 p n 1 .
Similarly proceeding like above, the problem associated with p n 1 , is similar as in (34) and (35). By applying the estimates, the bound on the solution is obtained for 0 k 3 , y n 1 ( k ) C ( 1 + ( ϵ n ) 1 k ) ,   z n 1 ( k ) C ( ϵ n ) k ,   q n 1 ( k ) C ( ϵ n ) k 1 . Then, by applying Theorem 3.1, and q n 1 ( k ) yields | p ( n 1 ) ( n 1 ) ( k ) | C ϵ n 3 ϵ n 1 k . Therefore, | r n 1 ( k ) | C ( 1 + ( ϵ n 1 ) 3 k ) , similarly for [ 1 , 2 ] , the bounds implies | s n 1 ( k ) | C ( 1 + ( ϵ n 1 ) 3 k ) . In an analogous way, singularly perturbed systems of l equations is derived, where l = n 2 , n 3 , , 2 , 1 ,
L ˜ 1 , ( l + 1 ) p ˜ l + 1 E l p ˜ l + 1 ( ϰ ) + A l ( ϰ ) p ˜ l + 1 ( ϰ ) B l ( ϰ ) p ˜ l + 1 ( ϰ ) = g l ( ϰ ) ,
L ˜ 2 , ( l + 1 ) p ˜ l + 1 E l t ˜ l + 1 ( x ) + A l ( x ) t ˜ l + 1 ( x ) B l ( x ) t ˜ l + 1 ( x ) + D l ( ϰ ) t ˜ l + 1 ( ϰ 1 ) = f l ( x ) .
Using similar decomposition, it can be found that | r l ( k ) | C ( 1 + ( ϵ l ) 3 k ) , similarly for [ 1 , 2 ] , it follows that | s l ( k ) | C ( 1 + ( ϵ l ) 3 k ) . Thus, satisfies the following bound for 1 i n and 0 k 3 ,
| r i ( k ) | C ( 1 + ( ϵ i ) 3 k ) , | s i ( k ) | C ( 1 + ( ϵ i ) 3 k ) .
Case (ii): α μ 2 γ ϵ j
Establishing the bounds of the regular components r and s , it is broken down as
r = y n + ϵ n z n + ϵ n 2 q n + ϵ n 3 p n , s = n + ϵ n g n + ϵ n 2 n + ϵ n 3 t n .
Furthermore, the maximum principle for a linear first-order operator is established and demonstrated within the framework of a terminal value problem. Define the operators
L 3 μ A ( ϰ ) u ( ϰ ) B ( ϰ ) u ( ϰ ) ,
L 4 μ A ( ϰ ) u ( ϰ ) B ( ϰ ) u ( ϰ ) + D ( ϰ ) u ( ϰ 1 ) .
Decompose y n , z n , q n , n , g n , n individually and similarly proceeding like case(i), the following bounds are established for all i in the range 1 i n and k in the range 0 k 3 ,
| r i ( k ) | C ( 1 + ϵ i 3 k μ k 3 ) , | s i ( k ) | C ( 1 + ϵ i 3 k μ k 3 ) .

6. Layer Functions

The functions for the layers are denoted by B i l 1 ( ϰ ) and B i r 1 ( ϰ ) , 1 i n are specified over the interval [ 0 , 1 ] ,
B i l 1 ( ϰ ) = e θ i ϰ , α μ 2 γ ϵ i e λ i ϰ , α μ 2 γ ϵ j , B i r 1 ( ϰ ) = e θ i ( 1 ϰ ) , α μ 2 γ ϵ i e κ ( 1 ϰ ) , α μ 2 γ ϵ j .
The layer functions B i l 2 ( ϰ ) , B i r 2 ( ϰ ) , 1 i n are specified over [ 1 , 2 ] ,
B i l 2 ( ϰ ) = e θ i ( ϰ 1 ) , α μ 2 γ ϵ i e λ i ( ϰ 1 ) , α μ 2 γ ϵ j , B i r 2 ( ϰ ) = e θ i ( 2 ϰ ) , α μ 2 γ ϵ i e κ ( 2 ϰ ) , α μ 2 γ ϵ j ,
where θ i = γ α ϵ i ; λ i = α μ ϵ i ; κ = γ 2 μ , for 1 i n . Following the Lemma 5 presented in [13], the points ϰ s 0 , 1 2 which satisfy the conditions, B i l p , B j l p , B i r p , B j r p , 1 i < j n , p = 1 , 2 and for the case α μ 2 γ ϵ i , can be proved.
B i l 1 ( ϰ i , j ( k ) ) ϵ i k = B j l 1 ( ϰ i , j ( k ) ) ϵ j k on [ 0 , 1 ] , B i l 2 ( ϰ i , j ( k ) 1 ) ϵ i k = B j l 2 ( ϰ i , j ( k ) 1 ) ϵ j k on [ 1 , 2 ] ,
B i r 1 ( 1 ϰ i , j ( k ) ) ϵ i k = B j r 1 ( 1 ϰ i , j ( k ) ) ϵ j k on [ 0 , 1 ] , B i r 2 ( 2 ϰ i , j ( k ) ) ϵ i k = B j r 2 ( 2 ϰ i , j ( k ) ) ϵ j k on [ 1 , 2 ] , 0 k 2 .
Similarly for the case α μ 2 γ ϵ j , it can be shown that there exist points ϰ i , j ( s ) , s = 1 , 2 , 3 in 0 , 1 2 such that
B i l p ( ϰ i , j ( s ) ) ϵ i s = B j l p ( ϰ i , j ( s ) ) ϵ j s .

7. Bounds on the Singular Component and Its Derivatives

Theorem 7.1.
Let w L , w R satisfy problems (13), (16) and (22), (25) for the cases α μ 2 γ ϵ i and α μ 2 γ ϵ j respectively. Consequently, the components of w L p and w R p , p=1,2 satisfy the following bounds on ( 0 , 1 ) . For the case α μ 2 γ ϵ i , 1 i n ,
| w i L p ( ϰ ) | C B n l p ( ϰ ) , | w i L p , ( 1 ) ( ϰ ) | C ϵ i 1 / 2 B i l p ( ϰ ) + ϵ n 1 / 2 B n l p ( ϰ ) , | w i L p , ( 2 ) ( ϰ ) | C k = i n ϵ k 1 / 2 B k l p ( ϰ ) , | w i L p , ( 3 ) ( ϰ ) | C ϵ i 1 k = 1 i 1 ϵ k 1 / 2 B k l p ( ϰ ) + k = i n ϵ k 1 / 2 B k l p ( ϰ ) .
For the case α μ 2 γ ϵ j , 1 < i < j n ,
| w i L p ( ϰ ) | C B n l p ( ϰ ) , | w i L p , ( 1 ) ( ϰ ) | C μ ϵ i 1 B i l p ( ϰ ) + ϵ n 1 B n l p ( ϰ ) , | w i L p , ( 2 ) ( ϰ ) | C μ 2 k = j n ϵ k 2 B k l p ( ϰ ) , | w i L p , ( 3 ) ( ϰ ) | C μ 3 ϵ i 1 k = 1 j 1 ϵ k 2 B k l p ( ϰ ) + k = j n ϵ k 3 B k l p ( ϰ ) .
Moreover, the components satisfy the following bounds of w R 1 . For the case α μ 2 γ ϵ i ,
| w i R p ( ϰ ) | C B n r p ( ϰ ) , | w i R p , ( 1 ) ( ϰ ) | C ϵ i 1 / 2 B i r p ( ϰ ) + ϵ n 1 / 2 B n r p ( ϰ ) , | w i R p , ( 2 ) ( ϰ ) | C k = i n ϵ k 1 / 2 B k r p ( ϰ ) , | w i R p , ( 3 ) ( ϰ ) | C ϵ i 1 k = 1 i 1 ϵ k 1 / 2 B k r p ( ϰ ) + k = i n ϵ k 1 / 2 B k r p ( ϰ ) .
For the case α μ 2 γ ϵ j ,
| w i R p ( ϰ ) | C B n r p ( ϰ ) , | w i R p , ( k ) ( ϰ ) | C μ k B i r p ( ϰ ) , k = 1 , 2 , 3 .
Proof. For the case α μ 2 γ ϵ j , defining the barrier functions ψ ± = ψ 1 ± , ψ 2 ± , , ψ n ± , where ψ i ± = C B n l 1 ± w i L 1 for i = 1 , 2 , , n . It is evident that ψ i ± ( 0 ) 0 and ψ i ± ( 1 ) 0 . Additionally, ( L 1 ψ ± ) i ( ϰ ) 0 for all x in the interval ( 0 , 1 ) . Therefore, it can be demonstrated that | w i L 1 ( ϰ ) | C B n l 1 ( ϰ ) . Considering the equation of w i L 1 from (22),
ϵ i w i L 1 , ( ϰ ) + μ a i ( ϰ ) w i L 1 , ( ϰ ) + j = 1 n b i j ( ϰ ) w j L 1 ( ϰ ) = 0 .
This can also be written as, w i L 1 , ( ϰ ) + μ ϵ i a i ( ϰ ) w i L 1 , ( ϰ ) = 1 ϵ i j = 1 n b i j ( ϰ ) w j L 1 ( ϰ ) 1 ϵ i h i ( ϰ ) , where h i ( ϰ ) = j = 1 n b i j ( ϰ ) w j L 1 ( ϰ ) . Now, taking A i ( 0 ) = a i ,
w i L 1 , ( ϰ ) = w i L 1 , ( 0 ) e μ ϵ i ( A i ( ϰ ) ) + ϵ i 1 0 x h i ( s ) e μ ϵ j ( A i ( s ) A i ( ϰ ) ) d s ,
where A i ( ϰ ) is the indefinite integral of a i ( ϰ ) . Using the bounds on u , it is established that | w i L 1 , ( 0 ) | C ϵ i 1 . Using the inequality e μ ϵ i ( A i ( s ) A i ( ϰ ) ) e μ ϵ i β ( ϰ s ) and using integration by parts, it follows from the above that
| w i L 1 , ( ϰ ) | C ϵ i 1 B i l 1 ( ϰ ) + ϵ n 1 B n l 1 ( ϰ ) .
Using a similar argument, it can be | w i L 1 , ( ϰ ) | C k = j n ϵ k 2 B k l 1 ( ϰ ) . Differentiating the equation and using a similar procedure as above, it can be shown that
| w i L 1 , ( ϰ ) | C ϵ i 1 k = 1 j 1 ϵ k 2 B k l 1 ( ϰ ) + k = j n ϵ k 3 B k l 1 ( ϰ ) .
It has been established that | w i L 1 , ( ϰ ) | C ϵ i 1 B i l 1 ( ϰ ) + ϵ n 1 B n l 1 ( ϰ ) . Consequently,
| L 1 w i L 1 , ( ϰ ) | C μ ϵ i B i l 1 ( ϰ ) + μ ϵ n B n l 1 ( ϰ ) .
By introducing the barrier function ψ ± ( ϰ ) = C μ ϵ i B i l 1 ( ϰ ) + μ ϵ n B n l 1 ( ϰ ) ± w i L 1 , ( ϰ ) , it can be demonstrated that ψ ± 0 on [ 0 , 1 ] and L 1 ψ ± ( ϰ ) 0 on [ 0 , 1 ] , which implies
| w i L 1 , ( ϰ ) | C μ ϵ i B i l 1 ( ϰ ) + μ ϵ n B n l 1 ( ϰ ) .
By introducing another barrier function ϕ ± ( ϰ ) = C μ 2 k = 1 j 1 ϵ k 2 B k l 1 ( ϰ ) + μ 2 k = j n ϵ k 2 B k l 1 ( ϰ ) ± w 1 L 1 , ( ϰ ) , It can be derived that | w i L 1 , ( ϰ ) | C μ 2 k = j n ϵ k 2 B k l 1 ( ϰ ) . Differentiating the equations of w i L 1 once and applying the bounds of w i L 1 , and w i L 1 , yields
| w i L 1 , ( ϰ ) | C μ 3 ϵ i 1 k = 1 j 1 ϵ k 2 B k l 1 ( ϰ ) + k = j n ϵ k 3 B k l 1 ( ϰ ) .
Similarly, this leads to analogous results for w L 2 on [ 1 , 2 ] . Next, the bounds on w L 1 for the case α μ 2 γ ϵ i are derived. The bounds on w i L 1 ( ϰ ) can be derived by defining the barrier functions, ψ i ± ( ϰ ) = C B n l 1 ( ϰ ) ± w i L 1 ( ϰ ) , i = 1 , 2 , , n . To bound w i L 1 , ( ϰ ) , the argument proceeds by examining the equation of w i L 1 in (13) and proceed by applying the result from Theorem 3.1, | w i L 1 , ( ϰ ) | C ϵ i 1 / 2 B i l 1 ( ϰ ) . To improve the above bound on | w i L 1 ( ϰ ) | , this process continues by differentiating w i L 1 from (13) once leads to | ( L 1 w L 1 , ) i ( ϰ ) | C ϵ i 1 / 2 B i l 1 ( x ) . To establish the necessary bounds, the barrier functions is defined as follows,
ϕ i ± ( x ) = C ϵ i 1 / 2 B i l 1 ( ϰ ) + ϵ n 1 / 2 B n l 1 ( ϰ ) ± w i L 1 , ( ϰ ) , i = 1 , 2 , , n .
The bound on w i L 1 , ( ϰ ) is obtained from the equation of w i L 1 in (13). To derive the bounds w i L 1 , ( ϰ ) , by differentiating twice and thrice respectively and an argument analogous to the one used to bound w i L 1 , ( ϰ ) leads to the bounds, which are required. The bound on w i L 1 , ( ϰ ) is obtained by differentiating equation of w i L 1 in (13) once and using the bounds of w i L 1 , ( ϰ ) and w i L 1 , ( ϰ ) , this approach allows us to obtain
| w i L 1 , ( ϰ ) | C ϵ i 1 k = 1 i 1 ϵ k 1 / 2 B k l 1 ( ϰ ) + k = i n ϵ k 1 / 2 B k l 1 ( ϰ ) .
Similarly, it can be derived for w L 2 on [1,2]. The bounds on w R 1 is established and its derivatives for the case α μ 2 γ ϵ j . In this scenario, w R 1 is decomposed over the interval ( 0 , 1 ) .
w R 1 = q n R 1 + ϵ n p n R 1 + ϵ n 2 z n R 1 + ϵ n 3 y n R 1 ,
μ A ( ϰ ) q n R 1 , B ( ϰ ) q n R 1 = 0 ,
μ A ( ϰ ) p n R 1 , B ( ϰ ) p n R 1 = ϵ n 1 E q n ,
μ A ( ϰ ) z n R 1 , B ( ϰ ) z n R 1 = ϵ n 1 E p n ,
L 1 y n R 1 ( ϰ ) = ϵ n 1 E z n .
Since ϵ n 1 E is a matrix with bounded entries and hence for 0 k 3 , it implies that
q n R 1 , ( k ) C μ k , p n R 1 , ( k ) C μ ( k + 2 ) , z n R 1 , ( k ) C μ ( k + 4 ) .
Now using Theorem 3.1 and (50), for the choice of y n n R 1 ( 0 ) = 0 , then
| y n n R 1 , ( k ) | C ϵ n k μ k 6 .
Then from (51) and (52), it leads to
| w n R 1 , ( k ) | C μ k .
The decomposition for each component with respect to ϵ n is given , then apply ϵ n 1 to the first n 1 components, followed by ϵ n 2 for the first n 2 components, and so on. Thus, the following bound for all i in the range 1 i n and k in the range 0 k 3 is as follows | w i R 1 , ( k ) | C μ k . Utilizing the previous derivations, deriving | w R 1 ( ϰ ) | C B i r 1 ( ϰ ) C e γ μ ( 1 ϰ ) , similarly finding | w i R 1 , ( k ) ( ϰ ) | C μ k B i r 1 ( ϰ ) for k = 1 , 2 , 3 . Similarly derive for w R 2 , for the case α μ 2 γ ϵ i , the bounds on w R are obtained by the above analogous arguments to those applied in bounding w L . The proof of the theorem is complete. □

8. Sharper Bounds for w i L 1 and w i L 2

To achieve sharper bounds on the derivatives of the singular components w i L 1 and w i L 2 , These components are further decomposed for the intervals [0,1] and [1,2]. This refinement will help in demonstrating nearly first-order convergence of the proposed method. Now, the case α μ 2 γ ϵ j is focused. In addition to that, the following ordering holds
ϰ i , j ( k ) < ϰ i + 1 , j ( k ) , if i + 1 < j and ϰ i , j ( k ) < ϰ i , j + 1 ( k ) , if i < j .
For w L 1 and w L 2 , it is decomposed as follows, w i L 1 = ρ = 1 n w i , ρ L 1 on [0,1] and w i L 2 = ρ = 1 n w i , ρ L 2 on [1,2], where the components w i , ρ L 1 are defined, on the interval [ 0 , 1 ] , by
w i , n L 1 ( ϰ ) = k = 0 3 ( ( ϰ ϰ n 1 , n ( 2 ) ) k k ! ) w i L 1 , ( k ) ( ϰ n 1 , n ( 2 ) ) , on [ 0 , ϰ n 1 , n ( 2 ) ) , w i L 1 ( ϰ ) on ( ϰ n 1 , n ( 2 ) , 1 ] ,
for n 1 ρ i ,
w i , ρ L 1 ( ϰ ) = k = 0 3 ( ( ϰ ϰ ρ 1 , ρ ( 2 ) ) k k ! ) ϱ i , ρ ( k ) ( ϰ ρ 1 , ρ ( 2 ) ) , on [ 0 , ϰ ρ 1 , ρ ( 2 ) ) , ϱ i , ρ ( ϰ ) on ( ϰ ρ 1 , ρ ( 2 ) , 1 ] ,
and for i 1 ρ 2 ,
w i , ρ L 1 ( ϰ ) = k = 0 3 ( ( ϰ ϰ ρ 1 , ρ ( 1 ) ) k k ! ) ϱ i , ρ ( k ) ( ϰ ρ 1 , ρ ( 1 ) ) , on [ 0 , ϰ ρ 1 , ρ ( 1 ) ) , ϱ i , ρ ( ϰ ) on ( ϰ ρ 1 , ρ ( 1 ) , 1 ] .
On the interval [ 1 , 2 ] ,
w i , n L 2 ( ϰ ) = k = 0 3 ( ( ϰ ( ϰ n 1 , n ( 2 ) 1 ) ) k k ! ) w i L 2 , ( k ) ( ϰ n 1 , n ( 2 ) 1 ) , on [ 1 , ϰ n 1 , n ( 2 ) 1 ) , w i L 2 ( ϰ ) on ( ϰ n 1 , n ( 2 ) 1 , 2 ] ,
for n 1 ρ i ,
w i , ρ L 2 ( ϰ ) = k = 0 3 ( ( ϰ ( ϰ ρ 1 , ρ ( 2 ) 1 ) ) k k ! ) ϱ i , ρ ( k ) ( ϰ ρ 1 , ρ ( 2 ) 1 ) , on [ 1 , ϰ ρ 1 , ρ ( 2 ) 1 ) , ϱ i , ρ ( ϰ ) on ( ϰ ρ 1 , ρ ( 2 ) 1 , 2 ] ,
and for i 1 ρ 2 ,
w i , ρ L 2 ( ϰ ) = k = 0 3 ( ( ϰ ( ϰ ρ 1 , ρ ( 1 ) 1 ) k k ! ) ) ϱ i , ρ ( k ) ( ϰ ρ 1 , ρ ( 1 ) 1 ) , on [ ϰ ρ 1 , ρ ( 1 ) 1 , 1 ) , ϱ i , ρ ( ϰ ) on ( ϰ ρ 1 , ρ ( 1 ) 1 , 2 ] ,
with ϱ i , ρ = k = ρ + 1 n w i , k L p and w i , 1 L p = w i L p k = 2 n w i , k L p , p = 1 , 2 on [ 0 , 2 ] .
Lemma 8.1.
Given the decompositions of component w i , ρ L 1 for each ρ and i , 1 i n , 1 ρ n , satisfy the following estimates hold on [0,1]
| w i , ρ L 1 , ( ϰ ) | C μ 3 ϵ i 1 ϵ ρ 2 B ρ l 1 ( ϰ ) , if i ρ , | w i , ρ L 1 , ( ϰ ) | C μ 3 ϵ i 2 ϵ ρ 1 B ρ l 1 ( ϰ ) , if i > ρ ,
| w i , ρ L 1 , ( ϰ ) | C μ 2 ϵ i 1 ϵ ρ 1 B ρ l 1 ( ϰ ) , if i ρ < n , | w i , ρ L 1 , ( ϰ ) | C μ 2 ϵ i 2 B ρ l 1 ( ϰ ) , if i > ρ ,
| w i , ρ L 1 , ( ϰ ) | C μ ϵ i 1 B ρ l 1 ( ϰ ) , if i < ρ ,
and the following estimates hold on [1,2]
| w i , ρ L 2 , ( ϰ ) | C μ 3 ϵ i 1 ϵ ρ 2 B ρ l 2 ( ϰ ) , if i ρ , | w i , ρ L 2 , ( ϰ ) | C μ 3 ϵ i 2 ϵ ρ 1 B ρ l 2 ( ϰ ) , if i > ρ ,
| w i , ρ L 2 , ( ϰ ) | C μ 2 ϵ i 1 ϵ ρ 1 B ρ l 2 ( ϰ ) , if i ρ < n , | w i , ρ L 2 , ( ϰ ) | C μ 2 ϵ i 2 B ρ l 2 ( ϰ ) , if i > ρ ,
| w i , ρ L 2 , ( ϰ ) | C μ ϵ i 1 B ρ l 2 ( ϰ ) , if i < ρ .
Proof. For the interval [0,1], differentiating (54) thrice,
| w i , n L 1 , ( ϰ ) | = | w i L 1 , ( ϰ n 1 , n ( 2 ) ) | , on [ 0 , ϰ n 1 , n ( 2 ) ) , | w i L 1 , ( ϰ ) | on ( ϰ n 1 , n ( 2 ) , 1 ] ,
Then for ϰ [ 0 , ϰ n 1 , n ( 2 ) ) , using Theorem 7.1,
| w i , n L 1 , ( ϰ ) | C μ 3 ϵ i 1 k = 1 i 1 ϵ k 2 B k l 1 ( ϰ n 1 , n ( 2 ) ) + k = i n ϵ k 2 B k l 1 ( ϰ n 1 , n ( 2 ) ) .
Since ϰ k , n ( 2 ) ϰ n 1 , n ( 2 ) for k < n , ϵ k 2 B k l 1 ( ϰ n 1 , n ( 2 ) ) ϵ n 2 B n l 1 ( ϰ n 1 , n ( 2 ) ) , and hence
| w i , n L 1 , ( ϰ ) | C μ 3 ϵ i 1 ϵ n 2 B n l 1 ( ϰ n 1 , n ( 2 ) ) C μ 3 ϵ i 1 ϵ n 2 B n l 1 ( ϰ ) .
For ϰ [ ϰ n 1 , n ( 2 ) , 1 ] , | w i , n L 1 , ( ϰ ) | = | w i L 1 , ( x ) | C μ 3 ϵ i 1 k = 1 i 1 ϵ k 2 B k l 1 ( ϰ ) + k = i n ϵ k 3 B k l 1 ( ϰ ) . As ϰ ϰ n 1 , n ( 2 ) , ϵ k 2 B k l 1 ( ϰ ) ϵ n 2 B n l 1 ( ϰ ) , and hence for ϰ [ ϰ n 1 , n ( 2 ) , 1 ] , | w i , n L 1 , ( ϰ ) | C μ 3 ϵ i 1 ϵ n 2 B n l 1 ( ϰ ) . From (55) and (56), it is evident that for each ρ , n 1 ρ i and ϰ [ ϰ ρ , ρ + 1 ( 2 ) , 1 ]
w i , ρ L 1 ( ϰ ) = ϱ i , ρ ( ϰ ) = w i L 1 ( x ) k = ρ + 1 n w i , k L 1 ( ϰ ) = w i L 1 ( ϰ ) w i L 1 ( ϰ ) = 0 .
Differentiating (55) thrice on ϰ [ 0 , ϰ ρ 1 , ρ ( 2 ) ) , | w i , ρ L 1 , ( ϰ ) | = | ϱ i , ρ ( ϰ ρ 1 , ρ ( 2 ) ) | C μ 3 ϵ i 1 ϵ ρ 2 B ρ l 1 ( ϰ ) . For ϰ [ ϰ ρ 1 , ρ ( 2 ) , ϰ ρ , ρ + 1 ( 2 ) ) ,
| w i , ρ L 1 , ( ϰ ) | C μ 3 ϵ i 1 ϵ ρ 2 B ρ l 1 ( ϰ ) .
From (55) and (56), it is evident that for each ρ , i 1 ρ 2 , and ϰ [ ϰ ρ , ρ + 1 ( 1 ) , 1 ] , w i , ρ L 1 ( ϰ ) = 0 . Differentiating (56) thrice on ϰ [ 0 , ϰ ρ 1 , ρ ( 1 ) ) ,
| w i , ρ L 1 , ( ϰ ) | = | ϱ i , ρ ( ϰ ρ 1 , ρ ( 1 ) ) | C μ 3 ϵ i 1 k = 1 i 1 ϵ k 2 B k l 1 ( ϰ ρ 1 , ρ ( 1 ) ) + k = i n ϵ k 2 B k l 1 ( ϰ ρ 1 , ρ ( 1 ) )
C μ 3 ϵ i 2 ϵ ρ 1 B ρ l 1 ( ϰ ρ 1 , ρ ( 1 ) ) C μ 3 ϵ i 2 ϵ ρ 1 B ρ l 1 ( ϰ ) .
For ϰ [ ϰ ρ 1 , ρ ( 1 ) , ϰ ρ , ρ + 1 ( 1 ) ) ,
| w i , ρ L 1 , ( ϰ ) | C μ 3 ϵ i 2 ϵ ρ 1 B ρ l 1 ( ϰ ) .
From (56) and w i , 1 L 1 = w i L 1 k = 2 n w i , k L 1 , it is evident that w i , 1 L 1 ( ϰ ) = 0 for ϰ [ ϰ 1 , 2 ( 1 ) , 1 ] , and for ϰ [ 0 , ϰ 1 , 2 ( 1 ) ) , | w i , 1 L 1 , ( x ) | | w i L 1 , ( ϰ ) | C μ 3 ϵ i 1 k = 1 i 1 ϵ k 2 B k l 1 ( ϰ ) + k = i n ϵ k 2 B k l 1 ( ϰ )   C μ 3 ϵ i 2 ϵ 1 1 B 1 l 1 ( ϰ ) . Since w i , ρ L 1 , ( 1 ) = 0 for ρ < n , it can be concluded that for any ϰ [ 0 , 1 ] and i > ρ ,
| w i , ρ L 1 , ( ϰ ) | = | ϰ 1 w i , ρ L 1 , ( s ) d s | C μ 2 ϰ 1 ϵ i 2 ϵ ρ 1 B ρ l 1 ( s ) d s C μ 2 ϵ i 2 B ρ l 1 ( ϰ ) .
Hence, | w i , ρ L 1 , ( ϰ ) | C μ 2 ϵ i 2 B ρ l 1 ( ϰ ) , for i > ρ . Similar arguments lead to | w i , ρ L 1 , ( ϰ ) | C μ 2 ϵ i 1 ϵ ρ 1   B ρ l 1 ( ϰ ) , for i ρ , and | w i , ρ L 1 , ( ϰ ) | C μ ϵ i 1 B ρ l 1 ( ϰ ) , 1 i n , 1 ρ n . Similarly, it is not hard to find for the interval [1,2]. The proof of the Lemma is complete. □
Lemma 8.2.
Given the decompositions of component w i , ρ L 1 for each ρ and i , 1 i n , 1 ρ n , for α μ 2 γ ϵ i satisfy the following estimates hold on [ 0 , 2 ] , p = 1 for [ 0 , 1 ] and p = 2 for [ 1 , 2 ]
| w i , ρ L p , ( ϰ ) | C ϵ ρ 3 / 2 B ρ l p ( ϰ ) , if i ρ , | w i , ρ L p , ( ϰ ) | C ϵ ρ 3 / 2 B ρ l p ( ϰ ) , if i > ρ ,
| w i , ρ L 1 , ( ϰ ) | C ϵ ρ 1 B ρ l 1 ( ϰ ) , if i ρ < n , | w i , ρ L 1 , ( ϰ ) | C μ 2 ϵ ρ 1 B ρ l 1 ( ϰ ) , if i > ρ ,
Proof. The proof follows the same logic as Lemma 8.1. Analogously, the decompositions can be made for w R 1 and w R 2 in both case . The corresponding bounds for these components can be demonstrated in a similar manner.

9. Numerical Method

This section explains the numerical method proposed for (1).

9.1. Shishkin Mesh

For these cases α μ 2 γ ϵ i and α μ 2 γ ϵ j , appropriate Shishkin meshes are developed over interval [ 0 , 2 ] .
Case (i): α μ 2 γ ϵ i
A piecewise uniform Shishkin mesh is constructed over the interval [ 0 , 2 ] , the interval is divided into subintervals based on transition points as follows, [ 0 , 1 ] [ 1 , 2 ] [ n 1 , n ] [ n , 1 n ] [ 1 n , 1 n 1 ] ( 1 2 , 1 1 ] [ 1 1 , 1 ] [ 1 , 1 + 1 ] [ 1 + 1 , 1 + 2 ] [ 1 + n 1 , 1 + n ] [ 1 + n , 2 n ] [ 2 n , 2 n 1 ] ( 2 2 , 2 1 ] [ 2 1 , 2 ] . The transition points q for 1 q n are defined as
n = min 1 4 , 2 ϵ n γ α ln N , q = min q q + 1 q + 1 , 2 ϵ q γ α ln N
for q = n 1 , , 1 , ensuring finer mesh density near layer regions. The intervals are populated with points as follows, N 4 n points on all inner regions and for [ n , 1 n ] , a uniform mesh of N 2 is placed . If each q takes the left choice in its definition, the mesh becomes a classical uniform mesh, with q = q 4 n and a constant step size h j = N 1 . The step sizes in the intervals are defined as H 1 = 4 n N 1 , H q = 4 n N ( q q 1 ) for 2 q n , and H n + 1 = 2 N ( 1 2 * n ) . At each transition point q , the change in step size from h q to h q + is given by h q + h q = 4 n N q + 1 q ( d q d q 1 ) , where d q = q q + 1 q + 1 q , with d n = 0 when n = 1 4 . The mesh Ω N becomes a classical uniform mesh when d q = 0 for all q = 1 , , n , ensuring uniformly spaced transition points and a constant step size throughout the interval.Then, from (60), q C ϵ q ln N , 1 q n and also q = q m m , d q = = d m = 0 , 1 q m n .
Case (ii): α μ 2 γ ϵ j
A piecewise uniform Shishkin mesh is constructed over the interval [ 0 , 2 ] , the interval is divided into subintervals based on transition points as follows, [ 0 , 1 ] [ 1 , 2 ] [ n 1 , n ] [ n , 1 σ 1 ] [ 1 σ 1 , 1 ] [ 1 , 1 + 1 ] [ 1 + 1 , 1 + 2 ] [ 1 + n 1 , 1 + n ] [ 1 + n , 2 σ 1 ] [ 2 σ 1 , 2 ] . The transition points q for 1 q n are defined as
n = min 1 4 , 2 ϵ n μ α ln N , q = min q q + 1 q + 1 , 2 ϵ q μ α ln N , σ 1 = min 1 4 , μ γ ln N
for q = n 1 , , 1 , ensuring finer mesh density near layer regions. The intervals are populated with points as follows, N 4 n points on all inner regions, for [ 1 σ 1 , 1 ] a mesh of N 4 is placed and for [ n , 1 σ 1 ] a mesh of N 2 is placed . If each q takes the left choice in its definition, the mesh becomes a classical uniform mesh, with q = q 2 n and a constant step size h j = N 1 . The step sizes in the intervals are defined as H 1 = 4 n N 1 , H q = 4 n N ( q q 1 ) for 2 q n , H n + 1 = 2 N ( 1 σ 1 n ) and H s = 4 N σ 1 for [ 1 σ 1 , 1 ] . At each transition point q , the change in step size from h q to h q + is given by h q + h q = 4 n N q + 1 q ( d q d q 1 ) , where d q = q q + 1 q + 1 q , with d n = 0 when n = 1 4 . The mesh Ω N becomes a classical uniform mesh when d q = 0 for all q = 1 , , n , ensuring uniformly spaced transition points throughout the interval. Then, from (61), q C ϵ q ln N , 1 q n and also q = q m m , d q = = d m = 0 , 1 q m n .

10. The Discrete Problem

The discrete problem is defined as follows,
L N U ( ϰ j ) E δ 2 U ( ϰ j ) + μ A ( ϰ j ) D + U ( ϰ j ) B ( ϰ j ) U ( ϰ j ) + D ( ϰ j ) U ( ϰ j 1 ) = f ( ϰ j ) on Ω N ,
0 j N 1 , with boundary conditions specified as follows, U ( ϰ j 1 ) = φ ( ϰ j 1 ) , for 0 j N 2 , U ( ϰ N ) = u ( ϰ N ) , where, U = ( U 1 , U 2 2 , , U n ) T . Let
L 1 N U ( ϰ j ) = E δ 2 U ( ϰ j ) + μ A ( ϰ j ) D + U ( ϰ j ) B ( ϰ j ) U ( ϰ j ) = g ( ϰ j ) on Ω 1 N ,
L 2 N U ( ϰ j ) = E δ 2 U ( ϰ j ) + μ A ( ϰ j ) D + U ( ϰ j ) B ( ϰ j ) U ( ϰ j ) + D ( ϰ j ) U ( ϰ j 1 ) = f ( ϰ j ) on Ω 2 N .
The discrete derivatives are defined as follows
D + U ( ϰ j ) = U ( ϰ j + 1 ) U ( ϰ j ) h j + 1 , D U ( ϰ j ) = U ( ϰ j ) U ( ϰ j 1 ) h j , δ 2 U ( ϰ j ) = 1 h ¯ j ( D + U ( ϰ j ) D U ( ϰ j ) ) ,
with h j = ϰ j ϰ j 1 , h ¯ j = h j + h j + 1 2 , ϰ j Ω ¯ N .

11. Numerical Results

This section derives a discrete minimum principle, demonstrates a discrete stability analysis of the proposed numerical method, and proves its first-order convergence.
Lemma 11.1.(Discrete Minimum Principle) Assume that the mesh function Ψ ( ϰ j ) = ( Ψ 1 ( ϰ j ) , Ψ 2 ( ϰ j ) , , Ψ n ( ϰ j ) ) T satisfies Ψ ( ϰ 0 ) 0 and Ψ ( ϰ N ) 0 . Then, if L 1 N Ψ ( ϰ j ) 0 for 1 j N 2 1 and L 2 N Ψ ( ϰ j ) 0 for N 2 j N 1 and D + Ψ ( ϰ N 2 ) D Ψ ( ϰ N 2 ) 0 , it implies that Ψ ( ϰ j ) 0 for all 0 j N .
Proof. Let i * and j * be such that Ψ i * ( ϰ j * ) = min i , j Ψ i ( ϰ j ) and suppose Ψ i * ( ϰ j * ) < 0 . Then, j * { 0 , N } , Ψ i * ( ϰ j * ) Ψ i * ( ϰ j * + 1 ) , and Ψ i * ( ϰ j * ) Ψ i * ( ϰ j * 1 ) . Therefore, δ 2 Ψ i * ( ϰ j * ) 0 . Consider the two cases, for 1 j * N 2 1 , if ϰ j * Ω 1 N , then
( L 1 N Ψ ) i * ( ϰ j * ) = ϵ i * δ 2 Ψ i * ( ϰ j * ) + μ a i * ( ϰ j * ) D + Ψ i * ( ϰ j * ) j = 1 n b i * j ( ϰ j * ) Ψ j ( ϰ j * ) > 0
which is a contradiction, which gives ( L 1 N Ψ ) i * ( ϰ j * ) 0 . For N 2 j * N 1 ,if ϰ j * Ω 2 N , then
( L 2 N Ψ ) i * ( ϰ j * ) = ϵ i * δ 2 Ψ i * ( ϰ j * ) + μ a i * ( ϰ j * ) D + Ψ i * ( ϰ j * ) j = 1 n b i * j ( ϰ j * ) Ψ j ( ϰ j * ) d i * ( x j * ) Ψ i * ( ϰ j * 1 ) > 0
which is a contradiction, which gives ( L 2 N Ψ ) i * ( ϰ j * ) 0 . The only remaining possibility is that ϰ j * = ϰ N 2 . Thus, by hypothesis, it follows that D Ψ i * ( ϰ N 2 ) 0 and D + Ψ i * ( ϰ N 2 ) 0 . This implies D + Ψ i * ( ϰ N 2 ) D Ψ i * ( ϰ N 2 ) 0 , and since D Ψ i * ( ϰ N 2 ) 0 D + Ψ i * ( ϰ N 2 ) , it holds that D + Ψ i * ( ϰ N 2 ) D Ψ i * ( ϰ N 2 ) . Consequently, it follows that Ψ i * ( ϰ N 2 1 ) = Ψ i * ( ϰ N 2 ) = Ψ i * ( ϰ N 2 + 1 ) < 0 . Now, consider the operator acting on the solution at ϰ N / 2 1
( L 1 N Ψ ) i * ( ϰ N 2 1 ) = ϵ i * δ 2 Ψ i * ( ϰ N 2 1 ) + a i * ( ϰ N 2 1 ) Ψ i * ( ϰ N 2 1 ) j = 1 n b i * j ( ϰ N 2 1 ) Ψ j ( ϰ N 2 1 ) = ϵ i * δ 2 Ψ i * ( ϰ N 2 1 ) j = 1 , j i * n b i * j ( ϰ N 2 1 ) Ψ i * ( ϰ N / 2 1 ) > 0 ,
leading to a contradiction, implying that Ψ ( ϰ j ) 0 for all 0 j N . The proof of the lemma is complete. □
Lemma 11.2.(Discrete Stability Result) If Ψ ( ϰ j ) = ( Ψ 1 ( ϰ j ) , Ψ 2 ( ϰ j ) , , Ψ n ( ϰ j ) ) T is any mesh function, then
| Ψ i ( ϰ j ) | max | Ψ ( ϰ 0 ) | , | Ψ ( ϰ N ) | , max 1 j N 2 1 | L 1 N Ψ ( ϰ j ) | , max N 2 j N 1 | L 2 N Ψ ( ϰ j ) | .

11.1. Error Estimate

Analogous to the continuous case, the discrete solution U can be decomposed into V and W as defined below.
L 1 N V 1 ( ϰ j ) = g ( ϰ j ) , for 0 < j < N 2 1 , V 1 ( ϰ 0 ) = r ( 0 ) , V 1 ϰ N 2 = r ( 1 ) ,
L 2 N V 2 ( ϰ j ) = f ( ϰ j ) , for N 2 < j < N 1 , V 2 ϰ N 2 = s ( 1 ) , V 2 ( ϰ N ) = s ( 2 ) ,
V 2 ( ϰ j 1 ) = V 1 ϰ j N 2 , for N 2 j < N ,
L 1 N W L 1 ( ϰ j ) = 0 , for 0 < j < N 2 1 , W L 1 ( ϰ 0 ) = w L 1 ( 0 ) , W L 1 ϰ N 2 = w L 1 ( 1 ) ,
L 2 N W L 2 ( ϰ j ) = 0 , for N 2 < j < N 1 , W L 2 ϰ N 2 = w L 2 ( 1 ) , W L 2 ( ϰ N ) = w L 2 ( 2 ) ,
W L 1 ( ϰ j 1 ) = W L 1 ϰ j N 2 , for N 2 j < N .
It is clear that
V ( ϰ j ) = V 1 ( ϰ j ) , for 0 j N 2 1 V 2 ( ϰ j ) , for N 2 j N , W L ( ϰ j ) = W L 1 ( ϰ j ) , for 0 j N 2 1 W L 2 ( ϰ j ) , for N 2 j N .
Lemma 11.3.
If v is the solution of (9), (10) and (19) and V is the solution of (65) and (66) , then
| ( V v ) ( ϰ j ) | C N 1 , for 0 j N .

Proof.

| ( V v ) ( ϰ j ) | = | ( V 1 r ) ( ϰ j ) | , for 0 j N 2 1 | ( V 2 s ) ( ϰ j ) | , for N 2 j N ,
L 1 N ( V r ) ( x j ) = g ( ϰ j ) L 1 N r ( ϰ j ) = ( L 1 L 1 N ) r ( ϰ j )
= E d 2 d ϰ 2 δ 2 r ( ϰ j ) + μ d d ϰ D + A ( ϰ j ) r ( x j ) .
Determining the local truncation error
ϵ i d 2 d ϰ 2 δ 2 r i ( ϰ j ) + μ a i ( ϰ j ) d d ϰ D + r i ( ϰ j ) C ( ϰ j + 1 ϰ j 1 ) ( ϵ i r i + μ r i ) ,
for i = 1 , 2 , , n . It is established that ( ϰ j + 1 ϰ j 1 ) C N 1 . In this case α μ 2 γ ϵ i , from (39) and (42), | L 1 N ( V 1 r ) ( ϰ j ) | C N 1 , | L 2 N ( V 1 s ) ( ϰ j ) | C N 1 . Using Lemma 11.2, consider the mesh functions, Ψ ± ( ϰ j ) = C N 1 ( e T ) ± ( V 1 r ) ( ϰ j ) , 0 j N 2 . Provided that the value of C is sufficiently large, it follows that L 1 N Ψ ± ( ϰ j ) 0 , for 1 j N 2 1 , Ψ ± ( ϰ 0 ) 0 and Ψ ± ( ϰ N 2 ) 0 . Thus, | ( V 1 r ) ( ϰ j ) | C N 1 , for 0 j N 2 . Similarly, | ( V 1 s ) ( ϰ j ) | C N 1 , for N 2 j N . For the case α μ 2 γ ϵ j , | ( V 2 r ) ( ϰ j ) | C N 1 , for 0 j N 2 . Similarly, | ( V 2 s ) ( ϰ j ) | C N 1 , for N 2 j N . Thus,
| ( V v ) ( ϰ j ) | C N 1 , for 0 j N .
The proof of the lemma is complete. □
The bounds on the error in the singular components w L 1 and w L 2 are estimated for the case α μ 2 γ ϵ j . These estimates are derived utilizing the mesh functions B i ( l p , N ) ( ϰ j ) , where 1 i n , which are defined over Ω ¯ N ,
B i ( l p , N ) ( ϰ j ) = k = 1 j 1 + α μ h k 2 ϵ i 1 with B i ( l p , N ) ( ϰ 0 ) = 1 .
Lemma 11.4.
For the case α μ 2 γ ϵ j , the layer components W i L 1 and W i L 2 , 1 i n satisfy the following bounds on Ω ¯ N ,
| W i L 1 ( ϰ j ) | C B n ( l 1 , N ) ( ϰ j ) , | W i L 2 ( ϰ j ) | C B n ( l 2 , N ) ( ϰ j ) .
Proof. This result can be demonstrated by defining the appropriate mesh functions ψ i ± ( ϰ j ) = C B n ( l p , N ) ( ϰ j ) ± W i L p ( ϰ j ) , i = 1 , 2 , , n and p = 1 , 2 and noticing that ψ i ± ( ϰ 0 ) 0 and ψ i ± ( ϰ N ) 0 . Furthermore, ( L 1 N ψ ± ) i ( ϰ j ) 0 , and ( L 2 N ψ ± ) i ( ϰ j ) 0 , j = 1 , 2 , , N 1 . Consequently, the discrete minimum principle yields the expected result. The proof of the lemma is complete. □
Lemma 11.5.
Assume that d q = 0 , for q = 1 , 2 , , n . Let w L 1 and w L 2 satisfy (22), W L 1 and W L 2 satisfy (67) and (68). Then,
W L 1 w L 1 C N 1 ln N , W L 2 w L 2 C N 1 ln N .
Proof. The local truncation error is given by
| L 1 N ( W L 1 w L 1 ) ( ϰ j ) | C ( ϰ j + 1 ϰ j 1 ) ϵ i w i L 1 , D + μ w i L 1 , D
| L 2 N ( W L 2 w L 2 ) ( ϰ j ) | C ( ϰ j + 1 ϰ j 1 ) ϵ i w i L 2 , D + μ w i L 2 , D
where D = [ ϰ j 1 , ϰ j + 1 ] . Since d q = 0 ,the mesh Ω ¯ N is uniform, then the value of h = N 1 . In this instance, μ ϵ k 1 C ln N and μ 1 C ln N .
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = 1 i 1 ϵ k 1 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 B k l 1 ( ϰ j 1 ) ,
similarly, | ( L 2 N ( W L 2 w L 2 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = 1 i 1 ϵ k 1 B k l 2 ( ϰ j 1 ) + k = i n ϵ k 1 B k l 2 ( ϰ j 1 ) . Let the barrier function ϕ ( ϰ j ) = ( ϕ 1 ( ϰ j ) , ϕ 2 ( ϰ j ) , , ϕ n ( ϰ j ) ) T given by
ϕ i ( ϰ j ) = C N 1 ln N ν ( α ν ) ( k = 1 i 1 e 2 ν μ h ϵ k Y k ( ϰ j ) + k = i n e 2 ν μ h ϵ k Z k ( ϰ j ) ) ,
on Ω N , where ν is a constant and it satisfies 0 < ν < α , Y k ( ϰ j ) = λ k N j 1 λ k N 1 with λ k = 1 + ν μ h ϵ k , 1 k n , Z k ( ϰ j ) = ζ k N j 1 ζ k N 1 with ζ k = 1 + ν μ h ϵ k . The mesh functions described above is inspired by those constructed in [14]. Now, that 0 Y k ( ϰ j ) , Z k ( ϰ j ) 1 , ( ϵ k δ 2 + μ ν D + ) Y k ( ϰ j ) = 0 , ( ϵ k δ 2 + μ ν D + ) Z k ( ϰ j ) = 0 , D + Y k ( ϰ j ) ν μ ϵ k 1 exp ν μ ϰ j + 1 ϵ 1 1 and D + Z k ( ϰ j ) ν μ ϵ 2 1 exp ν μ ϰ j + 1 ϵ 2 1 . Then, define ψ ± ( ϰ j ) = ϕ ( ϰ j ) ± ( W L w L ) ( ϰ j ) . It is easy to observe that ψ ( ϰ j ) 0 , j = 0 , . . . . . , N and L 1 N ψ ( ϰ j ) 0 , L 2 N ψ ( ϰ j ) 0 , 1 j N 1 . Hence, by applying minimum principle, | ( W L 1 w L 1 ) ( ϰ j ) | C N 1 ln N . Simillarly, | ( W L 2 w L 2 ) ( ϰ j ) | C N 1 ln N . The proof of the Lemma is complete. □
Lemma 11.6.
Let w L 1 and w L 2 satisfy (22) , W L 1 and W L 2 satisfy (67) and (68). Then,
W L 1 w L 1 C N 1 ln N , W L 2 w L 2 C N 1 ln N .
Proof. The required result is established for each mesh point ϰ j ( 0 , 1 ) by partitioning the interval ( 0 , 1 ) as shown in figure Preprints 149368 i001 for 2 m n 1 . In each of these scenarios, first an estimate for the local truncation error is derived. This is followed by the formulation of a suitable barrier function, designed to capture the essential properties of the solution within a specified domain. By utilizing these barrier functions, the desired estimate is obtained.
Case (a): ϰ j ( 0 , 1 ) .
Clearly ϰ j + 1 ϰ j 1 C ϵ 1 μ 1 N 1 ln N . Then, utilizing the standard approach to local truncation through Taylor expansions, error estimates is obtained, which are valid for ϰ j ( 0 , 1 ) and 1 i n ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = 1 i 1 ϵ k 1 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 B k l 1 ( ϰ j 1 ) ,
| ( L 2 N ( W L 2 w L 2 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = 1 i 1 ϵ k 1 B k l 2 ( ϰ j 1 ) + k = i n ϵ k 1 B k l 2 ( ϰ j 1 ) .
Let the mesh functions be defined for ϰ j ( 0 , 1 ) , where 1 i n and p = 1 , 2 ,
ϕ i ( ϰ j ) = C N 1 ln N k = 1 i 1 e 2 ν μ H 1 ϵ k B k ( l p , N ) ( ϰ j ) + k = i n e 2 ν μ H 1 ϵ k B k ( l p , N ) ( ϰ j ) + k = 1 n B k ( l p , N ) ( k ) .
Utilizing the minimum principle and barrier function Ψ ± ( ϰ j ) = ϕ ( ϰ j ) ± ( W L 1 w L 1 ) ( ϰ j ) , it has been established that | ( W L 1 w L 1 ) ( ϰ j ) | C N 1 ln N , similarly for the interval (1,2), | ( W L 2 w L 2 ) ( ϰ j ) | C N 1 ln N .
Case (b): ϰ j [ 1 , 2 ) .
The two scenarios considered are Case (b1):  d 1 = 0 and Case (b2): d 1 > 0 . Case (b1):  d 1 = 0 , this case, the mesh is uniform within the interval ( 0 , 2 ) . Consequently, for any ϰ j + 1 ϰ j 1 C ϵ 1 μ 1 N 1 ln N , for ϰ j [ 1 , 2 ) . Then,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = 1 i 1 ϵ k 1 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 B k l 1 ( ϰ j 1 ) .
Now for ϰ j [ 1 , 2 ) and 1 i n , define,
ϕ i ( x j ) = C N 1 ln N k = 1 i 1 e 2 ν μ H 2 ϵ k B k ( l p , N ) ( ϰ j ) + k = i n e 2 ν μ H 2 ϵ k B k ( l p , N ) ( ϰ j ) + k = 2 n B k ( l p , N ) ( k ) .
Utilizing the minimum principle and barrier function Ψ ± ( ϰ j ) = ϕ ( ϰ j ) ± ( W L 1 w L 1 ) ( ϰ j ) , it has been derived that | ( W L 1 w L 1 ) ( ϰ j ) | C N 1 ln N , similarly for the interval (1,2), | ( W L 2 w L 2 ) ( ϰ j ) | C N 1 ln N .  Case (b2):  d 1 > 0 , for this case, ϰ j + 1 ϰ j 1 C ϵ 2 N 1 μ 1 ln N , and hence for ϰ j [ 1 , 2 ) , by utilizing the standard approach to local truncation errors in Taylor series expansions, h ¯ = ϰ j + 1 ϰ j 1 then,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C ϵ i μ | w i , 1 L 1 , ( 2 ) ( ϰ j 1 ) | + C ( ϰ j + 1 ϰ j 1 ) ϵ i k = 2 n | w i , k L 1 , ( 3 ) ( ϰ j 1 ) | + C μ | w i , 1 L 1 , ( 1 ) ( ϰ j 1 ) | + C ( ϰ j + 1 ϰ j 1 ) k = 2 n | w i , k L 1 , ( 2 ) ( ϰ j 1 ) | .
Now using Lemma 8.1 , it is not hard to derive that
| ( L 1 N ( W L 1 w L 1 ) ) 1 ( ϰ j ) | C N 1 ln N μ 2 k = 2 n ϵ k 1 B k l 1 ( ϰ j 1 ) + C μ 2 ϵ 1 1 B 1 l 1 ( ϰ j 1 ) ,
and for 2 i n ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 ln N μ 2 k = i n ϵ k 1 B k l 1 ( ϰ j 1 ) + C μ 2 ϵ i 1 B 1 l 1 ( ϰ j 1 ) .
Specify
ϕ 1 ( x j ) = C N 1 ln N k = 2 n exp ( 2 α μ H 2 / ϵ k ) B k ( l p , N ) ( x j ) + C B 1 ( l p , N ) ( x j ) + k = 2 n B k ( l p , N ) ( k )
and for 2 i n ,
ϕ i ( x j ) = C N 1 ln N k = i n exp ( 2 α μ H 2 / ϵ k ) B k ( l p , N ) ( x j ) + C B 1 ( l p , N ) ( x j ) + k = 2 n B k ( l p , N ) ( k ) .
Case (c): ϰ j [ m , m + 1 ) .
Here are the three scenarios Case (c1): d 1 = d 2 = = d m = 0 , Case (c2): d q > 0 and d q + 1 = = d m = 0 for some q, 1 q m 1 , Case (c3): d m > 0 . Case (c1):  d 1 = d 2 = = d m = 0 , since 1 = C m + 1 and the mesh remains uniform over the interval ( 0 , m + 1 ) , it can be concluded that for ϰ j ( m , m + 1 ] , ϰ j + 1 ϰ j 1 C ϵ 1 μ 1 N 1 ln N and hence
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = 1 i 1 ϵ k 1 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 B k l 1 ( ϰ j 1 ) ,
| ( L 2 N ( W L 2 w L 2 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = 1 i 1 ϵ k 1 B k l 2 ( ϰ j 1 ) + k = i n ϵ k 1 B k l 2 ( ϰ j 1 ) .
For 1 i n ,
ϕ i ( ϰ j ) = C N 1 ln N k = 1 i 1 e 2 ν μ H m + 1 ϵ k B k ( l p , N ) ( ϰ j ) + k = i n e 2 ν μ H m + 1 ϵ k B k ( l p , N ) ( ϰ j ) + k = m + 1 n B k ( l p , N ) ( k ) .
Utilizing the minimum principle and barrier function Ψ ± ( ϰ j ) = ϕ ( ϰ j ) ± ( W L 1 w L 1 ) ( ϰ j ) , it has been derived that | ( W L 1 w L 1 ) ( ϰ j ) | C N 1 ln N , similarly for the interval (1,2), | ( W L 2 w L 2 ) ( ϰ j ) | C N 1 ln N .  Case (c2):  d q > 0 and d q + 1 = = d m = 0 for some q, 1 q m 1 . In this case, since q + 1 = C m + 1 , the region ( q , m + 1 ) , exhibits a uniform mesh points in this region satisfies ϰ j + 1 ϰ j 1 C ϵ q + 1 N 1 μ 1 ln N , for any point ϰ j ( m , m + 1 ] . By utilizing the approach to local truncation is derived from Taylor expansions,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C ϵ i μ k = 1 q | w i , k L 1 , ( 2 ) ( ϰ j 1 ) | + C ( ϰ j + 1 ϰ j 1 ) ϵ i k = q + 1 n | w i , k L 1 , ( 3 ) ( ϰ j 1 ) | + C μ k = 1 q | w i , k L 1 , ( 1 ) ( ϰ j 1 ) | + C ( ϰ j + 1 ϰ j 1 ) k = q + 1 n | w i , k L 1 , ( 2 ) ( ϰ j 1 ) | .
Now, utilizing Lemma 8.1, it is evident that for i q ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = q + 1 n ϵ k 1 B k l 1 ( ϰ j 1 ) + C k = i q ϵ k 1 B k l 1 ( ϰ j 1 )
and for i > q ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C 1 N μ 2 ln N k = i n ϵ k 1 B k l 1 ( ϰ j 1 ) + C ϵ i 1 B q l 1 ( ϰ j 1 ) .
Now specify, for i q ,
ϕ i ( ϰ j ) = C 1 N ln N k = q + 1 n exp 2 α μ H m + 1 ϵ k B k ( l p , N ) ( ϰ j ) + C k = i q B k ( l p , N ) ( ϰ j ) + C k = m + 1 n B k ( l p , N ) ( k )
and for i > q ,
ϕ i ( ϰ j ) = C 1 N ln N k = i n exp 2 α μ H m + 1 ϵ k B k ( l p , N ) ( ϰ j ) + C B q ( l p , N ) ( ϰ j ) + C k = m + 1 n B k ( l p , N ) ( k ) .
Case (c3): d m > 0 . In the previous arguments of the case (c2), replacing q by m and applying the inequality x j + 1 x j 1 C ϵ m + 1 N μ 1 ln N , the estimates are valid for ϰ j ( m , m + 1 ] . For i m ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = m + 1 n ϵ k 1 B k l 1 ( ϰ j 1 ) + C k = i m ϵ k 1 B k l 1 ( ϰ j 1 )
and for i > m ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C 1 N μ 2 ln N k = i n ϵ k 1 B k l 1 ( ϰ j 1 ) + C ϵ i 1 B m l 1 ( ϰ j 1 ) .
For i m , define
ϕ i ( ϰ j ) = C 1 N ln N k = m + 1 n exp 2 α μ H m + 1 ϵ k B k ( l p , N ) ( ϰ j ) + C k = i m B k ( l p , N ) ( ϰ j ) + C k = m + 1 n B k ( l p , N ) ( k )
and for i > m ,
ϕ i ( ϰ j ) = C 1 N ln N k = i n exp 2 α μ H m + 1 ϵ k B k ( l p , N ) ( ϰ j ) + C B m ( l p , N ) ( ϰ j ) + C k = m + 1 n B k ( l p , N ) ( k ) .
Case (d):
There are three distinct cases to consider, Case (d1):  d 1 = = d n = 0 , Case (d2): d q > 0 and d q + 1 = = d n = 0 for some q, 1 q n 1 and Case (d3): d n > 0 . Case (d1): d 1 = = d n = 0 . In this case, the mesh is uniformly distributed over the interval [ 0 , 1 ] . The corresponding result for this situation is derived in Lemma 11.5. Case (d2):  d q > 0 and d q + 1 = = d n = 0 for some q, 1 q n 1 , for this scenario, based on the definition of n , it can be shown that ϰ j + 1 ϰ j 1 C ϵ q + 1 N 1 μ 1 ln N and by applying analogous arguments similar to Case (c2) lead to the estimates for ϰ j ( n , 1 ] . For i q ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 μ 2 ln N k = q + 1 n ϵ k 1 B k l 1 ( ϰ j 1 ) + C k = i q ϵ k 1 B k l 1 ( ϰ j 1 )
and for i > q ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C 1 N μ 2 ln N k = i n ϵ k 1 B k l 1 ( ϰ j 1 ) + C ϵ i 1 B q l 1 ( ϰ j 1 ) .
Now define, for i q ,
ϕ i ( ϰ j ) = C 1 N ln N k = q + 1 n exp 2 α μ H n + 1 ϵ k B k ( l p , N ) ( ϰ j ) + C k = i q B k ( l p , N ) ( ϰ j )
and for i > q ,
ϕ i ( ϰ j ) = C 1 N ln N k = i n exp 2 α μ H n + 1 ϵ k B k ( l p , N ) ( ϰ j ) + C B q ( l p , N ) ( ϰ j )
respectively. Case (d3): d n > 0 , let n be defined as n = ϵ n μ α ln N . Then, considering the interval ( n , 1 ] ,
| ( W i L 1 w i L 1 ) ( ϰ j ) | | W i L 1 ( ϰ j ) | + | w i L 1 ( ϰ j ) | C B n ( l 1 , N ) ( ϰ j ) + C B n l 1 ( ϰ j ) C B n ( l 1 , N ) ( n ) + C B n l 1 ( n ) C N 1 .
Hence, | ( W i L 1 w i L 1 ) ( ϰ j ) | C N 1 , and | ( W i L 2 w i L 2 ) ( ϰ j ) | C N 1 . Thus, for each of the cases, the barrier function is constructed and using minimum principle, it has been derived that | ( W L 1 w L 1 ) ( ϰ j ) | C N 1 ln N and | ( W L 2 w L 2 ) ( ϰ j ) | C N 1 ln N . Therefore,
| ( W L w L ) ( ϰ j ) | C N 1 ln N .
The proof of the lemma is complete.□
The bounds on the error in the singular components w L 1 and w L 2 are estimated for the case α μ 2 γ ϵ i . These estimates are derived utilizing the mesh functions B i ( l p , N ) ( ϰ j ) , where 1 i n , which are defined over Ω ¯ N ,
B i ( l p , N ) ( ϰ j ) = k = 1 j 1 + γ α ϵ i h k 1 ,
with B i ( l p , N ) ( ϰ 0 ) = 1 , for Ω 1 , p = 1 , for Ω 2 , p = 2 .
Lemma 11.7.
Let w L 1 and w L 2 satisfy (13), W L 1 and W L 2 satisfy (67) and (68). Then,
W L 1 w L 1 C N 1 ln N , W L 2 w L 2 C N 1 ln N .
Proof. Assume that d q = 0 , for q = 1 , 2 , , n , the local truncation error is given by
| L 1 N ( W L 1 w L 1 ) ( ϰ j ) | C ( ϰ j + 1 ϰ j 1 ) ϵ i w i L 1 , D + μ w i L 1 , D
| L 2 N ( W L 2 w L 2 ) ( ϰ j ) | C ( ϰ j + 1 ϰ j 1 ) ϵ i w i L 2 , D + μ w i L 2 , D
where D = [ ϰ j 1 , ϰ j + 1 ] . Since d q = 0 , the mesh Ω ¯ N is uniform, then the value of h = N 1 . In this instance, ϵ k 1 / 2 C ln N ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C ( ϰ j + 1 ϰ j 1 ) k = 1 i 1 ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) C N 1 ln N .
This is established for each mesh point ϰ j ( 0 , 1 ) by partitioning the interval ( 0 , 1 ) as follows Preprints 149368 i002 for 2 m n 1 . In each of these scenarios, first an estimate for the local truncation error. is derived. This is followed by the formulation of a suitable barrier function, designed to capture the essential properties of the solution within a specified domain. By utilizing these barrier functions, the desired estimate is obtained.
Case (a): ϰ j ( 0 , 1 ) .
Clearly ϰ j + 1 ϰ j 1 C ϵ 1 N 1 ln N then, utilizing the approach to local truncation in Taylor expansions error estimates is obtained, which are valid for ϰ j ( 0 , 1 ) and 1 i n ,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 ϵ 1 ln N k = 1 i 1 ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) C N 1 ln N ,
| ( L 2 N ( W L 2 w L 2 ) ) i ( ϰ j ) | C N 1 ϵ 1 ln N k = 1 i 1 ϵ k 1 / 2 B k l 2 ( ϰ j 1 ) + k = i n ϵ k 1 / 2 B k l 2 ( ϰ j 1 ) C N 1 ln N .
Case (b): ϰ j [ 1 , 2 ) .
The two scenarios considered are Case (b1):  d 1 = 0 and Case (b2): d 1 > 0 . Case (b1):  d 1 = 0 , this case, the mesh is uniform within the interval ( 0 , 2 ) . Consequently, for any ϰ j + 1 ϰ j 1 C ϵ 1 N 1 ln N , for ϰ j [ 1 , 2 ) . Then,
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 ϵ 1 ln N k = 1 i 1 ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) C N 1 ln N .
Case (b2): d 1 > 0 , for this case, ϰ j + 1 ϰ j 1 C ϵ 2 N 1 μ 1 ln N and hence for ϰ j [ 1 , 2 ) , by utilizing the standard approach to local truncation errors in Taylor series expansions, the term h ¯ = ϰ j + 1 ϰ j 1 then, using Lemma 8.2
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 ln N .
Case (c): x j [ m , m + 1 ) .
Here are the three scenarios Case (c1): d 1 = d 2 = = d m = 0 ,Case (c2): d q > 0 and d q + 1 = = d m = 0 for some q, 1 q m 1 and Case (c3): d m > 0 . Case (c1):  d 1 = d 2 = = d m = 0 , since 1 = C m + 1 and the mesh remains uniform over the interval ( 0 , m + 1 ) ,it can be concluded that for ϰ j ( m , m + 1 ] , ϰ j + 1 ϰ j 1 C ϵ 1 N 1 ln N and hence
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 ϵ 1 ln N k = 1 i 1 ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) C N 1 ln N .
Case (c2): d q > 0 and d q + 1 = = d m = 0 for some q, 1 q m 1 . Since q + 1 = C m + 1 , the mesh is uniform in ( q , m + 1 ) , which implies that ϰ j + 1 ϰ j 1 C ϵ q + 1 N 1 ln N , for ϰ j ( m , m + 1 ] . By utilizing the approach to local truncation is derived from Taylor expansions, outlined in Lemma 8.2
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 ln N .
Case (c3): d m > 0 . In the previous arguments of the case (c2), replacing q by m and applying the inequality ϰ j + 1 ϰ j 1 C ϵ m + 1 N 1 μ 1 ln N , the estimates are valid for ϰ j ( m , m + 1 ] .
| ( L 1 N ( W L 1 w L 1 ) ) i ( ϰ j ) | C N 1 ϵ m + 1 ln N k = 1 i 1 ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) + k = i n ϵ k 1 / 2 B k l 1 ( ϰ j 1 ) C N 1 ln N .
Case (d): There are three distinct cases to consider, Case (d1):  d 1 = = d n = 0 , Case (d2): d q > 0 and d q + 1 = = d n = 0 for some q, 1 q n 1 and Case (d3): d n > 0 . Case (d1): d 1 = = d n = 0 . In this case, the mesh is uniformly distributed over the interval [ 0 , 1 ] . The corresponding result for this situation is derived in Lemma 11.5. Case (d2):  d q > 0 and d q + 1 = = d n = 0 for some q, 1 q n 1 , for this scenario, based on the definition of n , it can be shown that ϰ j + 1 ϰ j 1 C ϵ q + 1 N 1 ln N and by employing arguments analogous to Case (c2), it leads to the estimates for x j ( n , 1 ] . Case (d3): d n > 0 . Let n be defined as n = ϵ n μ α ln N on the interval ( n , 1 ] . Hence, | ( W L 1 w L 1 ) ( ϰ j ) | C N 1 ln N and similarly, | ( W L 2 w L 2 ) ( ϰ j ) | C N 1 ln N . Therefore,
| ( W L w L ) ( ϰ j ) | C N 1 ln N .
The proof of the lemma is complete.□
To establish the bounds on the error | ( W R 1 w R 1 ) ( ϰ j ) | , the mesh function is defined over Ω ¯ N
B ( r p , N ) ( ϰ j ) = i = j + 1 N 1 + γ h i 2 μ 1 , B ( r p , N ) ( ϰ N ) = 1 , where B ( r p , N ) , p = 1 , 2 .
Lemma 11.8.
For the case α μ 2 γ ϵ j , the layer components W i R 1 and W i R 2 , 1 i n satisfy the following bounds on Ω ¯ N ,
| W i R 1 ( ϰ j ) | C B n ( r 1 , N ) ( ϰ j ) , | W i R 2 ( ϰ j ) | C B n ( r 2 , N ) ( ϰ j ) .
Proof. This result can be demonstrated by defining the mesh functions Ψ R 1 ( ϰ j ) = C B ( r 1 , N ) ( ϰ j ) ± W R 1 and Ψ R 2 ( ϰ j ) = C B ( r 2 , N ) ( ϰ j ) ± W R 2 . Also, since W R 1 ( 0 ) e γ μ , W R 1 ( 0 ) B ( r 1 , N ) ( ϰ 0 ) . Hence, Ψ R 1 ( 0 ) 0 . Also, for an appropriate choice of C, it follows that Ψ R 1 ( ϰ N ) 0 . Further, L 1 N Ψ R 1 ( ϰ j ) 0 and L 2 N Ψ R 2 ( ϰ j ) 0 . Hence, by the minimum principle , Ψ R 1 ( ϰ j ) 0 and Ψ R 2 ( ϰ j ) 0 , for 0 j N . Hence, | W R 1 ( ϰ j ) | C B ( r 1 , N ) ( ϰ j ) and | W R 2 ( ϰ j ) | C B ( r 2 , N ) ( ϰ j ) on Ω ¯ N . The proof of the lemma is complete.□
Lemma 11.9.
At each point ϰ j Ω ¯ N , | ( W R w R ) ( ϰ j ) | C N 1 ln N , for the case α μ 2 γ ϵ j .
Proof. The local truncation error is given by
| L 1 N ( W R 1 w R 1 ) ( ϰ j ) | C ( ϰ j + 1 ϰ j 1 ) ϵ i w i R 1 , D + μ w i R 1 , D
where D = [ ϰ j 1 , ϰ j + 1 ] , μ 1 C ln N . Consider the case d 1 = 0 then, ϰ j + 1 ϰ j 1 C μ N 1 ln N . Hence, | L 1 N ( W R 1 w R 1 ) ( ϰ j ) | C N 1 ln N ,   | L 2 N ( W R 2 w R 2 ) ( ϰ j ) | C N 1 ln N . Consider the case d 1 > 0 , ϰ j ( 0 , 1 σ 1 ] . Hence,
| ( W R 1 w R 1 ) i ( ϰ j ) | | W i R 1 ( ϰ j ) | + | w i R 1 ( ϰ j ) | C B ( r 1 , N ) ( ϰ j ) + C B i r 1 ( ϰ j ) C B ( r 1 , N ) ( σ 1 ) + C B i r 1 ( σ 1 ) C N 1 ,
for ϰ j ( 1 , 2 σ 1 ] , similarly like above
| ( W R 2 w R 2 ) i ( ϰ j ) | | W i R 2 ( ϰ j ) | + | w i R 2 ( ϰ j ) | C B ( r 2 , N ) ( σ 1 ) + C B 2 r 2 ( σ 1 ) C N 1 .
Examine the mesh region ( 1 σ 1 , 1 ] . It is known that h ¯ = ϰ j + 1 ϰ j 1 , then, ϰ j + 1 ϰ j 1 C μ N 1 ln N , | L 1 N ( W R 1 w R 1 ) ( ϰ j ) | C N 1 ln N . For ( 2 σ 1 , 2 ] , | L 2 N ( W R 2 w R 2 ) ( ϰ j ) | C N 1 ln N . The proof of the lemma is complete.□
Theorem 11.1.
Let u be the solution of (1) and U be the solution of (62)-(64). Then, for each mesh point ϰ j Ω ¯ N ,
U u Ω ¯ N C N 1 ln N ,
for both of the cases α μ 2 γ ϵ i and α μ 2 γ ϵ j .
Proof. The proof follows Lemmas 11.3, 11.5,11.7 and 11.9.

12. Numerical Illustration

12.1. Example

The numerically approximation of the solution to the following system on the interval ( 0 , 2 ) is obtained by applying the proposed method to both cases where α μ 2 γ ϵ i and α μ 2 γ ϵ j
E u ( ϰ ) + μ A ( ϰ ) u ( ϰ ) B ( ϰ ) u ( ϰ ) + D ( ϰ ) u ( ϰ 1 ) = f ( ϰ ) for all ϰ Ω = ( 0 , 2 ) ,
where, A ( ϰ ) = d i a g ( 0.5 , 0.5 , 0.5 ) , f ( ϰ ) = ( 1.0 , 2.0 , 1.5 ) T , B ( ϰ ) = 6.0 1.0 1.0 1.0 6.0 1.0 1.0 1.0 6.0 ,   D ( ϰ ) = d i a g ( 0.8 , 0.8 , 0.8 ) . To evaluate the order of convergence, maximum pointwise errors and error constants, a modified two-mesh algorithm was utilized. The results are summarized in Tables 1 and 2. As the parameter η decreases, the error stabilizes for each N, while the maximum pointwise error D N decreases and the observed order of convergence p N improves with increasing N, confirming the theoretical predictions. Figure 1 and 2 display the solution profiles for the n- system over the interval ( 0 , 1 ) . In Figure 1, corresponding to the condition μ 2 ϵ i γ α , boundary layers are observed for the components u i ( i = 1 , 2 , , n ) near ϰ = 0 , ϰ = 1 and ϰ = 2 consistent with theoretical expectations. On the other hand, Figure 2 illustrates the case where μ 2 ϵ j γ α . Here, layers are observed for u i near ϰ = 0 , while boundary layers emerge near ϰ = 2 and delay near ϰ = 1 . Log-log plots effectively illustrate the relationship between the number of mesh points N and the maximum pointwise errors, offering a clear depiction of convergence behavior. Figure 3 presents the maximum pointwise errors for different η values in cases 1 and 2, demonstrating how the error decreases as N increases. These plots validate theoretical predictions and emphasize the impact of η on the accuracy of the numerical method.
Table 1. Values of D ε N , D N , p N , p * and C p * N when ϵ 1 = η 32 , ϵ 2 = η 16 , ϵ 3 = η 8 , μ = η 8 for α μ 2 γ ϵ i
Table 1. Values of D ε N , D N , p N , p * and C p * N when ϵ 1 = η 32 , ϵ 2 = η 16 , ϵ 3 = η 8 , μ = η 8 for α μ 2 γ ϵ i
η Number of mesh points N
96 192 384 768
0.1E+00 0.4050E-02 0.2201E-02 0.1123E-02 0.5639E-03
0.1E-01 0.3832E-02 0.2522E-02 0.1509E-02 0.8641E-03
0.1E-02 0.3832E-02 0.2522E-02 0.1509E-02 0.8641E-03
0.1E-03 0.3832E-02 0.2522E-02 0.1509E-02 0.8641E-03
0.1E-04 0.3832E-02 0.2522E-02 0.1509E-02 0.8641E-03
D N 0.4050E-02 0.2522E-02 0.1509E-02 0.8641E-03
p N 0.6828E+00 0.7408E+00 0.8049E+00
C p N 0.2424E+00 0.2424E+00 0.2329E+00 0.2140E+00
The order of convergence p * = 0.6828 E + 00
Computed error constant, C p * N = 0.2424 E + 00
Table 2. Values of D ε N , D N , p N , p * and C p * N when ϵ 1 = η 64 , ϵ 2 = η 32 , ϵ 3 = η 16 , μ = η 4 for α μ 2 γ ϵ j
Table 2. Values of D ε N , D N , p N , p * and C p * N when ϵ 1 = η 64 , ϵ 2 = η 32 , ϵ 3 = η 16 , μ = η 4 for α μ 2 γ ϵ j
η Number of mesh points N
96 192 384 768
0.625E-01 0.1369E-01 0.6205E-02 0.1707E-02 0.4365E-03
0.156E-01 0.1420E-01 0.1353E-02 0.6080E-02 0.1669E-02
0.391E-02 0.1258E-01 0.1401E-01 0.1317E-01 0.5811E-02
0.977E-03 0.4006E-01 0.2536E-01 0.1383E-01 0.1242E-01
0.244E-03 0.8087E-01 0.6491E-01 0.4528E-01 0.2625E-01
D N 0.8087E-01 0.6491E-01 0.4528E-01 0.2625E-01
p N 0.3170E+00 0.5195E+00 0.7862E+00
C p N 0.1742E+01 0.1742E+01 0.1514E+01 0.1093E+00
The order of convergence p * = 0.3170 E + 00
Computed error constant, C p * N = 0.1742 E + 01

13. Conclusions

This paper presented a robust fitted mesh finite difference method for solving a system of two-parameter n singularly perturbed delay differential equations of convection-reaction-diffusion type. The method leverages a piecewise uniform Shishkin mesh to address the intricate challenges posed by small perturbation parameters and delay terms across multiple equations. Our theoretical analysis demonstrates that the proposed scheme attains nearly first-order convergence in the maximum norm, uniformly with respect to the perturbation parameters. Numerical experiments confirm the method’s robustness and accuracy, demonstrating its capability to resolve boundary layers with precision across a system of equations. This work marks a significant advancement in numerical techniques for SPDDEs, emphasizing the critical importance of developing parameter-uniform methods to address the unique challenges posed by systems of equations with multiple layers and delays. Future investigations could focus on extending these methods to enhance computational efficiency, improve convergence rates and handle more intricate systems encountered in real-world applications.

References

  1. Bhatti, M. M. , Alamri, S. Z., Ellahi, R., & others. (2021). Intra - uterine particle - fluid motion through a compliant asymmetric tapered channel with heat transfer. Journal of Thermal Analysis and Calorimetry, 144, 2259–2267. [CrossRef]
  2. Glizer, V. (2003). Asymptotic analysis and solution of a finite-horizon H control problem for singularly-perturbed linear systems with small state delay. Journal of Optimization Theory and Applications, 117, 295–325. [CrossRef]
  3. Miller, J. J. H. , O’Riordan, E., & Shishkin, G. I. (2012). Fitted numerical methods for singular perturbation problems: error estimates in the maximum norm for linear problems in one and two dimensions. World Scientific.
  4. Doolan, E. P. , Miller, J. J. H., & Schilders, W. H. A. (1980). Uniform numerical methods for problems with initial and boundary layers, Boole Press.
  5. Cen, Z. (2005). Parameter-uniform finite difference scheme for a system of coupled singularly perturbed convection - diffusion equations. International Journal of Computer Mathematics, 82, 177–192. [CrossRef]
  6. Gracia, J. L. , O’Riordan, E., & Pickett, M. L. (2006). A parameter robust second order numerical method for a singularly perturbed two-parameter problem. Applied Numerical Mathematics, 56, 962–980. [CrossRef]
  7. O’Malley, R. E. (1967). Two-parameter singular perturbation problems for second-order equations. Journal of Mathematics and Mechanics, 16, 1143–1164.
  8. O’Riordan, E. , Pickett, M. L., & Shishkin, G. I. (2003). Singularly perturbed problems modeling reaction-convection-diffusion processes. Computational Methods in Applied Mathematics, 3, 424–442.
  9. Selvaraj, D. , & Mathiyazhagan, J. P. (2021). A parameter uniform convergence for a system of two singularly perturbed initial value problems with different perturbation parameters and Robin initial conditions. Malaya Journal of Matematik, 9, 498–505.
  10. Kalaiselvan, S. S. , Miller, J. J. H., & Sigamani, V. (2019). A parameter uniform numerical method for a singularly perturbed two-parameter delay differential equation. Applied Numerical Mathematics, 145, 90–110. [CrossRef]
  11. Nagarajan, S. (2022). A parameter robust fitted mesh finite difference method for a system of two reaction-convection-diffusion equations. Applied Numerical Mathematics, 179, 87–104. [CrossRef]
  12. Arthur, J. , Chatzarakis, G. E., Panetsos, S. L., & Mathiyazhagan, J. P. (2025). A robust-fitted-mesh-based finite difference approach for solving a system of singularly perturbed convection–diffusion delay differential equations with two parameters. Symmetry, 17, 68. [CrossRef]
  13. Mathiyazhagan, P. , Sigamani, V., & Miller, J. J. H. (2010). Second order parameter-uniform convergence for a finite difference method for a singularly perturbed linear reaction-diffusion system. Mathematical Communications, 15, 587–612.
  14. Farrell, P. Robust computational techniques for boundary layers (1st ed.). Chapman and Hall/CRC. [CrossRef]
Figure 1. Graphical representation of Numerical solutions for the case: α μ 2 γ ϵ i
Figure 1. Graphical representation of Numerical solutions for the case: α μ 2 γ ϵ i
Preprints 149368 g001
Figure 2. Graphical representation of Numerical solutions for the case: α μ 2 γ ϵ j
Figure 2. Graphical representation of Numerical solutions for the case: α μ 2 γ ϵ j
Preprints 149368 g002
Figure 3. Graphical representation of maximum pointwise errors for different η values for the cases 1 and 2
Figure 3. Graphical representation of maximum pointwise errors for different η values for the cases 1 and 2
Preprints 149368 g003
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

Disclaimer

Terms of Use

Privacy Policy

Privacy Settings

© 2025 MDPI (Basel, Switzerland) unless otherwise stated