António Camunga Tchikoko^1*, David Benje Mucenje², Euclides Faustino da Costa Fernando^3
1Department of Basic Sciences, Institute of Information and Communication Technologies, University of Luanda, Luanda, Angola.
²Department of Exact and Natural Sciences, Higher Institute of Educational Sciences of Benguela, Benguela, Angola.
³Department of Exact and Natural Sciences, Higher Institute of Educational Sciences of Huambo, Huambo, Angola.
DOI: 10.4236/jamp.2026.1410188   PDF    HTML   XML   12 Downloads   96 Views   Citations

Abstract

We study a singularly perturbed Cauchy problem for a second order linear ordinary differential equation with constant coefficients, in which a small parameter $\epsilon >0$ multiplies the highest derivative. Because the reduced first order equation cannot generally satisfy both initial conditions, a layer of width $O\left(\epsilon \right)$ develops near the initial point. An exact Green’s function representation, combined with the separation of the characteristic roots into slow and fast modes, yields an $O\left(\epsilon \right)$ error bound for the solution in the uniform norm on the full interval. For the derivative, we isolate an exponentially localized term proportional to the initial-data mismatch and bound the remaining contribution by $O\left(\epsilon \right)$ . When the initial slope is compatible with the reduced equation, the order one fast component is absent and convergence holds in $C^1$ with order $O\left(\epsilon \right)$ . A closed form example proves that the uniform solution rate and the order one derivative discrepancy at the initial point are optimal. We also state the additional kernel and spectral assumptions needed for analogous results for variable coefficients, higher order equations, and systems.

Keywords

Singular Perturbations, Initial Layer, Error Estimates, Uniform Convergence, Green’s Function

Share and Cite:

1. Introduction

Singularly perturbed differential equations arise in mathematical models in which processes evolving on widely separated time or spatial scales interact. Their defining feature is the presence of a small positive parameter multiplying the highest-order derivative. As this parameter approaches zero, the order of the differential equation is reduced, and the limiting problem can no longer accommodate all the conditions imposed on the original equation. This loss of order generally produces thin transition regions, known as boundary or initial layers, where the solution or its derivatives vary rapidly.

These problems occur in fluid mechanics, reaction kinetics, control theory, electrical circuits, transport phenomena, and other areas involving multiple timescales. From the analytical point of view, the principal difficulty is that convergence to the reduced solution is usually nonuniform near the layer. Although the solution itself may converge uniformly, its derivatives may retain an order-one discrepancy within a region whose thickness tends to zero with the perturbation parameter.

For linear ordinary differential equations, the qualitative convergence theory is well established [1]-[3]. Nevertheless, qualitative convergence alone does not determine the rate of approximation, the dependence of the estimates on the data, or the precise contribution of the rapidly decaying characteristic mode. These questions are essential for rigorous asymptotic analysis and for the construction of numerical methods whose accuracy remains uniform with respect to the perturbation parameter.

The purpose of this paper is to derive explicit and parameter-uniform error estimates for a second-order singularly perturbed Cauchy problem with constant coefficients. Our approach combines an exact Green’s function representation with a spectral decomposition of the characteristic roots into slow and fast modes. This makes it possible to identify the initial-layer contribution explicitly and to determine the compatibility condition under which that contribution disappears.

The main contributions of the paper are as follows:

(i) an explicit and rigorous $O\left(\epsilon \right)$ uniform estimate for the difference between the perturbed and reduced solutions;

(ii) a pointwise decomposition of the derivative error into an exponentially localized initial-layer term and an $O\left(\epsilon \right)$ remainder;

(iii) an exact compatibility condition ensuring uniform $C^1$ convergence;

(iv) numerical experiments comparing exact errors with the explicit theoretical bounds;

(v) a closed-form example demonstrating the sharpness of the estimates.

The remainder of the paper is organized as follows. Section 2 introduces the singularly perturbed problem and its reduced counterpart. Section 3 establishes the functional setting and regularity of the reduced solution. Section 4 analyzes the characteristic roots and the separation of timescales, while Section 5 derives the error equation and its Green’s-function representation. Sections 6 and 7 establish the uniform solution and derivative estimates. Section 8 states precise prototypes for extensions, Section 9 presents numerical experiments, and Section 10 gives an explicit example confirming sharpness. Conclusions are collected in Section 11.

2. Problem Formulation

2.1. Problem Statement

Fix $b>0$ and consider the family of linear, constant-coefficient, second-order Cauchy problems

$\epsilon a_2{y^{″}}_{\epsilon }\left(x\right)+a_1{y^{\prime }}_{\epsilon }\left(x\right)+a_0{y}_{\epsilon }\left(x\right)=f\left(x\right),0<x\le b,$ (2.1)

$y_{\epsilon }\left(0\right)=A_0,{y^{\prime }}_{\epsilon }\left(0\right)=A_1,$ (2.2)

under the standing hypotheses

$a_0>0,a_1>0,a_2>0,f\in C^2\left(\left[0,b\right]\right),$ (2.3)

Define once and for all the admissible parameter threshold

${\epsilon }_{\text{adm}}:=\text{min}\left\{1,\frac{a_1^2}{8a_0{a}_2}\right\}.$ (2.4)

We take $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ . This single threshold guarantees real, distinct, negative characteristic roots, the quantitative root separation used in the Green’s-function estimates, and constants independent of $\epsilon$ . As $\epsilon \to 0^{+}$ , one formally discards the term $\epsilon a_2{y^{″}}_{\epsilon }$ and is left with the reduced first-order problem

$a_1{y}_{0\text{'}}+a_0{y}_0=f,y_0\left(0\right)=A_0.$ (2.5)

Problem (2.5) is of first order, so it can only absorb one of the two initial data in (0.2); the datum $y_0\left(0\right)=A_0$ is imposed, but in general

${y^{\prime }}_0\left(0\right)=\frac{f\left(0\right)-a_0{A}_0}{a_1}\ne A_1.$

This mismatch generates an initial layer of width $O\left(\epsilon \right)$ near $x=0$ , where the derivative of the perturbed solution differs substantially from that of the reduced solution. We address the following questions:

(Q1) Uniform solution error. At what rate does $y_{\epsilon }$ converge to $y_0$ on the full interval $\left[0,b\right]$ ?

(Q2) Derivative behavior. Can ${y^{\prime }}_{\epsilon }-{y^{\prime }}_0$ be decomposed into an exponentially localized layer and a uniformly small remainder?

(Q3) Compatibility. Which condition on $\left(A_0,A_1,f\right)$ removes the order-one fast component and restores convergence in $C^1\left(\left[0,b\right]\right)$ ?

(Q4) Optimality. Are the resulting orders best possible?

2.2. State of the Art

The qualitative theory of (2.1)-(2.2) belongs to the classical literature on singular perturbations [4]. O’Malley’s analysis of singular initial-value problems gives uniform asymptotic expansions on finite intervals [5]; more recently, parameter-uniform discretizations and initial-layer correctors have been developed for stiff initial-value problems [6] [7]. Green’s-function estimates with explicit dependence on the perturbation parameter also play a central role in the analysis of layer problems [8]. Closely related is the analysis of incompatible initial and boundary data for singularly perturbed parabolic problems, where an interior or initial layer of order one similarly arises from a data mismatch and is treated by an explicit corrector function [9]. These peer-reviewed results provide the closest analytical and computational context for the present estimates.

In this paper, a problem is called regular if the solutions of the full Cauchy problems converge in $C\left(\left[0,b\right]\right)$ to the solution of the reduced problem as $\epsilon \to 0^{+}$ ; this is precisely the qualitative uniform-convergence property referred to in Question 1 above. For the scalar family considered here, the uniform-convergence criterion quoted below is used under the precise assumptions that the coefficients and initial data are independent of $\epsilon$ , $a_2\ne 0$ , $a_1\ne 0$ , $f\in C\left(\left[0,b\right]\right)$ , the reduced problem has a unique solution, and the root lost in the reduction is uniformly stable, namely $\mathrm{Re}{r}_{\epsilon }^{-}\le -c/\epsilon$ for some $c>0$ . The standing assumptions (2.3), together with $0<\epsilon \le {\epsilon }_{\text{adm}}$ from (2.4), imply all these conditions.

• Batard [10] [11] studies the complementary characteristic equation and, using a continuous-convergence theorem attributed to Elsgoltz [12], obtains a qualitative criterion for uniform convergence. In the present positive-coefficient setting, this criterion is consistent with the stability of the characteristic roots:

$y_{\epsilon }\to y_0\text{uniformly}\text{on}\left[0,b\right]\text{as}\epsilon \to 0^{+}\iff \text{the}\text{problem}\text{is}\text{regular}$

under the terminology and hypotheses used in those references. The result establishes convergence but does not provide the parameter-uniform bounds derived below.

• Comet [13] [14] develops approximate error formulae for simplified second- and third-order linear equations. Those formulae are useful for asymptotic and numerical interpretation. The present work takes a different route: it derives inequalities directly from the Green’s kernel. This distinction leads to three specific additions:

(L1) a non-asymptotic upper bound for the solution error;

(L2) a separate estimate for the derivative that identifies its fast component;

(L3) a $C^1$ estimate under the compatibility condition $\text{Δ}=0$ .

These observations define the scope of the contribution without suggesting that the broader theory of initial layers is new. The novelty claimed here is the particular set of parameter-uniform estimates, explicit constants, and the accompanying optimality calculation for problem (2.1)-(2.2).

2.3. Main Contributions

Using the Green’s function of ${ℒ}_{\epsilon }=\epsilon a_2{D}^2+a_{1}D+a_0$ and the separation of the characteristic roots into fast and slow modes, we obtain:

(C1) (Theorem 6.1) the non-asymptotic bound

${\Vert y_{\epsilon }-y_0\Vert }_{\infty }\le C\epsilon \left(\left|\text{Δ}\right|+{\Vert {y^{″}}_0\Vert }_{\infty }\right),\forall \epsilon \in \left(0,{\epsilon }_{\text{adm}}\right],$

valid on the full interval, with $C$ independent of $\epsilon$ ;

(C2) (Theorem 7.1) the derivative estimate

$\left|{y^{\prime }}_{\epsilon }\left(x\right)-{y^{\prime }}_0\left(x\right)\right|\le C\left|\Delta \right|{\text{e}}^{-\kappa x/\epsilon }+C\epsilon \left(\left|\Delta \right|+{\Vert {y^{″}}_0\Vert }_{\infty }+{\Vert {y^{‴}}_0\Vert }_{\infty }\right),$

which separates the order-one fast contribution from an $O\left(\epsilon \right)$ remainder and yields uniform derivative convergence away from $x=0$ ;

(C3) (Corollary 7.4) the compatibility implication

$\Delta :=A_1-{y^{\prime }}_0\left(0\right)=0\Rightarrow {\Vert y_{\epsilon }-y_0\Vert }_{C^1\left(\left[0,b\right]\right)}=O\left(\epsilon \right);$

(C4) (Example 10.1) a closed-form calculation showing that the $O\left(\epsilon \right)$ uniform solution error and the order-one derivative discrepancy at $x=0$ are optimal in general.

(C5) A precise discussion of the additional kernel and spectral hypotheses required to extend the framework to variable coefficients, higher-order equations, and systems (Section 8).

All estimates follow from the exact representation. Every main result below states explicitly which quantities its constants may depend on; none of those constants depends on $\epsilon$ .

3. Functional Setting and Reduced Problem

Let $C^k\left(\left[0,b\right]\right)$ denote the Banach space of $k$ -times continuously differentiable functions with the sup-norm ${\Vert \cdot \Vert }_{\infty }$ . Constants are defined locally in each result and are independent of $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ .

Proposition 3.1 (Reduced Solution Regularity) Under $f\in C^2\left(\left[0,b\right]\right)$ , the reduced problem (2.5) admits a unique solution $y_0\in C^3\left(\left[0,b\right]\right)$ given explicitly by

$y_0\left(x\right)=A_0{\text{e}}^{-\frac{a_0}{a_1}x}+\frac{1}{a_1}{\displaystyle {\int }_0^x{\text{e}}^{-\frac{a_0}{a_1}\left(x-s\right)}f\left(s\right)\text{d}s}.$ (3.1)

Moreover, $y_0^{\left(k\right)}=\left(f^{\left(k-1\right)}-a_0{y}_0^{\left(k-1\right)}\right)/a_1$ for $k=1,2,3$ .

Proof. Equation (2.5) is linear first-order with constant coefficient $a_0/a_1$ . Multiplying by the integrating factor ${\text{e}}^{\left(a_0/a_1\right)x}$ ,

$\frac{\text{d}}{\text{d}x}\left({\text{e}}^{\frac{a_0}{a_1}x}{y}_0\left(x\right)\right)={\text{e}}^{\frac{a_0}{a_1}x}\frac{f\left(x\right)}{a_1},$

and integrating from 0 to $x$ , using $y_0\left(0\right)=A_0$ ,

${\text{e}}^{\frac{a_0}{a_1}x}{y}_0\left(x\right)-A_0=\frac{1}{a_1}{\displaystyle {\int }_0^x{\text{e}}^{\frac{a_0}{a_1}s}f\left(s\right)\text{d}s}\Rightarrow y_0\left(x\right)=A_0{\text{e}}^{-\frac{a_0}{a_1}x}+\frac{1}{a_1}{\displaystyle {\int }_0^x{\text{e}}^{-\frac{a_0}{a_1}\left(x-s\right)}f\left(s\right)\text{d}s},$

which is (3.1). Since $f\in C^2$ , differentiating (2.5) once gives $a_1{y^{″}}_0+a_0{y^{\prime }}_0=f^{\prime }$ , so ${y^{″}}_0=\left(f^{\prime }-a_0{y^{\prime }}_0\right)/a_1\in C^1\left(\left[0,b\right]\right)$ , i.e. $y_0\in C^3\left(\left[0,b\right]\right)$ ; differentiating once more gives ${y^{‴}}_0=\left(f^{″}-a_0{y^{″}}_0\right)/a_1\in C\left(\left[0,b\right]\right)$ , confirming the stated recurrence for $k=1,2,3$ . □

Theorem 3.2 (Well-Posedness of the Full Problem) For every fixed $\epsilon >0$ , the problem (2.1)-(2.2) possesses a unique solution $y_{\epsilon }\in C^2\left(\left[0,b\right]\right)$ .

Proof. Writing $Y={\left(y_{\epsilon },{y^{\prime }}_{\epsilon }\right)}^{\text{T}}$ , equation (2.1) is equivalent to the linear first-order system

$Y^{\prime }=\left(\begin{array}{cc}0& 1\\ -\frac{a_0}{\epsilon a_2}& -\frac{a_1}{\epsilon a_2}\end{array}\right)Y+\left(\begin{array}{c}0\\ f\left(x\right)/\left(\epsilon a_2\right)\end{array}\right),Y\left(0\right)={\left(A_0,A_1\right)}^{\text{T}}.$

The matrix is constant and $f/\left(\epsilon a_2\right)\in C\left(\left[0,b\right]\right)$ , so the right-hand side is globally Lipschitz in $Y$ and continuous in $x$ . Standard existence and uniqueness theory for linear systems therefore gives a unique solution on $\left[0,b\right]$ , with $Y\in C^1\left(\left[0,b\right];{ℝ}^2\right)$ and hence $y_{\epsilon }\in C^2\left(\left[0,b\right]\right)$ . □

4. Spectral Analysis and Timescale Separation

The characteristic polynomial of ${\mathcal{G}}_{\epsilon }=\epsilon a_2{D}^2+a_{1}D+a_0$ is

$p_{\epsilon }\left(r\right)=\epsilon a_2{r}^2+a_{1}r+a_0,$

with roots

$r_{\epsilon }^{\pm }=\frac{-a_1\pm \sqrt{a_1^2-4\epsilon a_2{a}_0}}{2\epsilon a_2}.$ (4.1)

Since $\epsilon \le {\epsilon }_{\text{adm}}\le a_1^2/\left(8a_0{a}_2\right)$ , the discriminant is positive. Thus both roots are real, negative (by Vieta: $r_{\epsilon }^{+}+r_{\epsilon }^{-}=-a_1/\left(\epsilon a_2\right)<0$ and $r_{\epsilon }^{+}{r}_{\epsilon }^{-}=a_0/\left(\epsilon a_2\right)>0$ ), and distinct.

Proposition 4.1 (Rigorous Root Bounds and Asymptotic Expansion) Let ${\kappa }_1:=a_0/a_1$ and ${\kappa }_2:=a_1/\left(2a_2\right)$ . Then, for every $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ ,

$-\frac{2a_0}{a_1}\le r_{\epsilon }^{+}\le -{\kappa }_1,r_{\epsilon }^{-}\le -\frac{{\kappa }_2}{\epsilon }.$ (4.2)

Moreover, as $\epsilon \to 0^{+}$ , the roots admit the uniformly valid asymptotic expansions

$r_{\epsilon }^{+}=-\frac{a_0}{a_1}-\epsilon \frac{a_2{a}_0^2}{a_1^3}+O\left({\epsilon }^2\right),$ (4.3)

$r_{\epsilon }^{-}=-\frac{a_1}{\epsilon a_2}+\frac{a_0}{a_1}+\epsilon \frac{a_2{a}_0^2}{a_1^3}+O\left({\epsilon }^2\right).$ (4.4)

Proof. Step 1 (an exact closed form for $r_{\epsilon }^{+}$ , and the bound (4.2)). Rationalize the numerator of $r_{\epsilon }^{+}$ in (4.1) by multiplying and dividing by $\left(a_1+\sqrt{a_1^2-4\epsilon a_2{a}_0}\right)$ :

$\begin{array}{c}{r}_{\epsilon }^{+}=\frac{-a_1+\sqrt{a_1^2-4\epsilon a_2{a}_0}}{2\epsilon a_2}=\frac{{\left(\sqrt{a_1^2-4\epsilon a_2{a}_0}\right)}^2-a_1^2}{2\epsilon a_2\left(a_1+\sqrt{a_1^2-4\epsilon a_2{a}_0}\right)}\\ =\frac{-4\epsilon a_2{a}_0}{2\epsilon a_2\left(a_1+\sqrt{a_1^2-4\epsilon a_2{a}_0}\right)},\end{array}$

which simplifies, after cancelling $2\epsilon a_2\ne 0$ , to the exact, $\epsilon$ -uniform identity

$r_{\epsilon }^{+}=-\frac{2a_0}{a_1+\sqrt{a_1^2-4\epsilon a_2{a}_0}} .$ (4.5)

For $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ the quantity ${\sigma }_{\epsilon }:=\sqrt{a_1^2-4\epsilon a_2{a}_0}$ satisfies $a_1/\sqrt{2}\le {\sigma }_{\epsilon }<a_1$ , hence

$a_1\le a_1+{\sigma }_{\epsilon }\le 2a_1\Rightarrow -\frac{2a_0}{a_1}\le -\frac{2a_0}{a_1+{\sigma }_{\epsilon }}\le -\frac{a_0}{a_1},$

which is exactly the bound on $r_{\epsilon }^{+}$ in (4.2), with no asymptotic approximation involved.

Step 2 (bound on $r_{\epsilon }^{-}$ ). By Vieta’s formula $r_{\epsilon }^{+}+r_{\epsilon }^{-}=-a_1/\left(\epsilon a_2\right)$ , so

$r_{\epsilon }^{-}=-\frac{a_1}{\epsilon a_2}-r_{\epsilon }^{+}\le -\frac{a_1}{\epsilon a_2}+\frac{2a_0}{a_1},$

using the lower bound $r_{\epsilon }^{+}\ge -2a_0/a_1$ from Step 1 (so $-r_{\epsilon }^{+}\le 2a_0/a_1$ ). The admissible threshold (2.4) implies

$\frac{2a_0}{a_1}\epsilon \le \frac{a_1}{2a_2}\iff \epsilon \le \frac{a_1^2}{4a_0{a}_2}.$

For such $\epsilon$ ,

$r_{\epsilon }^{-}\le -\frac{a_1}{\epsilon a_2}+\frac{2a_0}{a_1}=-\frac{1}{\epsilon }\left(\frac{a_1}{a_2}-\frac{2a_0\epsilon }{a_1}\right)\le -\frac{1}{\epsilon }\left(\frac{a_1}{a_2}-\frac{a_1}{2a_2}\right)=-\frac{a_1}{2a_2\epsilon }=-\frac{{\kappa }_2}{\epsilon },$

which is the second bound in (4.2).

Step 3 (asymptotic expansions (4.3)-(4.4)). Write $u:=4\epsilon a_2{a}_0/a_1^2\to 0$ and use the binomial (Taylor) expansion $\sqrt{1-u}=1-\frac{1}{2}u-\frac{1}{8}{u}^2+O\left(u^3\right)$ for $\left|u\right|<1$ :

$\begin{array}{c}{\sigma }_{\epsilon }=\sqrt{a_1^2-4\epsilon a_2{a}_0}=a_1\sqrt{1-u}\\ =a_1\left(1-\frac{2\epsilon a_2{a}_0}{a_1^2}-\frac{2{\epsilon }^2{a}_2^2{a}_0^2}{a_1^4}+O\left({\epsilon }^3\right)\right)\\ =a_1-\frac{2\epsilon a_2{a}_0}{a_1}-\frac{2{\epsilon }^2{a}_2^2{a}_0^2}{a_1^3}+O\left({\epsilon }^3\right).\end{array}$

Substituting into $r_{\epsilon }^{+}=\left(-a_1+{\sigma }_{\epsilon }\right)/\left(2\epsilon a_2\right)$ :

$\begin{array}{c}{r}_{\epsilon }^{+}=\frac{-a_1+a_1-\frac{2\epsilon a_2{a}_0}{a_1}-\frac{2{\epsilon }^2{a}_2^2{a}_0^2}{a_1^3}+O\left({\epsilon }^3\right)}{2\epsilon a_2}\\ =\frac{-\frac{2\epsilon a_2{a}_0}{a_1}-\frac{2{\epsilon }^2{a}_2^2{a}_0^2}{a_1^3}+O\left({\epsilon }^3\right)}{2\epsilon a_2}\\ =-\frac{a_0}{a_1}-\epsilon \frac{a_2{a}_0^2}{a_1^3}+O\left({\epsilon }^2\right),\end{array}$

which is (4.3). For $r_{\epsilon }^{-}=\left(-a_1-{\sigma }_{\epsilon }\right)/\left(2\epsilon a_2\right)$ ,

$\begin{array}{c}{r}_{\epsilon }^{-}=\frac{-a_1-a_1+\frac{2\epsilon a_2{a}_0}{a_1}+\frac{2{\epsilon }^2{a}_2^2{a}_0^2}{a_1^3}+O\left({\epsilon }^3\right)}{2\epsilon a_2}\\ =\frac{-2a_1+\frac{2\epsilon a_2{a}_0}{a_1}+\frac{2{\epsilon }^2{a}_2^2{a}_0^2}{a_1^3}+O\left({\epsilon }^3\right)}{2\epsilon a_2}\\ =-\frac{a_1}{\epsilon a_2}+\frac{a_0}{a_1}+\epsilon \frac{a_2{a}_0^2}{a_1^3}+O\left({\epsilon }^2\right),\end{array}$

which is (4.4). (Alternatively, (4.4) follows instantly from (4.3) and Vieta’s identity $r_{\epsilon }^{-}=-a_1/\left(\epsilon a_2\right)-r_{\epsilon }^{+}$ , without a second Taylor expansion.) □

Remark 4.2 Only the exact bounds (4.2), not the asymptotic expansions (4.3)-(4.4), are needed to prove the error estimates of Sections 5-6; the expansions are retained because they give the leading-order constants in the worked example of Section 1.

5. Error Equation and Green’s Function

Define the error $z_{\epsilon }:=y_{\epsilon }-y_0$ . Subtracting (2.5) from (2.1) and using $a_1{y^{\prime }}_0+a_0{y}_0=f$ ,

$\epsilon a_2\left({z^{″}}_{\epsilon }+{y^{″}}_0\right)+a_1\left({z^{\prime }}_{\epsilon }+{y^{\prime }}_0\right)+a_0\left(z_{\epsilon }+y_0\right)=f\Rightarrow \epsilon a_2{z^{″}}_{\epsilon }+a_1{z^{\prime }}_{\epsilon }+a_0{z}_{\epsilon }=-\epsilon a_2{y^{″}}_0.$

Therefore the exact error equation is

$\epsilon a_2{z^{″}}_{\epsilon }+a_1{z^{\prime }}_{\epsilon }+a_0{z}_{\epsilon }=-\epsilon a_2{y^{″}}_0,z_{\epsilon }\left(0\right)=0,{z^{\prime }}_{\epsilon }\left(0\right)=\Delta ,$ (5.1)

where $\text{Δ}:=A_1-{y^{\prime }}_0\left(0\right)=A_1-\left(f\left(0\right)-a_0{A}_0\right)/a_1$ , since $z_{\epsilon }\left(0\right)=y_{\epsilon }\left(0\right)-y_0\left(0\right)=A_0-A_0=0$ and ${z^{\prime }}_{\epsilon }\left(0\right)={y^{\prime }}_{\epsilon }\left(0\right)-{y^{\prime }}_0\left(0\right)=A_1-{y^{\prime }}_0\left(0\right)=\Delta$ .

The Green’s function (impulse-response fundamental solution) of ${ℒ}_{\epsilon }=\epsilon a_2{D}^2+a_{1}D+a_0$ , i.e. the solution of ${ℒ}_{\epsilon }{\mathcal{G}}_{\epsilon }=0$ with ${\mathcal{G}}_{\epsilon }\left(0\right)=0$ , ${{\mathcal{G}}^{\prime }}_{\epsilon }\left(0\right)=1/\left(\epsilon a_2\right)$ , is obtained from the two exponential modes ${\text{e}}^{r_{\epsilon }^{+}t},{\text{e}}^{r_{\epsilon }^{-}t}$ :

${\mathcal{G}}_{\epsilon }\left(t\right):=\frac{{\text{e}}^{r_{\epsilon }^{+}t}-{\text{e}}^{r_{\epsilon }^{-}t}}{\epsilon a_2\left(r_{\epsilon }^{+}-r_{\epsilon }^{-}\right)},t\ge 0.$ (5.2)

Indeed ${\mathcal{G}}_{\epsilon }\left(0\right)=0$ and ${{\mathcal{G}}^{\prime }}_{\epsilon }\left(0\right)=\left(r_{\epsilon }^{+}-r_{\epsilon }^{-}\right)/\left(\epsilon a_2\left(r_{\epsilon }^{+}-r_{\epsilon }^{-}\right)\right)=1/\left(\epsilon a_2\right)$ , as required.

Lemma 5.1 (Uniform Bounds for the Green’s Function and Its Derivative) For every $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ and $t\in \left[0,b\right]$ , estimates (5.3)-(5.5) hold with

$C_G=\frac{2}{a_1},{C^{\prime }}_G=\text{max}\left\{\frac{4a_0}{a_1^2},\frac{2}{a_2}+\frac{4a_0}{a_1^2}\right\},\kappa =\text{min}\left\{\frac{a_0}{a_1},\frac{a_1}{2a_2}\right\}.$

These constants depend only on $a_0,a_1,a_2$ and are independent of $\epsilon$ , $f$ , and the initial data.

$\left|{\mathcal{G}}_{\epsilon }\left(t\right)\right|\le C_G,$ (5.3)

$\left|{{\mathcal{G}}^{\prime }}_{\epsilon }\left(t\right)\right|\le {C^{\prime }}_G\left({\text{e}}^{-\kappa t}+\frac{1}{\epsilon }{\text{e}}^{-\kappa t/\epsilon }\right),$ (5.4)

${\displaystyle {\int }_0^b\left|{\mathcal{G}}_{\epsilon }\left(t\right)\right|\text{d}t}\le C_{G}b,{\displaystyle {\int }_0^b\left|{{\mathcal{G}}^{\prime }}_{\epsilon }\left(t\right)\right|\text{d}t}\le {C^{\prime }}_G\left(\frac{1}{\kappa }+\frac{1}{\kappa }\right).$ (5.5)

Proof. Step 1 (lower bound for the denominator). By (4.1), $\epsilon a_2\left(r_{\epsilon }^{+}-r_{\epsilon }^{-}\right)=\sqrt{a_1^2-4\epsilon a_2{a}_0}={\sigma }_{\epsilon }$ . Definition (2.4) gives $4\epsilon a_2{a}_0\le a_1^2/2$ , hence ${\sigma }_{\epsilon }\ge a_1/\sqrt{2}>a_1/2$ . Equivalently,

$a_1^2-4\epsilon a_2{a}_0\ge \frac{a_1^2}{2}>\frac{a_1^2}{4}.$

Thus both (4.2) and ${\sigma }_{\epsilon }\ge a_1/2$ hold throughout $\left(0,{\epsilon }_{\text{adm}}\right]$ ; no additional parameter threshold is introduced.

Step 2 (bound (5.3)). Since $r_{\epsilon }^{-}<r_{\epsilon }^{+}<0$ by Proposition 4.1 and $t\ge 0$ , both ${\text{e}}^{r_{\epsilon }^{+}t}$ and ${\text{e}}^{r_{\epsilon }^{-}t}$ lie in $\left(0,1\right]$ , so

$\left|{\text{e}}^{r_{\epsilon }^{+}t}-{\text{e}}^{r_{\epsilon }^{-}t}\right|\le \text{max}\left\{{\text{e}}^{r_{\epsilon }^{+}t},{\text{e}}^{r_{\epsilon }^{-}t}\right\}\le 1.$

Hence, using $\epsilon a_2\left(r_{\epsilon }^{+}-r_{\epsilon }^{-}\right)={\sigma }_{\epsilon }\ge a_1/2$ ,

$\left|{\mathcal{G}}_{\epsilon }\left(t\right)\right|=\frac{\left|{\text{e}}^{r_{\epsilon }^{+}t}-{\text{e}}^{r_{\epsilon }^{-}t}\right|}{{\sigma }_{\epsilon }}\le \frac{1}{a_1/2}=\frac{2}{a_1}=:C_G.$

This proves (5.3) with the explicit constant $C_G=2/a_1$ .

Step 3 (bound (5.4)). Differentiating (5.2),

${{\mathcal{G}}^{\prime }}_{\epsilon }\left(t\right)=\frac{r_{\epsilon }^{+}{\text{e}}^{r_{\epsilon }^{+}t}-r_{\epsilon }^{-}{\text{e}}^{r_{\epsilon }^{-}t}}{{\sigma }_{\epsilon }}=\underset{\left(\text{I}\right)}{\underbrace{\frac{r_{\epsilon }^{+}}{{\sigma }_{\epsilon }}{\text{e}}^{r_{\epsilon }^{+}t}}}-\underset{\left(\text{II}\right)}{\underbrace{\frac{r_{\epsilon }^{-}}{{\sigma }_{\epsilon }}{\text{e}}^{r_{\epsilon }^{-}t}}}.$

Term (I). By (4.2), $\left|r_{\epsilon }^{+}\right|\le 2a_0/a_1$ and ${\sigma }_{\epsilon }\ge a_1/2$ , so

$\left|\frac{r_{\epsilon }^{+}}{{\sigma }_{\epsilon }}\right|\le \frac{2a_0/a_1}{a_1/2}=\frac{4a_0}{a_1^2},{\text{e}}^{r_{\epsilon }^{+}t}\le {\text{e}}^{-{\kappa }_{1}t}$ (${\kappa }_1=a_0/a_1$ , by (4.2),

hence $\left|\left(\text{I}\right)\right|\le \frac{4a_0}{a_1^2}{\text{e}}^{-{\kappa }_{1}t}$ .

Term (II). Write $\left|r_{\epsilon }^{-}\right|/{\sigma }_{\epsilon }$ . From Vieta, $r_{\epsilon }^{-}=-a_1/\left(\epsilon a_2\right)-r_{\epsilon }^{+}$ , so $\left|r_{\epsilon }^{-}\right|\le a_1/\left(\epsilon a_2\right)+2a_0/a_1$ ; since ${\sigma }_{\epsilon }\ge a_1/2$ ,

$\left|\frac{r_{\epsilon }^{-}}{{\sigma }_{\epsilon }}\right|\le \frac{2}{a_1}\left(\frac{a_1}{\epsilon a_2}+\frac{2a_0}{a_1}\right)=\frac{2}{\epsilon a_2}+\frac{4a_0}{a_1^2}\le \frac{C_2}{\epsilon }$

for $C_2:=2/a_2+4a_0/a_1^2$ , because $\epsilon \le {\epsilon }_{\text{adm}}\le 1$ . By (0.8), ${\text{e}}^{r_{\epsilon }^{-}t}\le {\text{e}}^{-{\kappa }_{2}t/\epsilon }$ with ${\kappa }_2=a_1/\left(2a_2\right)$ . Hence $\left|\left(\text{II}\right)\right|\le \frac{C_2}{\epsilon }{\text{e}}^{-{\kappa }_{2}t/\epsilon }$ .

Combining (I) and (II) with $\kappa :=\text{min}\left\{{\kappa }_1,{\kappa }_2\right\}$ and ${C^{\prime }}_G:=\text{max}\left\{4a_0/a_1^2,C_2\right\}$ gives

$\left|{{\mathcal{G}}^{\prime }}_{\epsilon }\left(t\right)\right|\le \left|\left(\text{I}\right)\right|+\left|\left(\text{II}\right)\right|\le {C^{\prime }}_G\left({\text{e}}^{-\kappa t}+\frac{1}{\epsilon }{\text{e}}^{-\kappa t/\epsilon }\right),$

which is (5.4).

Step 4 (integral bounds (5.5)). Directly, ${\displaystyle {\int }_0^b\left|{\mathcal{G}}_{\epsilon }\left(t\right)\right|\text{d}t}\le C_G{\displaystyle {\int }_0^b\text{d}t}=C_{G}b$ . For the second,

$\begin{array}{c}{\displaystyle {\int }_0^b\left|{{\mathcal{G}}^{\prime }}_{\epsilon }\left(t\right)\right|\text{d}t}\le {C^{\prime }}_G{\displaystyle {\int }_0^b{\text{e}}^{-\kappa t}\text{d}t}+\frac{{C^{\prime }}_G}{\epsilon }{\displaystyle {\int }_0^b{\text{e}}^{-\kappa t/\epsilon }\text{d}t}\\ ={C^{\prime }}_G\frac{1-{\text{e}}^{-\kappa b}}{\kappa }+{C^{\prime }}_G\frac{1-{\text{e}}^{-\kappa b/\epsilon }}{\kappa }\\ \le \frac{2{C^{\prime }}_G}{\kappa },\end{array}$

using ${\displaystyle {\int }_0^b{\text{e}}^{-\kappa t/\epsilon }\text{d}t}=\epsilon \left(1-{\text{e}}^{-\kappa b/\epsilon }\right)/\kappa \le \epsilon /\kappa$ and bounding both terms by ${C^{\prime }}_G/\kappa$ .

This is (5.5) (with the constant $2{C^{\prime }}_G/\kappa$ in place of the schematic ${C^{\prime }}_G\left(1/\kappa +1/\kappa \right)$ ). □

By the variation-of-parameters formula for the linear operator ${ℒ}_{\epsilon }$ , the unique solution of (5.1) is

$z_{\epsilon }\left(x\right)=\epsilon a_2\Delta {\mathcal{G}}_{\epsilon }\left(x\right)-\epsilon a_2{\displaystyle {\int }_0^x{\mathcal{G}}_{\epsilon }\left(x-s\right){y^{″}}_0\left(s\right)\text{d}s}.$ (5.6)

Remark 5.2 (Derivation of (5.6)) Write $z_{\epsilon }=z_h+z_p$ . The homogeneous part $z_h$ solves ${ℒ}_{\epsilon }{z}_h=0$ with $z_h\left(0\right)=0$ , ${z^{\prime }}_h\left(0\right)=\Delta$ ; since ${\mathcal{G}}_{\epsilon }$ itself already satisfies ${\mathcal{G}}_{\epsilon }\left(0\right)=0$ ,${{\mathcal{G}}^{\prime }}_{\epsilon }\left(0\right)=1/\left(\epsilon a_2\right)$ , the multiple $z_h\left(x\right)=\epsilon a_2\Delta {\mathcal{G}}_{\epsilon }\left(x\right)$ satisfies $z_h\left(0\right)=0$ and ${z^{\prime }}_h\left(0\right)=\epsilon a_2\Delta \cdot \frac{1}{\epsilon a_2}=\Delta$ , as required. The particular part $z_p$ solves ${ℒ}_{\epsilon }{z}_p=-\epsilon a_2{y^{″}}_0$ with zero initial data; by Duhamel’s principle for the impulse response ${\mathcal{G}}_{\epsilon }$ ,

$z_p\left(x\right)={\displaystyle {\int }_0^x{\mathcal{G}}_{\epsilon }\left(x-s\right) \left(-\epsilon a_2{y^{″}}_0\left(s\right)\right)\text{d}s}=-\epsilon a_2{\displaystyle {\int }_0^x{\mathcal{G}}_{\epsilon }\left(x-s\right){y^{″}}_0\left(s\right)\text{d}s},$

because the forcing term of (5.1) is exactly $-\epsilon a_2{y^{″}}_0$ and ${\mathcal{G}}_{\epsilon }$ is, by definition, the response of ${ℒ}_{\epsilon }$ to a unit impulse. Summing $z_h+z_p$ gives (5.6).

6. Uniform Convergence of the Solution

Theorem 6.1 (Optimal $O\left(\epsilon \right)$ Uniform Error Bound) Under (2.3), for every $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ ,

${\Vert y_{\epsilon }-y_0\Vert }_{\infty }\le C\epsilon \left(\left|\Delta \right|+{\Vert {y^{″}}_0\Vert }_{\infty }\right),C:=a_2{C}_G\max \left\{1,b\right\},$ (6.1)

with $C_G=2/a_1$ as in Lemma 5.1. In particular $y_{\epsilon }\to y_0$ uniformly on $\left[0,b\right]$ as $\epsilon \to 0^{+}$ .

Proof. Fix $x\in \left[0,b\right]$ . From the exact representation (5.6) and the triangle inequality,

$\left|z_{\epsilon }\left(x\right)\right|\le \epsilon a_2\left|\Delta \right|\left|{\mathcal{G}}_{\epsilon }\left(x\right)\right|+\epsilon a_2{\displaystyle {\int }_0^x\left|{\mathcal{G}}_{\epsilon }\left(x-s\right)\right|\left|{y^{″}}_0\left(s\right)\right|\text{d}s}.$

By (5.3), $\left|{\mathcal{G}}_{\epsilon }\left(x\right)\right|\le C_G$ for the first term. For the integral term, bound $\left|{y^{″}}_0\left(s\right)\right|\le {\Vert {y^{″}}_0\Vert }_{\infty }$ (constant, taken out of the integral) and $\left|{\mathcal{G}}_{\epsilon }\left(x-s\right)\right|\le C_G$ pointwise for $s\in \left[0,x\right]$ , so

$\begin{array}{c}\epsilon a_2{\displaystyle {\int }_0^x\left|{\mathcal{G}}_{\epsilon }\left(x-s\right)\right|\left|{y^{″}}_0\left(s\right)\right|\text{d}s}\le \epsilon a_2{C}_G{\Vert {y^{″}}_0\Vert }_{\infty }{\displaystyle {\int }_0^{x}1\text{d}s}\\ =\epsilon a_2{C}_G{\Vert {y^{″}}_0\Vert }_{\infty }x\\ \le \epsilon a_2{C}_{G}b{\Vert {y^{″}}_0\Vert }_{\infty },\end{array}$

since $x\le b$ . Adding the two bounds,

$\begin{array}{c}\left|z_{\epsilon }\left(x\right)\right|\le \epsilon a_2{C}_G\left|\Delta \right|+\epsilon a_2{C}_{G}b{\Vert {y^{″}}_0\Vert }_{\infty }\\ =\epsilon a_2{C}_G\left(\left|\Delta \right|+b{\Vert {y^{″}}_0\Vert }_{\infty }\right)\\ \le \epsilon a_2{C}_G\max \left\{1,b\right\}\left(\left|\Delta \right|+{\Vert {y^{″}}_0\Vert }_{\infty }\right).\end{array}$

Since the right-hand side does not depend on $x$ , taking the supremum over $x\in \left[0,b\right]$ gives (6.1) with $C=a_2{C}_G\max \left\{1,b\right\}$ . □

Remark 6.2 Theorem 6.1 strengthens the qualitative convergence statement discussed in Section 2.2 by providing a parameter-uniform $O\left(\epsilon \right)$ bound valid for every $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ .

7. Initial Layer Analysis for the Derivative

Differentiating (5.6) in $x$ (the boundary term at $s=x$ inside the integral is handled by Leibniz’s rule, and vanishes because ${\mathcal{G}}_{\epsilon }\left(0\right)=0$ ):

${z^{\prime }}_{\epsilon }\left(x\right)=\epsilon a_2\Delta {{\mathcal{G}}^{\prime }}_{\epsilon }\left(x\right)-\epsilon a_2{\displaystyle {\int }_0^x{{\mathcal{G}}^{\prime }}_{\epsilon }\left(x-s\right){y^{″}}_0\left(s\right)\text{d}s}.$

Theorem 7.1 (Derivative Estimate with an Initial-Layer Term) For all $x\in \left[0,b\right]$ and $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ , estimate (7.1) holds with

$C=\max \left\{\frac{3}{2},\frac{4a_0{a}_2}{a_1^2},a_2{C}_G\max \left\{1,b\right\}\right\},\kappa =\frac{a_1}{2a_2}.$

Thus $C$ and $\kappa$ depend only on the fixed coefficients and $b$ , and are independent of $\epsilon$ , $f$ , and the initial data.

$\left|{y^{\prime }}_{\epsilon }\left(x\right)-{y^{\prime }}_0\left(x\right)\right|\le C\left|\Delta \right|{\text{e}}^{-\kappa x/\epsilon }+C\epsilon \left(\left|\Delta \right|+{\Vert {y^{″}}_0\Vert }_{\infty }+{\Vert {y^{‴}}_0\Vert }_{\infty }\right).$ (7.1)

Proof. We treat the two terms of ${z^{\prime }}_{\epsilon }\left(x\right)$ separately.

Step 1: the homogeneous (initial-condition) term $\epsilon a_2\Delta {{\mathcal{G}}^{\prime }}_{\epsilon }\left(x\right)$ . From (5.2),

$\epsilon a_2{{\mathcal{G}}^{\prime }}_{\epsilon }\left(x\right)=\epsilon a_2\cdot \frac{r_{\epsilon }^{+}{\text{e}}^{r_{\epsilon }^{+}x}-r_{\epsilon }^{-}{\text{e}}^{r_{\epsilon }^{-}x}}{\epsilon a_2\left(r_{\epsilon }^{+}-r_{\epsilon }^{-}\right)}=\frac{r_{\epsilon }^{+}}{r_{\epsilon }^{+}-r_{\epsilon }^{-}}{\text{e}}^{r_{\epsilon }^{+}x}-\frac{r_{\epsilon }^{-}}{r_{\epsilon }^{+}-r_{\epsilon }^{-}}{\text{e}}^{r_{\epsilon }^{-}x}.$

For the slow component, $r_{\epsilon }^{+}-r_{\epsilon }^{-}={\sigma }_{\epsilon }/\left(\epsilon a_2\right)\ge a_1/\left(2\epsilon a_2\right)$ and $\left|r_{\epsilon }^{+}\right|\le 2a_0/a_1$ . Therefore

$\left|\frac{r_{\epsilon }^{+}}{r_{\epsilon }^{+}-r_{\epsilon }^{-}}\right|\le \frac{2a_0/a_1}{a_1/\left(2\epsilon a_2\right)}=\frac{4a_0{a}_2}{a_1^2}\epsilon ,$

and $\left|{\text{e}}^{r_{\epsilon }^{+}x}\right|\le 1$ . For the fast component, the exact formulae for the roots give

$\left|\frac{r_{\epsilon }^{-}}{r_{\epsilon }^{+}-r_{\epsilon }^{-}}\right|=\frac{a_1+{\sigma }_{\epsilon }}{2{\sigma }_{\epsilon }}\le \frac{3}{2},$

where we used ${\sigma }_{\epsilon }\ge a_1/2$ . Since ${\text{e}}^{r_{\epsilon }^{-}x}\le {\text{e}}^{-{\kappa }_{2}x/\epsilon }$ , the two components satisfy

$\left|\epsilon a_2\Delta {{\mathcal{G}}^{\prime }}_{\epsilon }\left(x\right)\right|\le \frac{3}{2}\left|\Delta \right|{\text{e}}^{-{\kappa }_{2}x/\epsilon }+\frac{4a_0{a}_2}{a_1^2}\epsilon \left|\Delta \right|$ (7.2)

for every $x\in \left[0,b\right]$ and $\epsilon \in \left(0,{\epsilon }_{\text{adm}}\right]$ .

Step 2: the convolution (particular) term. Integrate by parts in

$I\left(x\right):={\displaystyle {\int }_0^x{{\mathcal{G}}^{\prime }}_{\epsilon }\left(x-s\right){y^{″}}_0\left(s\right)\text{d}s}$

by setting $u={y^{″}}_0\left(s\right)$ , $\text{d}v={{\mathcal{G}}^{\prime }}_{\epsilon }\left(x-s\right)\text{d}s$ , so $\text{d}u={y^{‴}}_0\left(s\right)\text{d}s$ and, since $\frac{\text{d}}{\text{d}s}\left[-{\mathcal{G}}_{\epsilon }\left(x-s\right)\right]={{\mathcal{G}}^{\prime }}_{\epsilon }\left(x-s\right)$ , we take $v=-{\mathcal{G}}_{\epsilon }\left(x-s\right)$ :

$I\left(x\right)={\left[-{\mathcal{G}}_{\epsilon }\left(x-s\right){y^{″}}_0\left(s\right)\right]}_{s=0}^{s=x}+{\displaystyle {\int }_0^x{\mathcal{G}}_{\epsilon }\left(x-s\right){y^{‴}}_0\left(s\right)\text{d}s}.$

At $s=x$ : ${\mathcal{G}}_{\epsilon }\left(0\right)=0$ , so this boundary contribution vanishes. At $s=0$ : the term is $-\left(-{\mathcal{G}}_{\epsilon }\left(x\right){y^{″}}_0\left(0\right)\right)={\mathcal{G}}_{\epsilon }\left(x\right){y^{″}}_0\left(0\right)$ (mind the overall sign from evaluating the bracket at the lower limit with a minus sign in front). Hence

$I\left(x\right)={\mathcal{G}}_{\epsilon }\left(x\right){y^{″}}_0\left(0\right)+{\displaystyle {\int }_0^x{\mathcal{G}}_{\epsilon }\left(x-s\right){y^{‴}}_0\left(s\right)\text{d}s}.$

Using $\left|{\mathcal{G}}_{\epsilon }\left(t\right)\right|\le C_G$ from (5.3) for both the boundary term and inside the integral,

$\left|I\left(x\right)\right|\le C_G\left|{y^{″}}_0\left(0\right)\right|+C_G{\displaystyle {\int }_0^x\left|{y^{‴}}_0\left(s\right)\right|\text{d}s}\le C_G{\Vert {y^{″}}_0\Vert }_{\infty }+C_{G}b{\Vert {y^{‴}}_0\Vert }_{\infty }.$

Multiplying by $\epsilon a_2$ ,

$\epsilon a_2\left|I\left(x\right)\right|\le C_5\epsilon \left({\Vert {y^{″}}_0\Vert }_{\infty }+{\Vert {y^{‴}}_0\Vert }_{\infty }\right),C_5:=a_2{C}_G\max \left\{1,b\right\}.$ (7.3)

Step 3: conclusion. Since ${z^{\prime }}_{\epsilon }\left(x\right)=\epsilon a_2\text{Δ}{{\mathcal{G}}^{\prime }}_{\epsilon }\left(x\right)-\epsilon a_{2}I\left(x\right)$ , combining (7.2) and (7.3) via the triangle inequality gives, with $C:=\text{max}\left\{3/2,4a_0{a}_2/a_1^2,C_5\right\}$ and $\kappa :={\kappa }_2$ ,

$\left|{y^{\prime }}_{\epsilon }\left(x\right)-{y^{\prime }}_0\left(x\right)\right|=\left|{z^{\prime }}_{\epsilon }\left(x\right)\right|\le C\left|\Delta \right|{\text{e}}^{-\kappa x/\epsilon }+C\epsilon \left(\left|\Delta \right|+{\Vert {y^{″}}_0\Vert }_{\infty }+{\Vert {y^{‴}}_0\Vert }_{\infty }\right),$

which is (7.1). □

Remark 7.2 Estimate (7.1) separates the possible order-one contribution, whose amplitude is controlled by $\left|\Delta \right|$ , from a remainder of order $\epsilon$ . Once the fast term has decayed, or when $\Delta =0$ , the derivative error is uniformly $O\left(\epsilon \right)$ .

Corollary 7.3 (Uniform Convergence Away from the Layer) For every fixed $\delta \in \left(0,b\right)$ , ${y^{\prime }}_{\epsilon }\to {y^{\prime }}_0$ uniformly on $\left[\delta ,b\right]$ as $\epsilon \to 0^{+}$ ; more precisely,

$\underset{x\in \left[\delta ,b\right]}{\sup }\left|{y^{\prime }}_{\epsilon }\left(x\right)-{y^{\prime }}_0\left(x\right)\right|\le C\left|\Delta \right|{\text{e}}^{-\kappa \delta /\epsilon }+C\epsilon \left(\left|\Delta \right|+{\Vert {y^{″}}_0\Vert }_{\infty }+{\Vert {y^{‴}}_0\Vert }_{\infty }\right)\underset{\epsilon \to 0^{+}}{\to }0.$

Proof. For $x\ge \delta$ , ${\text{e}}^{-\kappa x/\epsilon }\le {\text{e}}^{-\kappa \delta /\epsilon }$ since $t\mapsto {\text{e}}^{-\kappa t/\epsilon }$ is decreasing; substituting this into (7.1) gives the displayed bound. As $\epsilon \to 0^{+}$ , ${\text{e}}^{-\kappa \delta /\epsilon }\to 0$ (in fact faster than any power of $\epsilon$ , since ${\text{e}}^{-\kappa \delta /\epsilon }/{\epsilon }^n\to 0$ for every $n$ , by L’Hôpital or by the standard exponential-dominates-polynomial estimate) and the remainder term $C\epsilon (\cdots )\to 0$ linearly; hence the whole right-hand side tends to 0. □

The Compatibility Condition

Corollary 7.4 (Strong $C^1$ -Convergence under Compatibility) If $\Delta =0$ (equivalently $A_1={y^{\prime }}_0\left(0\right)=\left(f\left(0\right)-a_0{A}_0\right)/a_1$ ), then

${\Vert y_{\epsilon }-y_0\Vert }_{C^1\left(\left[0,b\right]\right)}:={\Vert y_{\epsilon }-y_0\Vert }_{\infty }+{\Vert {y^{\prime }}_{\epsilon }-{y^{\prime }}_0\Vert }_{\infty }\le C\epsilon$ (7.4)

for a constant $C$ depending only on $a_0,a_1,a_2,b,{\Vert {y^{″}}_0\Vert }_{\infty },{\Vert {y^{‴}}_0\Vert }_{\infty }$ .

Proof. If $\text{Δ}=0$ , the order-one term in (7.1) is absent, leaving

$\left|{y^{\prime }}_{\epsilon }\left(x\right)-{y^{\prime }}_0\left(x\right)\right|\le C\epsilon \left({\Vert {y^{″}}_0\Vert }_{\infty }+{\Vert {y^{‴}}_0\Vert }_{\infty }\right),\forall x\in \left[0,b\right].$

Also, setting $\text{Δ}=0$ in (6.1) gives ${\Vert y_{\epsilon }-y_0\Vert }_{\infty }\le C\epsilon {\Vert {y^{″}}_0\Vert }_{\infty }$ . Adding the two bounds (and enlarging $C$ if necessary to absorb both) gives (7.4). □

Remark 7.5 The compatibility condition removes the order-one fast contribution associated with the initial slope mismatch and upgrades uniform convergence of the solution to $C^1$ convergence with the same order $O\left(\epsilon \right)$ .

8. Extensions to Variable Coefficients and Higher-Order Systems

8.1. Variable Coefficients

Consider the prototype

$\epsilon a_2\left(x\right){y^{″}}_{\epsilon }+a_1\left(x\right){y^{\prime }}_{\epsilon }+a_0\left(x\right)y_{\epsilon }=f\left(x\right),y_{\epsilon }\left(0\right)=A_0,{y^{\prime }}_{\epsilon }\left(0\right)=A_1.$ (8.1)

A direct extension of the present argument would assume $a_2,a_1\in C^2\left(\left[0,b\right]\right)$ , $a_0\in C^1\left(\left[0,b\right]\right)$ , $f\in C^2\left(\left[0,b\right]\right)$ , and constants ${\alpha }_2,{\alpha }_1,M_0>0$ such that

$0<{\alpha }_2\le a_2\left(x\right),0<{\alpha }_1\le a_1\left(x\right),{\displaystyle \sum }_{j=0}^2{\Vert a_j\Vert }_{C^1}+{\Vert a_2\Vert }_{C^2}+{\Vert a_1\Vert }_{C^2}\le M_0.$

The reduced problem $a_1\left(x\right){y^{\prime }}_0+a_0\left(x\right)y_0=f$ , $y_0\left(0\right)=A_0$ , must possess a unique $C^3$ solution. In addition, the first-order formulation must admit a slow-fast exponential dichotomy: its evolution operator has a one-dimensional fast component $U_{\epsilon }^f\left(x,s\right)$ satisfying

$\Vert U_{\epsilon }^f\left(x,s\right)\Vert \le M{\text{e}}^{-\kappa \left(x-s\right)/\epsilon },0\le s\le x\le b,$

with $M,\kappa$ independent of $\epsilon$ . Finally, the Green kernel $G_{\epsilon }\left(x,s\right)$ must satisfy analogues of (5.3)-(5.5), uniformly on the triangle $0\le s\le x\le b$ . The natural compatibility condition remains $A_1={y^{\prime }}_0\left(0\right)$ . Under these assumptions the representation and estimates of Sections 5 - 7 can be repeated, with constants depending on $b,{\alpha }_1^{-1},{\alpha }_2^{-1},M_0,M$ , and $\kappa$ , but not on $\epsilon$ . Proving this dichotomy and the kernel bounds is a separate theorem and is not claimed here.

8.2. Higher-Order Equations

For the prototype $\epsilon a_m\left(x\right)y_{\epsilon }^{\left(m\right)}+{\displaystyle {\sum }_{j=0}^{m-1}{a}_j\left(x\right)y_{\epsilon }^{\left(j\right)}}=f$ , assume $a_m$ is bounded away from zero, the coefficients and $f$ have sufficient bounded regularity to produce the derivatives of the reduced solution used by the residual, and the reduced $\left(m-1\right)$ st-order problem is well posed. A usable spectral hypothesis is that the first-order companion system has exactly one fast spectral subspace with real parts bounded above by $-\kappa /\epsilon$ , while its remaining spectrum stays in $\left\{\mathrm{Re}z\le -{\kappa }_s\right\}$ and is separated from the fast spectrum by at least $c/\epsilon$ . The initial datum associated with the lost derivative must match the corresponding derivative of the reduced solution if uniform convergence of that derivative is sought. Without this spectral separation, neither the number nor the scale of the layers is determined by the differential order alone.

8.3. Systems

Consider $\epsilon E\left(x\right){Y^{\prime }}_{\epsilon }\left(x\right)+A\left(x\right)Y_{\epsilon }\left(x\right)=F\left(x\right)$ , where $E$ may lose rank in the limit. Assume $E,A\in C^1\left(\left[0,b\right];{ℝ}^{n\times n}\right)$ , $F\in C^1\left(\left[0,b\right];{ℝ}^n\right)$ , all coefficients and their stated derivatives are bounded, and the reduced differential-algebraic problem is well posed. Let the matrix pencil $\lambda E\left(x\right)+A\left(x\right)$ split uniformly into slow and fast spectral bundles with constant dimensions, uniformly bounded spectral projectors, a gap of size at least $c/\epsilon$ , and fast evolution satisfying $\Vert U_{\epsilon }^f\left(x,s\right)\Vert \le M{\text{e}}^{-\kappa \left(x-s\right)/\epsilon }$ . Compatibility requires the initial vector to have no leading-order component in the fast bundle relative to the reduced initial state. Together with a matrix Green-kernel bound uniform in $\epsilon$ , these assumptions are the system counterpart of those used in Lemma 5.1; see [15]. Their verification for particular systems lies beyond the scope of this paper.

9. Numerical Experiments

We test the explicit constants using the problem of Example 10.1 on $\left[0,1\right]$ . Thus $a_0=a_1=a_2=1$ , $f\equiv 1$ , $A_0=A_1=0$ , $y_0\left(x\right)=1-{\text{e}}^{-x}$ , and $\text{Δ}=-1$ . The exact solution is evaluated from (10.1); no numerical ODE solver or discretization error enters the comparison. The parameter values

$\epsilon \in \left\{2^{-3},2^{-4},2^{-5},2^{-6},2^{-7}\right\}$

all belong to $\left(0,{\epsilon }_{\text{adm}}\right]$ , since ${\epsilon }_{\text{adm}}=1/8$ . Supremum errors were calculated on a uniform grid of ${10}^6+1$ points. For this example, Theorem 6.1 gives

$E_0\left(\epsilon \right):={\Vert y_{\epsilon }-y_0\Vert }_{\infty }\le 4\epsilon ,$ (9.1)

whereas Theorem 7.1 gives the pointwise estimate

$E_1\left(x;\epsilon \right):=\left|{y^{\prime }}_{\epsilon }\left(x\right)-{y^{\prime }}_0\left(x\right)\right|\le 4{\text{e}}^{-x/\left(2\epsilon \right)}+12\epsilon .$ (9.2)

Table 1. Exact maximum-norm errors and explicit theoretical bounds on $\left[0,1\right]$ .

$\epsilon$	$E_0\left(\epsilon \right)$	$E_0\left(\epsilon \right)/\epsilon$	bound $4\epsilon$	${\Vert {y^{\prime }}_{\epsilon }-{y^{\prime }}_0\Vert }_{\infty }$
$1.2500\times {10}^{-1}$	$8.4402\times {10}^{-2}$	0.6752	$5.0000\times {10}^{-1}$	1.0000
$6.2500\times {10}^{-2}$	$4.8402\times {10}^{-2}$	0.7744	$2.5000\times {10}^{-1}$	1.0000
$3.1250\times {10}^{-2}$	$2.6598\times {10}^{-2}$	0.8511	$1.2500\times {10}^{-1}$	1.0000
$1.5625\times {10}^{-2}$	$1.4155\times {10}^{-2}$	0.9059	$6.2500\times {10}^{-2}$	1.0000
$7.8125\times {10}^{-3}$	$7.3643\times {10}^{-3}$	0.9426	$3.1250\times {10}^{-2}$	1.0000

Table 1 shows that $E_0\left(\epsilon \right)/\epsilon$ tends toward a nonzero limit, numerically confirming the sharp $O\left(\epsilon \right)$ rate. The bound $4\epsilon$ holds for every tested value but is deliberately conservative because it is obtained from coefficient-only kernel estimates. The derivative error remains exactly one in the maximum norm because $E_1\left(0;\epsilon \right)=\left|\text{Δ}\right|=1$ for every $\epsilon$ ; this confirms that uniform derivative convergence fails when the compatibility condition is violated.

Figure 1. Exact errors and theoretical bounds for the test problem. Top left: exact solution error on $\left[0,1\right]$ . Top right: $E_0\left(\epsilon \right)$ and the bound $4\epsilon$ on logarithmic axes. Bottom left: derivative error near the initial layer. Bottom right: derivative error and the pointwise bound (0.25), plotted against the stretched coordinate $x/\epsilon$ .

Figure 1 confirms both scales predicted by the analysis. The solution error decreases proportionally to $\epsilon$ , while the derivative discrepancy is concentrated in a region whose width is proportional to $\epsilon$ . In the stretched coordinate $x/\epsilon$ , the derivative profiles nearly collapse onto a common exponentially decaying curve. The explicit bounds remain above the exact errors, thereby validating their parameter-uniform character, while the nonzero initial derivative discrepancy confirms the sharpness of the layer amplitude.

10. Explicit Example and Sharpness

Example 10.1 (Fully Worked Example Validating Sharpness) Take $a_2=a_1=a_0=1$ , $f\equiv 1$ , $A_0=A_1=0$ :

$\epsilon {y^{″}}_{\epsilon }+{y^{\prime }}_{\epsilon }+y_{\epsilon }=1,y_{\epsilon }\left(0\right)=0,{y^{\prime }}_{\epsilon }\left(0\right)=0.$

Reduced problem. ${y^{\prime }}_0+y_0=1$ , $y_0\left(0\right)=0$ ; by (0.6), $y_0\left(x\right)=1-{\text{e}}^{-x}$ , so ${y^{\prime }}_0\left(x\right)={\text{e}}^{-x}$ and ${y^{\prime }}_0\left(0\right)=1$ . Hence $\text{Δ}=A_1-{y^{\prime }}_0\left(0\right)=0-1=-1\ne 0$ : the compatibility condition of Corollary 0.12 fails, and we expect a genuine $O\left(1\right)$ boundary layer in the derivative.

Exact roots. By (4.5) with $a_0=a_1=a_2=1$ : setting $\sigma :={\sigma }_{\epsilon }=\sqrt{1-4\epsilon }$ ,

$r_{\epsilon }^{+}=\frac{-1+\sigma }{2\epsilon },r_{\epsilon }^{-}=\frac{-1-\sigma }{2\epsilon },r_{\epsilon }^{+}-r_{\epsilon }^{-}=\frac{\sigma }{\epsilon }.$

Exact closed form of $y_{\epsilon }$ . Writing $y_{\epsilon }=1+C_1{\text{e}}^{r_{\epsilon }^{+}x}+C_2{\text{e}}^{r_{\epsilon }^{-}x}$ and imposing $y_{\epsilon }\left(0\right)=0$ , ${y^{\prime }}_{\epsilon }\left(0\right)=0$ gives the linear system

$\left(\begin{array}{cc}1& 1\\ r_{\epsilon }^{+}& r_{\epsilon }^{-}\end{array}\right)\left(\begin{array}{c}{C}_1\\ C_2\end{array}\right)=\left(\begin{array}{c}-1\\ 0\end{array}\right),$

whose determinant is $r_{\epsilon }^{-}-r_{\epsilon }^{+}=-\sigma /\epsilon \ne 0$ ; solving by Cramer’s rule,

$\begin{array}{l}{C}_1=\frac{\left(-1\right)r_{\epsilon }^{-}}{-\sigma /\epsilon }=\frac{r_{\epsilon }^{-}\epsilon }{\sigma }=\frac{-1-\sigma }{2\sigma },\\ C_2=\frac{r_{\epsilon }^{+}}{-\sigma /\epsilon }=-\frac{r_{\epsilon }^{+}\epsilon }{\sigma }=\frac{1-\sigma }{2\sigma }.\end{array}$

Thus, exactly (no approximation yet),

$y_{\epsilon }\left(x\right)=1-\frac{1+\sigma }{2\sigma }{\text{e}}^{r_{\epsilon }^{+}x}+\frac{1-\sigma }{2\sigma }{\text{e}}^{r_{\epsilon }^{-}x}.$ (10.1)

At $x=0$ , the expression gives $y_{\epsilon }\left(0\right)=1-\left(1+\sigma \right)/\left(2\sigma \right)+\left(1-\sigma \right)/\left(2\sigma \right)=0$ , as required.

Expansion in $\epsilon$ . Since $\sigma =\sqrt{1-4\epsilon }=1-2\epsilon -2{\epsilon }^2+O\left({\epsilon }^3\right)$ (Step 3 of Proposition 4.1 with $a_0=a_1=a_2=1$ ), we compute the two coefficients to $O\left({\epsilon }^2\right)$ :

$\begin{array}{c}\frac{1+\sigma }{2\sigma }=\frac{2-2\epsilon +O\left({\epsilon }^2\right)}{2\left(1-2\epsilon +O\left({\epsilon }^2\right)\right)}=\frac{1-\epsilon +O\left({\epsilon }^2\right)}{1-2\epsilon +O\left({\epsilon }^2\right)}\\ =\left(1-\epsilon \right)\left(1+2\epsilon \right)+O\left({\epsilon }^2\right)=1+\epsilon +O\left({\epsilon }^2\right),\end{array}$

$\frac{1-\sigma }{2\sigma }=\frac{2\epsilon +O\left({\epsilon }^2\right)}{2\left(1-2\epsilon +O\left({\epsilon }^2\right)\right)}=\epsilon \left(1+2\epsilon +O\left({\epsilon }^2\right)\right)=\epsilon +O\left({\epsilon }^2\right).$

Also $r_{\epsilon }^{+}=\left(-1+\sigma \right)/\left(2\epsilon \right)=\left(-2\epsilon -2{\epsilon }^2+O\left({\epsilon }^3\right)\right)/\left(2\epsilon \right)=-1-\epsilon +O\left({\epsilon }^2\right)$ (consistent with (4.3)), so

${\text{e}}^{r_{\epsilon }^{+}x}={\text{e}}^{-x}{\text{e}}^{-\epsilon x+O\left({\epsilon }^2\right)x}={\text{e}}^{-x}\left(1-\epsilon x+O\left({\epsilon }^2\right)\right).$

Substituting into (10.1) and keeping terms up to $O\left(\epsilon \right)$ :

$\begin{array}{c}\frac{1+\sigma }{2\sigma }{\text{e}}^{r_{\epsilon }^{+}x}=\left(1+\epsilon \right){\text{e}}^{-x}\left(1-\epsilon x\right)+O\left({\epsilon }^2\right)\\ ={\text{e}}^{-x}\left(1+\epsilon -\epsilon x\right)+O\left({\epsilon }^2\right)\\ ={\text{e}}^{-x}+\epsilon \left(1-x\right){\text{e}}^{-x}+O\left({\epsilon }^2\right),\end{array}$

so that

$\begin{array}{c}{y}_{\epsilon }\left(x\right)=1-{\text{e}}^{-x}-\epsilon \left(1-x\right){\text{e}}^{-x}+\epsilon {\text{e}}^{r_{\epsilon }^{-}x}+O\left({\epsilon }^2\right)\\ =y_0\left(x\right)-\epsilon \left(1-x\right){\text{e}}^{-x}+\epsilon {\text{e}}^{r_{\epsilon }^{-}x}+O\left({\epsilon }^2\right).\end{array}$ (10.2)

Since $r_{\epsilon }^{-}\le -{\kappa }_2/\epsilon$ with ${\kappa }_2=1/2$ (Proposition 4.1 with $a_1=a_2=1$ ), the last term satisfies $\epsilon {\text{e}}^{r_{\epsilon }^{-}x}=O\left(\epsilon {\text{e}}^{-x/\left(2\epsilon \right)}\right)$ ; combined with (10.2) this gives

$y_{\epsilon }\left(x\right)-y_0\left(x\right)=\underset{O\left(\epsilon \right),\text{slow correction}}{\underbrace{-\epsilon \left(1-x\right){\text{e}}^{-x}}}+\underset{\text{exponentially localized layer}}{\underbrace{O\left(\epsilon {\text{e}}^{-x/\left(2\epsilon \right)}\right)}}+O\left({\epsilon }^2\right),$

which has the form $O\left(\epsilon \right)+O\left(\epsilon {\text{e}}^{-x/\left(2\epsilon \right)}\right)$ . To verify optimality, fix any $x_{\text{*}}\in \left(0,\text{min}\left\{b,1\right\}\right)$ . Then (10.2) gives

$\frac{y_{\epsilon }\left(x_{\text{*}}\right)-y_0\left(x_{\text{*}}\right)}{\epsilon }\to -\left(1-x_{\text{*}}\right){\text{e}}^{-x_{\text{*}}}\ne 0.$

Hence ${\Vert y_{\epsilon }-y_0\Vert }_{\infty }$ is bounded below by a positive multiple of $\epsilon$ for sufficiently small $\epsilon$ . Together with Theorem 6.1, this proves ${\Vert y_{\epsilon }-y_0\Vert }_{\infty }=\Theta \left(\epsilon \right)$ for this example.

Derivative. Differentiating (10.1) exactly,

$y_{\epsilon \text{'}}\left(x\right)=-\frac{1+\sigma }{2\sigma }{r}_{\epsilon }^{+}{\text{e}}^{r_{\epsilon }^{+}x}+\frac{1-\sigma }{2\sigma }{r}_{\epsilon }^{-}{\text{e}}^{r_{\epsilon }^{-}x}.$

For the first term, using $\frac{1+\sigma }{2\sigma }=1+\epsilon +O\left({\epsilon }^2\right)$ and $r_{\epsilon }^{+}=-1-\epsilon +O\left({\epsilon }^2\right)$ ,

$-\frac{1+\sigma }{2\sigma }{r}_{\epsilon }^{+}=-\left(1+\epsilon \right)\left(-1-\epsilon \right)+O\left({\epsilon }^2\right)={\left(1+\epsilon \right)}^2+O\left({\epsilon }^2\right)=1+2\epsilon +O\left({\epsilon }^2\right),$

so, multiplying by ${\text{e}}^{r_{\epsilon }^{+}x}={\text{e}}^{-x}\left(1-\epsilon x\right)+O\left({\epsilon }^2\right)$ ,

$\begin{array}{c}-\frac{1+\sigma }{2\sigma }{r}_{\epsilon }^{+}{\text{e}}^{r_{\epsilon }^{+}x}=\left(1+2\epsilon \right){\text{e}}^{-x}\left(1-\epsilon x\right)+O\left({\epsilon }^2\right)\\ ={\text{e}}^{-x}\left(1+\epsilon \left(2-x\right)\right)+O\left({\epsilon }^2\right)\\ ={y^{\prime }}_0\left(x\right)+\epsilon \left(2-x\right){\text{e}}^{-x}+O\left({\epsilon }^2\right).\end{array}$

For the second term, $\frac{1-\sigma }{2\sigma }=\epsilon +O\left({\epsilon }^2\right)$ and, since $r_{\epsilon }^{-}=-1/\epsilon +1+O\left(\epsilon \right)$ (by (4.4) with $a_0=a_1=a_2=1$ ), we get $\epsilon r_{\epsilon }^{-}=-1+\epsilon +O\left({\epsilon }^2\right)$ , so

$\frac{1-\sigma }{2\sigma } r_{\epsilon }^{-}{\text{e}}^{r_{\epsilon }^{-}x}=\left(\epsilon r_{\epsilon }^{-}+O\left({\epsilon }^2\right)\right){\text{e}}^{r_{\epsilon }^{-}x}=\left(-1+O\left(\epsilon \right)\right){\text{e}}^{r_{\epsilon }^{-}x}=-{\text{e}}^{r_{\epsilon }^{-}x}+O\left(\epsilon {\text{e}}^{r_{\epsilon }^{-}x}\right).$

Adding the two contributions,

${y^{\prime }}_{\epsilon }\left(x\right)-{y^{\prime }}_0\left(x\right)=-{\text{e}}^{r_{\epsilon }^{-}x}+\epsilon \left(2-x\right){\text{e}}^{-x}+O\left({\epsilon }^2\right)+O\left(\epsilon {\text{e}}^{r_{\epsilon }^{-}x}\right).$ (10.3)

Since $r_{\epsilon }^{-}x=-x/\epsilon +x+O\left(\epsilon x\right)$ , we have ${\text{e}}^{r_{\epsilon }^{-}x}={\text{e}}^{-x/\epsilon }{\text{e}}^{x+O\left(\epsilon x\right)}$ . On the fixed interval $\left[0,b\right]$ , this is bounded by $C_b{\text{e}}^{-\kappa x/\epsilon }$ for any fixed $\kappa <1$ and all sufficiently small $\epsilon$ . Hence (10.3) reads

${y^{\prime }}_{\epsilon }\left(x\right)-{y^{\prime }}_0\left(x\right)=O\left(1\right)\cdot {\text{e}}^{-\kappa x/\epsilon }+O\left(\epsilon \right),$

which is the two-scale structure predicted by Theorem 7.1. Moreover, ${y^{\prime }}_{\epsilon }\left(0\right)-{y^{\prime }}_0\left(0\right)=\text{Δ}=-1$ for every $\epsilon >0$ . Consequently, the derivative cannot converge uniformly on $\left[0,b\right]$ when $\text{Δ}\ne 0$ , and the order-one layer amplitude in the theorem is optimal in general.

11. Conclusions

We derived parameter-uniform estimates for the singularly perturbed Cauchy problem by combining an exact Green’s-function representation with the separation of the slow and fast characteristic roots. The main conclusions are:

• ${\Vert y_{\epsilon }-y_0\Vert }_{\infty }=O\left(\epsilon \right)$ on the complete interval, with a constant independent of $\epsilon$ ;

• the derivative error is bounded by an exponentially localized term proportional to $\left|\text{Δ}\right|$ plus an $O\left(\epsilon \right)$ remainder;

• if $\text{Δ}=0$ , then $y_{\epsilon }\to y_0$ in $C^1\left(\left[0,b\right]\right)$ with order $O\left(\epsilon \right)$ ;

• the closed-form example proves that the uniform solution rate and the order-one derivative discrepancy at the initial point are optimal in general.

The discussion of variable coefficients, higher-order equations, and systems also clarifies which kernel estimates and spectral assumptions would be needed to extend the argument. These extensions require separate analysis and are not claimed by the present results.

Possible continuations include nonlinear perturbations and numerical schemes whose error bounds remain uniform as $\epsilon \to 0^{+}$ .

Acknowledgments

The author thanks the institutions of affiliation for their support. Special thanks to Dr. Lorgio Félix Batard Martínez for his foundational work on singular perturbations, and to Aisseli Comet Marrero for her detailed error estimates, which provided the starting point for our rigorous analysis.

Author Contributions

All authors contributed to the mathematical analysis and the revision of the manuscript. All authors read and approved the final version.