Numerical Method for Solving n Singularly Perturbed Mixed-Type Initial Value Problems with the Same Perturbation Parameter and Discontinuous Source Terms

1. Introduction

Singularly perturbed initial–boundary value problems occur frequently across applied mathematics and engineering disciplines. For example, systems of first-order, singularly perturbed ordinary differential equations are fundamental in chemical reactor modeling. The small magnitude of the perturbation parameters induces a multi-scale behavior that typically precludes closed-form solutions. Consequently, the exact solution often features sharply varying layers—initial, boundary, or interior—confined to narrow regions. Over the past several decades, a variety of numerical schemes that achieve uniform convergence with respect to these small parameters have been developed and analyzed.

Examine a Robin-type, singularly perturbed initial-value system with discontinuous sources over $\Omega =(0,1]$, assume the source term is discontinuous at a single point d within the domain $\Omega$. Let ${\Omega }^{-}=(0,d)$ and ${\Omega }^{+}=(d,1]$. The jump of a function $\omega$ at the point d is defined by

$$[\omega ](d)=\omega (d+)-\omega (d-)$$

where $\omega (d+)$ and $\omega (d-)$ denote the right and left limits of $\omega$ at d, respectively. The problem can be stated as finding the solution to the initial value problem $u_1,u_2,\dots ,u_n\in \mathbb{D}={\mathbb{C}}^0(\overline{\Omega })\cap {\mathbb{C}}^1({\Omega }^{-}\cup {\Omega }^{+})$, such that

$$\overrightarrow{L}\overrightarrow{u}(x)=E{\overrightarrow{u}}^{\prime }(x)+A(x)\overrightarrow{u}(x)=\overrightarrow{f}(x),x\in {\Omega }^{-}\cup {\Omega }^{+}$$

(1)

along with the initial conditions provided

$$\overrightarrow{\beta }\overrightarrow{u}(0)=\overrightarrow{u}(0)-\epsilon {\overrightarrow{u}}^{\prime }(0)=\overrightarrow{\varphi }$$

(2)

where, $E=diag(\epsilon ,\epsilon ,\dots ,\epsilon ),\overrightarrow{u}(x)={(u_1(x),u_2(x),\dots ,u_n(x))}^T,A(x)={(a_{ij}(x))}_{n\times n}$ and $\overrightarrow{f}(x)={(f_i(x))}_{n\times 1}$.

From (1) and (2), the problem can equivalently be expressed in operator form as follows:

$$\overrightarrow{L}\overrightarrow{u}=\overrightarrow{f}\mathrm{on}\Omega$$

(3)

with

$$\overrightarrow{\beta }\overrightarrow{u}(0)=\overrightarrow{\varphi }$$

(4)

where the definitions of the operators $\overrightarrow{L}and\overrightarrow{\beta }$ are given by

$$\overrightarrow{L}=ED+A,\overrightarrow{\beta }=I-ED.$$

Here, the operator I is the identity, and $D={\displaystyle \frac{d}{dx}}$ signifies the first derivative operator.

Assumption 1.

The functions $a_{ij},f_i\in {\mathbb{C}}^{(2)}(\overline{\Omega }),i,j=1,\dots ,n$ satisfy the following positivity conditions

$$\left.\begin{array}{cccc}\hfill & (i)\hfill & \hfill a_{ii}(x)& >\sum _{\begin{array}{c}j\ne i\\ j=1\end{array}}^n|a_{ij}(x)|\mathrm{for}i=1(1)n\hfill \\ \hfill & (ii)\hfill & \hfill a_{ij}(x)& \le 0\mathrm{for}i\ne j\mathrm{and}i=1(1)n\hfill \end{array}\right\}\forall x\in \overline{\Omega }.$$

(5)

Assumption 2.

The positive value α adheres to the inequality

$${\displaystyle 0<\alpha <\underset{\begin{array}{c}i=1(1)n\\ x\in \overline{\Omega }\end{array}}{min}\left\{\sum _{j=1}^n{a}_{ij}(x)\right\}.}$$

(6)

Assumption 3.

The parameter ε, representing the singular perturbation and distinct, is assumed to satisfy $0<\epsilon \le 1$.

The given problem is singularly perturbed in the sense that setting $\epsilon =0$ in system (1) yields the following reduced linear algebraic system.

$$A(x)\overrightarrow{v}(x)=\overrightarrow{f}(x),x\in {\Omega }^{-}\cup {\Omega }^{+}$$

(7)

where $A(x)={(a_{ij}(x))}_{n\times n},$$\overrightarrow{v}(x)={(v_1(x),\dots ,v_n(x))}^T$ and $\overrightarrow{f}(x)={(f_1(x),\dots ,f_n(x))}^T$.

The source terms $f_1(x),f_2(x),\dots ,f_n(x)$ are sufficiently smooth throughout $\overline{\Omega }$ except at the point d. The component of the solution $u_1,u_2,\dots ,u_n$ of problem (1) and (2) display initial layers that overlap at $x=0$ and interior layers that overlap just to the right of the discontinuity located at $x=d.$

For a detailed account of parameter-uniform numerical techniques used in the study of singular perturbation problems, see [3,4,6]. The seminal work of Shishkin [8] laid the groundwork for addressing singularly perturbed reaction-diffusion equations with discontinuous coefficients. In the context of scalar equations, Dunne and Riordan [2] investigated initial value problems exhibiting singular perturbations and discontinuous data. Numerical methods that maintain robustness with respect to the perturbation parameter for systems of such problems were discussed in [7]. In [12] a parameter-uniform numerical method was developed for systems where all singular perturbation parameters are equal, leading to solution components with initial layers of uniform width and thus simplifying the analysis. The case in which the singular perturbation parameter appears in only one of the equations was examined in [5]. The most general and challenging scenario, involving individual initial layers for each solution component that overlap and influence one another, has been studied in [1,9,13]. Notably, [9] introduced a uniformly convergent numerical method of first order accuracy (up to a logarithmic factor), while [1] proposed a hybrid finite difference scheme on a piecewise-uniform Shishkin mesh achieving nearly second-order accuracy uniformly with respect to both small parameters. In all these previous studies, the source term was assumed to be smooth. In contrast, the present work deals with discontinuous source terms, resulting in each solution component exhibiting both initial and interior layers.

Theorem 1.

Let $A(x)$ satisfy (5) and (6). The problem (1) - (2) admits a solution $\overrightarrow{u}\in \mathbb{D}$.

Proof. 

Examine $\overrightarrow{y}$ and $\overrightarrow{z}$ as the specific solutions to the differential equations

$$\epsilon y_i^{\prime }(x)+\sum _{j=1}^n{a}_{ij}(x)y_j(x)=f_i(x),i=1,2,\dots ,n,\mathrm{for}\mathrm{all}x\in {\Omega }^{-}$$

(8)

and

$$\epsilon z_i^{\prime }(x)+\sum _{j=1}^n{a}_{ij}(x)z_j(x)=f_i(x),i=1,2,\dots ,n,\mathrm{for}\mathrm{all}x\in {\Omega }^{+}.$$

(9)

Let us analyze the function

$$\begin{array}{c}\hfill u_i(x)=\left\{\begin{array}{c}{y}_i(x)+\beta (u_i(0)-y_i(0)){\varphi }_i(x),i=1,2,\dots ,n,x\in {\Omega }^{-}\hfill \\ z_i(x)+B_i{\varphi }_i(x),i=1,2,\dots ,n,x\in {\Omega }^{+}\hfill \end{array}\right.\end{array}$$

(10)

where ${\varphi }_i$ is the solution of

$$\begin{array}{c}\hfill \left.\begin{array}{c}\hfill \epsilon {\varphi }_i^{\prime }+\sum _{j=1}^n{a}_{ij}(x){\varphi }_j(x)=0\\ \hfill {\beta }_i{\varphi }_i(0)=1,\end{array}\right\}i=1,2,\dots ,n,\mathrm{for}\mathrm{all}x\in \Omega .\end{array}$$

Here $B_i,i=1(1)n$ is chosen so that $\overrightarrow{u}\in \mathbb{D}$. In $\Omega ,0<\overrightarrow{\varphi }\le 1,$ The function $\overrightarrow{\varphi }$ cannot attain an internal maximum or minimum, and therefore, ${\varphi }_i^{\prime }<0,i=1(1)n$ in $\Omega$. Choose the constants $B_i$ such that

$$\begin{array}{c}\hfill \overrightarrow{y}(d-)=\overrightarrow{z}(d+),\overrightarrow{u}(d-)=\overrightarrow{u}(d+).\end{array}$$

The existence of the constants $B_i$ requires that

$$\begin{array}{c}\hfill {\displaystyle \frac{[u_i(0)-y_i(0)]{\varphi }_i(d-)}{{\varphi }_i(d+)}}\ne 0\mathrm{for}i=1(1)n.\end{array}$$

Since ${\varphi }_i(d+)>0$, the existence of $\overrightarrow{B}$ and consequently $\overrightarrow{u}$ is guaranteed. □

Remark: In this section, C stands for a general vector of positive constants, which remain unaffected by the perturbation parameters and the discretization parameter N.

2. Mathematical Analysis

The operator $\overrightarrow{L}$ adheres to the following maximum principle.

Lemma 1.

Let $A(x)$ satisfy conditions (5) and (6). Suppose that a function $\overrightarrow{u}\in \mathbb{D}$ satisfies $\overrightarrow{\beta }\overrightarrow{u}(0)\ge \overrightarrow{0},[\overrightarrow{u}](d)=\overrightarrow{0},[\epsilon u_i^{\prime }](d)\le 0,$ and $\overrightarrow{L}\overrightarrow{u}(x)\ge \overrightarrow{0}$ for all $x\in {\Omega }^{-}\cup {\Omega }^{+}$. Then, $\overrightarrow{u}(x)\ge \overrightarrow{0}$ for all $x\in \overline{\Omega }$.

Proof. 

Assume that $A(x)$ meets the requirements specified in (5) and (6). Let $\overrightarrow{u}\in \mathbb{D}$ be a function such that

$$\overrightarrow{\beta }\overrightarrow{u}(0)\ge \overrightarrow{0},[\overrightarrow{u}](d)=\overrightarrow{0},[\epsilon u_i^{\prime }](d)\le 0,$$

and

$$\overrightarrow{L}\overrightarrow{u}(x)\ge \overrightarrow{0}\mathrm{for}\mathrm{every}x\in {\Omega }^{-}\cup {\Omega }^{+}.$$

Then it follows that

$$\overrightarrow{u}(x)\ge \overrightarrow{0}\mathrm{for}\mathrm{all}x\in \overline{\Omega }.$$

Case (i):  $p_1\in {\Omega }^{-}\cup {\Omega }^{+}$

$$\begin{array}{cc}\hfill \overrightarrow{\beta }\overrightarrow{u}(0)& =\overrightarrow{u}(0)-\epsilon {\overrightarrow{u}}^{\prime }(0)\hfill \\ \hfill & <0,\mathrm{a}\mathrm{contradiction}\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill {(\overrightarrow{L}\overrightarrow{u})}_1(p_1)& =\epsilon u_1^{\prime }(p_1)+\sum _{j=1}^n{a}_{1j}(p_1)u_j(p_1)\hfill \\ \hfill & =\epsilon u_1^{\prime }(p_1)+\sum _{j=1}^n{a}_{1j}(p_1)u_j(p_1)+\sum _{j=2}^n{a}_{1j}(p_1)u_1(p_1)-\sum _{j=2}^n{a}_{1j}(p_1)u_1(p_1)\hfill \\ \hfill & <0,\mathrm{thus}\mathrm{arriving}\mathrm{at}\mathrm{a}\mathrm{contradiction}.\hfill \end{array}$$

Case (ii):   $p_1=d$

Given that $\overrightarrow{u}\in \mathbb{C}(\Omega )$ and $u_1(d)<0$, there exists a neighborhood $N_h=(d-h,d+h)$ around d where $u_1(x)<0$ holds for every $x\in N_h$. Select a point $x_1\in N_h$, distinct from d, such that $u_1(x_1)>u_1(d)$. By the Mean Value Theorem, there exists some $x_2\in N_h$ for which

$$u_1^{\prime }(x_2)={\displaystyle \frac{u_1(d)-u_1(x_1)}{d-x_1}}<0.$$

Since $x_2$ lies within $N_h$, an argument analogous to that used in the first case implies that,

$$\begin{array}{cc}\hfill {(\overrightarrow{L}\overrightarrow{u})}_1(x_2)& =\epsilon u_1^{\prime }(x_2)+\sum _{j=1}^n{a}_{1j}(x_2)u_j(x_2)<0\hfill \end{array}$$

thus arriving at a contradiction. □

As a straightforward implication of the previous lemma, the following stability result holds.

Lemma 2.

Let $A(x)$ satisfies (5) and (6). Let $\overrightarrow{u}$ be the solution of (1) and (2). Then,

$$||\overrightarrow{u}(x){||}_{\overline{\Omega }}\le max\left\{||\overrightarrow{\beta }\overrightarrow{u}(0)||,{\displaystyle \frac{1}{\alpha }}||\overrightarrow{L}\overrightarrow{u}{||}_{{\Omega }^{-}\cup {\Omega }^{+}}\right\}.$$

Proof. 

Introduce the pair of functions

$$\begin{array}{cc}\hfill {\overrightarrow{\theta }}^{\pm }(x)& {\displaystyle =max\left\{\Vert \overrightarrow{\beta }\overrightarrow{u}(0)\Vert ,{\displaystyle \frac{1}{\alpha }}{\Vert \overrightarrow{L}\overrightarrow{u}\Vert }_{{\Omega }^{-}\cup {\Omega }^{+}}\right\}\pm \overrightarrow{u}(x),x\in \overline{\Omega }}\hfill \\ \hfill {\overrightarrow{\theta }}^{\pm }(x)& =M\pm \overrightarrow{u}(x)\hfill \end{array}$$

where $M=max\{||\beta \overrightarrow{u}(0)||,{\displaystyle \frac{1}{\alpha }}||\overrightarrow{L}\overrightarrow{u}|{|}_{{\Omega }^{-}\cup {\Omega }^{+}}\}$. Then, it is true that $\overrightarrow{\beta }{\overrightarrow{\theta }}^{\pm }(0)\ge \overrightarrow{0},[{\overrightarrow{\theta }}^{\pm }](d)\le \overrightarrow{0}$, by a proper choice of C and $\overrightarrow{L}{\overrightarrow{\theta }}^{\pm }(x)\ge \overrightarrow{0}$ on ${\Omega }^{-}\cup {\Omega }^{+}$. It follows from Lemma 1 that ${\overrightarrow{\theta }}^{\pm }(x)\ge \overrightarrow{0}$ on $\overline{\Omega }.$ Hence,

$$|\overrightarrow{u}(x)|\le max\left\{||\overrightarrow{\beta }\overrightarrow{u}(0)||,{\displaystyle \frac{1}{\alpha }}||\overrightarrow{L}\overrightarrow{u}{||}_{{\Omega }^{-}\cup {\Omega }^{+}}\right\}.$$

□

Lemma 3.

Assume that $A(x)$ fulfills the conditions specified in (5) and (6). Let $\overrightarrow{u}$ denote the solution of (1) and (2). Then, for each index $i=1,2,\dots ,n$ and for all $x\in {\Omega }^{-}\cup {\Omega }^{+}$, there exists a constant C such that

$$\begin{array}{cc}\hfill |u_i(x)|& \le C\left\{\Vert \overrightarrow{\varphi }\Vert +\Vert \overrightarrow{f}{\Vert }_{{\Omega }^{-}\cup {\Omega }^{+}}\right\}\hfill \\ \hfill |u_i^{\prime }(x)|& \le C{\epsilon }^{-1}\left\{\Vert \overrightarrow{\varphi }\Vert +\Vert \overrightarrow{f}{\Vert }_{{\Omega }^{-}\cup {\Omega }^{+}}\right\}\hfill \\ \hfill |u_i^{\prime \prime }(x)|& \le C{\epsilon }^{-2}\left\{\Vert \overrightarrow{\varphi }\Vert +\Vert \overrightarrow{f}{\Vert }_{{\Omega }^{-}\cup {\Omega }^{+}}+{\Vert \overrightarrow{f^{\prime }}\Vert }_{{\Omega }^{-}\cup {\Omega }^{+}}\right\}.\hfill \end{array}$$

Proof. 

The argument proceeds analogously to the proof of Lemma 2 in [11].