Higher-order convergence analysis for interior and boundary layers in a semi-linear reaction-diffusion system networked by a $k$-star graph with non-smooth source terms

1.   Introduction

Differential equations on connected graphs are an emerging research domain with several applications in physical, chemical, and biological models. In physical processes, such as lateral oscillations within a network of strings, beam vibrations connected by a network, and the stationary states of electrons within a molecule, the modeling has been done typically through boundary value problems on each edge of the network ^[1]. Pavlov et al. ^[2] gives a concrete example, where the free electrons encompass the interaction between molecules and their surrounding medium by a second-order boundary value problem connected by a graph. These problems are connected by a continuity condition at the junction nodes, which imposes the equilibrium of forces acting on each node from the adjacent edges. In ^[3], Nagatani formulates a traffic flow model on a star graph using a nonlinear diffusion equation. He examines traffic flow management by employing cell transmission on the graphs. The concept of discontinuity in traffic flow on each edge of the graph arises because the continuous flow may not be possible through the connecting vertices of the star, cycle, and complete graphs. It is important to note that a graph $\mathcal{G}$ is a finite set of vertices $V(\mathcal{G})$ connected by a finite set of edges $E(\mathcal{G})$ with a relation between each edge and corresponding two vertices (endpoints) (see ^[4]). More concrete applications on differential equations connected by a graph appear inside chemical models ^[5] where a neural network model predicts the chemical reactivity. One can also look into reaction-diffusion models in ecological and chemical contexts ^[6] and the blood flow model inside vessels in one-dimensional flow ^[7,8].

Since graph theoretical models are different than usual problems, there are limited attempts to consider the existence of solutions for problems connected by the star graphs. In ^[9], the authors find the non-constant solution of the reaction-diffusion problem under the continuity and Kirchhoff conditions at the junction vertices. After that, Iwasaki et al. investigated the stability and instability of wave solutions of the reaction-diffusion problems on a metric graph in ^[10]. A theoretical approach to nonlinear fractional differential equations is also presented in ^[11] to demonstrate the well-posedness with Ulam Hyers stability. Existence and uniqueness results for fractional differential equations on a $k$-star graph domain are explored in ^[12]. The numerical approximation of differential equations connected by star graphs has received attention from a few authors. Gordezian et al. ^[13] considered a second-order linear differential equation on each edge of the $k$-star graph (see Figure 1) and used the Dirichlet boundary conditions at the outer vertices and the continuity preserving Kirchhoff's condition at the junction vertex. Their numerical analysis uses a central difference scheme for second-order derivatives and a forward difference scheme for first-order derivatives to achieve linear accuracy across the domain.

Singularly perturbed differential equations are crucial due to their extensive applications in hydrodynamics and gas dynamics. The concepts of boundary and interior layers in these areas were introduced in the nineteenth century. Ongoing efforts aim to understand the boundary and interior layer simulations associated with various singularly perturbed differential equations, including the viscous Burgers' and Navier–Stokes equations. For a recent survey that references singularly perturbed boundary value problems, one can refer to Gie ^[14]. In the present analysis, we focus on the effect of unbounded solution derivatives due to the presence of the viscosity parameter and show its complicacy to obtain the approximate solutions in maximum norm due to the presence of the viscosity parameter.

In ^[15], the authors provide a higher-order error analysis for the boundary layer-oriented singularly perturbed reaction-diffusion coupled systems with mixed boundary conditions on an equidistributed mesh. The nonlinear singularly perturbed problems and their coupled version are also considered for optimal convergence analysis in ^{[16,17,18,19,20]} on the equidistributed meshes. Note that the computational cost gradually increases with the number of equations in coupled systems. Therefore, the reduction of computational cost is also an important feature, which is carried out in ^[21,22,23] for reaction-dominated problems having different diffusion parameters. However, the generation of this mesh is relatively more complicated than the piecewise uniform meshes used in this article. Analogously, the problems on the $k$-star graph are also relatively complex and behave like a coupled system since edges are connected at a junction point. We refer to ^[24,25,26] for the well-posedness of the nonlinear hyperbolic conservation law on the network. Numerical approximation based on finite volume discretization for the hyperbolic conservation problem has been developed in ^[24,27]. Moreover, validation of these schemes has been done in ^[28] by utilizing the exact solution of the Riemann problem.

Furthermore, the linear transport problem on the network has been considered for convergence analysis of the analytical solution based on parabolic approximation by utilizing the suitable transmission conditions in the inner node ^[29,30]. For the numerical approximation, we refer to ^[31,32], where the authors proposed a hybrid discontinuous Galerkin method to find the numerical solution on the network and provide the detailed analysis for the error estimate. Here, we are interested in the nonlinear singularly perturbed elliptic problem with non-smooth data on the star-shaped network for the convergence analysis of the higher-order scheme along with strong numerical evidence. However, some works are available on studying singular perturbation problems on the $k$-star graphs with smooth data. Kumar et al. ^[33] discuss a linear singularly perturbed reaction-diffusion problem connected by a star graph, with Dirichlet boundary conditions at the outer vertices, and the continuity and Kirchhoff's conditions at the junction point. This work utilizes a central difference scheme for the first and second-order derivatives over a piecewise uniform Shishkin mesh to obtain an almost second-order convergence. They consider the M-matrix criteria to derive the maximum principle based on the discretization of the graph Laplacian. However, it is unclear from the numerical point of view. Hence, there is a research gap on nonlinear singularly perturbed boundary value problems with non-smooth data, particularly in producing higher-order accurate uniformly stable solutions. In this regard, we mention the seminal monograph ^[34] having non-smooth source functions, which considers a linear singularly perturbed reaction-diffusion problem with Dirichlet boundary conditions on the fitted mesh to obtain parameter uniform linear accurate solutions. Moreover, a higher-order accuracy was obtained for convection-dominated problems having discontinuous convection terms to get an almost second-order accuracy based on a five-point scheme at the point of discontinuity ^[35]. One can also see more higher-order schemes ^[36,37] to obtain uniformly convergent solutions in space. The increasing need for real-life issues such as flow rate optimization ^[7,8], reaction-dominated processes ^[33], and convection-dominated processes ^[38] drive the study of elliptic problems on networks. Elliptic partial differential equations (PDEs) describe equilibrium states where the system is balanced and no temporal changes occur. While studying these equations on networks, the solutions must account for the network's specific structure and satisfy the equilibrium conditions at each node to ensure continuity and compatibility across the network's junction points. These complexities are not typically encountered while studying scalar elliptic problems defined on simple domains ^[34]. Therefore, we study a more general class of elliptic problems with non-smooth source terms, which are typically applicable in real-world scenarios.

Based on the above motivations, we observe that nonlinear steady-state singularly perturbed reaction-diffusion problems having discontinuous source terms with mixed boundary conditions are more general and realistic from applied sciences on a star graph domain. Our major goals are the higher-order accurate results for this class of problems and their mathematical analysis with numerical difficulties.

The structure of this paper is as follows: In Section 2, we state a nonlinear reaction-diffusion problem with non-smooth data on a $k$-star graph. In Section 3, we analyze the properties of the analytic solution and provide its derivative bounds. In Section 4, we present a discretization using a three-point scheme and a five-point scheme on the Shishkin mesh on the $k$-star graph. We also estimate the maximum error within the graph domain. In Section 5, we focus on the numerical experiments conducted on the $3$-star graph, investigate its errors, and analyze the convergence of experimental results of the test problem.

2.   Problem description

A $k$-star graph $\mathcal{G}$ is a collection of a finite set of vertices $V(\mathcal{G}) = \{v_i; \; i = 0, 1, \cdots, k\},$ and edges $E(\mathcal{G}) = \{e_i; \; i = 1, \cdots, k\}$, where all the $k$ vertices $v_i's\; \text{for}\; i = 1, \cdots, k$ (see Figure 1) have degree (no. of edges incident on that vertex) $1$, and one central vertex (see, $v_0$ in Figure 1) has degree $k$ ($v_0$ is incident with $k$ edges). Here $v_0$ is called the junction point, and $e_i$ defines the edge joining the vertex $v_i$ to the junction vertex $v_0$, i.e., $e_i = \overrightarrow{v_iv_0}$. We also consider $\ell_i = |e_i| = |\overrightarrow{v_iv_0}|,$ for all $i = 1, \cdots, k$. Then, the singularly perturbed boundary layer-originated problems are defined on each edge of the $k$-star graph $\mathcal{G}$ and connected at the junction point.

Hence, without loss of generality, we assume each outer vertex $v_{i}$, for $i = 1, \cdots, k,$ as an origin. We can identify each edge $e_i$ as $e_{i}: = (0, \ell_i)$. Then, the singularly perturbed semi-linear reaction-diffusion problem with discontinuous source term as given on the $k$-star graph, where we represent $u$ along the $i$th edge $u_{i}$, is defined by:

where $T_{\varepsilon_{i}}$ is a nonlinear operator, defined by $T_{\varepsilon_{i}}u_{i}(x): = -\varepsilon_{i}u^{\prime\prime}_{i}(x)+f_{i}(x, u_{i}(x))$, $T_{0_{i}}$ is the linear operator at the free boundary vertices, defined by $T_{0_{i}}u_{i}(0): = -\mu_{i, 1}u_{i}(0)+\mu_{i, 2}u^{\prime}_{i}(0)$ and $T_{\ell_i}\equiv I$ is the identity operator, defined by $T_{\ell_i}u_{i}(\ell_i): = u_{i}(\ell_i)$ and, $u^{J}$ is unknown where the superscript $J$ indicates the junction point, $\Omega^{-}_{i} = (0, d_{i}), \; \Omega^{+}_{i} = (d_{i}, \ell_i)$ and $\Omega_{i} = \Omega^{-}_{i} \cup \Omega^{+}_{i}$. We assume that source term $g_{i}(x)$ is discontinuous at the point $x = d_{i}$, i.e., $g_{i}(d_{i}-)\neq g_{i}(d_{i}+),$ where $d_{i}\in e_{i} = (0, \ell_i).$ We also assume that $\psi_{i}$, $\mu_{i, 1}$ are positive constants and $\mu_{i, 2}$ is either a constant or takes the value $\varepsilon_{i}$. The jump of any function $\delta$ at $d_{i}$ is denoted by $[\delta](d_{i}) = \delta(d_{i}+)-\delta(d_{i}-).$ For all $x \in \bar{\Omega}_{i},$ we assume that the semi-linear term $f_{i}(x, u_{i}(x))$ is sufficiently smooth, $f_i(x, 0) = 0$ and satisfy the following condition

Note that $u_{i}(x) \notin C^2(e_{i})$, as $g_{i}$ is discontinuous at $d_{i}$, and $\, u_{i}(d_{i}+) = u_{i}(d_{i}-).$

The reduced problem corresponding to Eqs (2.1)–(2.4) is defined by

Since the source term $g_i(x)$ is discontinuous at the point $x = d_{i}$, we expect an interior layer of width $O(\sqrt{\varepsilon})$ near the discontinuity point ^[34].

Notation: In this paper, $C$ denotes a generic positive constant independent of the diffusion parameters $\varepsilon_{i}$ and the number of step sizes $N_{i}$. Here, $\frac{\partial f}{\partial u}(x, .)$ denotes the partial derivative concerning the dependent variable. $\| \cdot \|_{\Omega^N}$ denotes the maximum norm over a given mesh function on the domain $\Omega^N$. We use the symbol $C^{p}(\Omega)$ to represent the space of $p$-th time continuously differentiable functions on the domain $\Omega$.

3.   Analytical properties

In this section, we discuss the analytical behavior of the solution on each edge of the $k$-star graph domain. Let us assume $e_{i} = (0, \ell_i) = (0, 1),$ on every edge of the $k$-star graph, for all $i = 1, 2, \cdots, k.$

Theorem 1. On each edge of the $k$-star graph domain, let $u_{i}$ be the solution of the $i$th edge for $1\leq i\leq k$. Then $u_{i} \in C^{0}(\bar e_{i}) \cap C^{1}(e_{i}) \cap C^{2}(\Omega^{-}_{i} \cup\Omega^{+}_{i}).$

Proof. Let $u_{i}$ be the solution of the $i$th edge, and $u_{i, 1}$ and $u_{i, 2}$ be the particular solutions of the following nonlinear differential equations

where $u_{i, 1}\in C^2(\Omega^{-}_{i})$ and $u_{i, 2} \in C^2(\Omega^{+}_{i})$.

Consider the function

where $\Upsilon_{i, 1}$ and $\Upsilon_{i, 2}$ are defined as follows

Here, we choose the constants $A$ and $B$ in such a way that $u_{i}\in C^{1}(e_{i})$ and $u_{i}\in C^{2}(\Omega^{-}_{i} \cup\Omega^{+}_{i})$. Note that the solutions of Eqs (3.2) and (3.3) on $(0, 1)$, $0 < \Upsilon_{i, 1}, \Upsilon_{i, 2} < 1,$. Thus, we cannot have maximum and minimum in the interior points, and hence $\Upsilon^{\prime}_{i, 1} < 0, \; \; \Upsilon^{\prime}_{i, 2} > 0$ on the interval $(0, 1).$

Now, we choose the constants $A$ and $B$ by imposing the following condition, which ensures the continuous differentiability at all the interior points, i.e., $u_{i}\in C^{1}(e_{i})$. Based on this, we can write

From the above conditions, we get

For the existence of $A$ and $B$, we need

This gives $\Upsilon^{\prime}_{i, 2}(d_{i}) \Upsilon_{i, 1}(d_{i})-\Upsilon^{\prime}_{i, 1}(d_{i})\Upsilon_{i, 2}(d_{i})\neq 0.$ Hence, $u_{i}\in C^{0}(\bar e_{i}) \cap C^{1}(e_{i}) \cap C^{2}(\Omega^{-}_{i} \cup\Omega^{+}_{i})$ exists.

The well-posedness of the problem is considered in the Hilbert space framework. From Eq (2.1) along with the homogeneous boundary conditions (2.2)–(2.4) (by assuming $\psi_i = 0$ and $u^J = 0$) on $i$th edge $e_i$, we define the pivotal space of the entire graph $\mathcal{G}$ as

Now, we define the Sobolev space on $e_i$:

and, the Sobolev space on the graph $\mathcal{G}$

Define the corresponding Sobolev norm ^[30] as follows:

Similarly,

Analogously, let us define the space on the edge $e_i$ as $\mathbb{V}_i^{*} = \Bigl\{u_i\in H^{1}_{e_i}; T_{0_i}u_i(0) = 0\Bigl\}$, and introduce the dense space $\mathbb{D}^{*}(T_{\varepsilon_{i}}): = \Bigl\{u_i\in \mathbb{V}_i^{*}; T_{\varepsilon_{i}}u_{i}\in L^2(e_i)\Bigl\}$. Let us also define the spaces of the entire graph as $\mathbb{V}^{*}: = \Bigl\{u\in H^1(\mathcal{G}); T_{\ell_i}u_i(\ell_i) = T_{\ell_k}u_k(\ell_k), \; \text{for all} \; i, k, \; \text{and} \; T_{0_i}u_i(0) = 0\Bigl\}$ and define the dense space on the entire graph $\mathcal{G}$ is $\mathbb{D}^{*}(L): = \Bigl\{u\in \mathbb{V}^{*}; u\in \prod_{\substack{i = 1}} ^{{ k}}\mathbb{D}^{*}(T_{\varepsilon_{i}}), \; \sum_{\substack{i = 1}}^{{k}}\varepsilon_i u_i^{\prime}(\ell_i) = 0\Bigl\}$ on which we define the operator $L$ as

Consider that $g\in{H(\mathcal{G})}$ such that $\left.g\right|_{e_i} = g_i.$ Multiplying the state equation (2.1) by an arbitrary test function $v\in \mathbb{D}^{*}(L)$ and integrating by parts, we get the weak formulation of the problems (2.1)–(2.4):

Here, the quadratic form $B(\cdot, \cdot): \mathbb{D}^{*}(L)\times \mathbb{D}^{*}(L)\longrightarrow \mathbb{R}$ is defined by

and

Now, we introduce the minimization problem corresponding to the weak formulation (3.4)

where

is the "energy functional" and $F_i(x, v_i) = \int_{0}^{v_i}f_i(x, s)ds$.

Lemma 1. The energy function $\mathbb{E}(\cdot)$ is coercive, continuous, strictly convex, Fréchet differentiable and it satisfies

Proof. The result follows from ^[39,40]. The rigorous procedure can be seen in ^[41].

Theorem 2. Assume that $g\in{H(\mathcal{G})}$. Then, the weak formulation (3.4) and the minimization problem (3.7) are equivalent, and both admit a unique solution.

Proof. The results follow from Lemma 1 by considering the analysis given in ^[40,41].

Moreover, considering the velocity vector zero in the convection-diffusion-reaction nonlinear problem given in ^[42], the stability of the solution of the current nonlinear reaction-diffusion problem follows for $f_{i}(x, u_i(x)) = k_i(u_i)u_i,$ where $k_i(u_i)$ satisfies the Lipschitz condition with respect to $u_i$. Note that we consider the test function $v = 0$ at the boundary.

3.1.   Bounds of solution and its derivatives

For singularly perturbed reaction-diffusion problems, it is well known that the solution can have boundary and interior layers in its domain of definition depending on the smoothness of the data and boundary conditions. In these layer regions, the solution varies abruptly and has large, steep gradients. In contrast, outside these regions, the solution varies smoothly and has bounded derivatives up to some finite order. A rigorous convergence analysis of numerical approximations; requires derivatives of the solution.

Hence, to simplify this convergence analysis, we decompose the solution into the regular component $r_i$ and the singular component $s_i$ such that $u_i = r_i+s_i,$ on each $i$th edge of the graph. The regular component $r_i,$ typically represents the smooth, well-behaved part of the solution that varies gradually across the domain and has bounded derivatives up to some finite order. This part is often associated with the domain's bulk or interior regions with minimal perturbation parameter effects. The convergence analysis for this component becomes more straightforward as the required derivatives for this analysis are bounded. The singular component $s_i$ captures the rapid variations and steep gradients near boundaries or interfaces. In reaction-diffusion problems, these regions often exhibit behaviors not captured by the regular part alone, such as boundary or interior layers where the solution changes abruptly. The singular component $s_i$ corrects the discrepancies between the regular part and the actual solution, especially near regions where the perturbation parameter has a significant impact. The analysis of this component is complicated, in general.

Let us assume $\mu_{i, 2}$ is a constant. Therefore, we decompose the solution into two parts on each $i$th edge of the $k$-star graph domain, such that

where $u_{i}(x)$ is a solution of the $i$th edge of the $k$-star graph domain.

Then, the regular component $r_{i}(x)$ satisfies the following problem

We further decompose the regular component solution $r_i$ of Eq (3.9) into two different parts of the domain, so that the analysis available for problems having two boundary layers can be utilized in both the part of the domains separately. This decomposition is done as follows $r_i = v_{i, 1}+v_{i, 2}$, where $v_{i, 1}$ is the solution of the following problem

and $v_{i, 2}$ is the solution of the following problem

Similarly, we can derive a similar expression from Eq (3.10) on the domain $\Omega^{+}_i$.

Note that the singular component $s_{i}(x)$ satisfies the following problem

Furthermore, we decompose the singular component $s_{i}(x)$, as follows

i.e., $s_{L_{i}}(x)$ is defined on $\Omega^{-}_{i}$ and $s_{R_{i}}(x)$ is defined on $\Omega^{+}_{i}.$ Then, from Theorem 1, singular component $s_{i}$ is well defined; as follows

where $\phi_{i, m}\; \text{for} \; m = 1, 2, 3, 4,$ satisfy the following boundary value problems

and $\beta_{1}$ and $\beta_{2}$ are constants that we need to choose in such a way that the jump condition at $x = d_{i}$ is satisfied, that is,

and

From above, we can write

where $b^1_{i}(x) = \int_{0}^{1} \frac{\partial f_{i}}{\partial u_{i}}(x, r_{i}+q s_{L_{i}}\, )dq$, and

where $b^2_{i}(x) = \int_{0}^{1} \frac{\partial f_{i}}{\partial u_{i}}(x, r_{i}+t s_{R_{i}}\, ) dt$.

Let us use the following notations for the layer functions

Theorem 3. On the $i$th edge of the $k$-star graph domain, the regular component satisfies the following derivatives bounds

where $C$ is independent of $\varepsilon_{i}.$

Proof. We can prove this theorem step by step. First, we establish the derivative bounds for the regular component on the domain $\Omega_i^-$ using Eqs (3.9), (3.11) and (3.12). We employ a similar argument as in Theorem 1; given in ^[43] to establish the bounds for the regular components and their derivatives. In a similar way, we can establish the derivative bounds of the regular components on $\Omega_i^+$ using an analogous decomposition of the regular component from Eq (3.10) on the domain $\Omega_i^+$. Finally, we obtain the bounds of the regular component and its derivatives by adding them over $e_i$.

Theorem 4. On the $i$th edge of the $k$-star graph domain, the singular component satisfies the following bounds

where $C$ is independent of $\varepsilon_{i}.$

Proof. The linear problems (3.15) and (3.16) together with the boundary conditions in Eqs (3.13) and (3.14) are similar to the problems considered in ^[43]. Hence, we obtain the bounds for $s_{L_{i}}$ and $s_{R_{i}}$ and their derivatives analogously.

Remark 1. Consider that $\mu_{i, 2} = \varepsilon_{i}$. In this case, the derivative bounds of the components will be sharp; compared to $\mu_{i, 2}$ as a constant. Therefore, the bounds for Dirichlet type problems; will be used for further analysis, which is as follows

and

4.   Discretization and convergence analysis

Assume that $e_{i} = (0, 1)$ which is obtained by setting $\ell_{i} = 1,$ for all $i = 1, 2, \cdots, k$, as shown in Figure 1. We discretize each edge using a fitted mesh approach with $N_{i}+1$ number of mesh points to discretize the domain $\bar\Omega_{i} = \bar\Omega^{-}_{i}\cup\bar\Omega^{+}_{i}.$

Now, partition the domain $\bar\Omega^{-}_{i} = \left[0, d_{i}\right]$ into following three sub-intervals

and the domain $\bar\Omega^{+}_{i} = [d_{i}, 1]$ into following three sub-intervals

where

On the sub-intervals $\left[\tau_{i, 1}, d_{i}-\tau_{i, 1}\right]$ and $\left[d_{i}+\tau_{i, 2}, 1-\tau_{i, 2}\right],$ a uniform mesh with $N_{i}/4$ mesh-intervals are placed, whereas on the sub-intervals $\left[0, \tau_{i, 1}\right], \left[d_{i}-\tau_{i, 1}, d_{i}\right]$, $\left[d_{i}, d_{i}+\tau_{i, 2}\right]$ and $\left[1-\tau_{i, 2}, 1\right]$, a uniform mesh with $N_{i}/8$ mesh-intervals are placed. So, we define the discrete domain as follows:

Now consider $\Omega^{N_{i}}_{i} = \Omega_{i, 1}^{N_{i}}\, \cup\, \Omega_{i, 2}^{N_{i}}, \, \mbox{where} \, \Omega_{i, 1}^{N_{i}} = \left\{x_{i, j}\right\}_{j = 1}^{\frac{N_{i}}{2}-1} \mbox {, } {\Omega}_{i, 2}^{N_{i}} = \left\{x_{i, j}\right\}_{j = \frac{N_{i}}{2}+1}^{N_{i}-1},$ $x_{i, 0} = 0$, $x_{i, N_{i}} = 1$ and $\Omega^{N_{i}}_{i*} = \Omega^{N_{i}}_{i}\cup \{x_{i, N_{i}}\}.$ It is clear that $x_{i, \frac{N_{i}}{2}} = d_{i}$. Therefore, we consider the whole discretize domain for $i$th edge denoted by $\bar {\mathbb{D}}_{h_i}^{N_{i}} = \left\{x_{i, j}\right\}_{j = 0}^{N_{i}}$, the step size is defined by $h_{i, j} = x_{i, j+1}-x_{i, j}.$

4.1.   Discrete problem

Consider that $U_{i, j} = U_{i}(x_{i, j})$ be the numerical solution on the $i$th edge of the $k$-star graph domain on the fitted mesh. The following expressions offer the discrete formulation corresponding to Eqs (2.1)–(2.4):

where $i, m = 1, 2, \cdots, k, \delta^2U_{i, j} = \dfrac{2}{h_{i, j+1}+h_{i, j}}\left(\frac{U_{i, j+1}-U_{i, j}}{x_{i, j+1}-x_{i, j}}-\frac{U_{i, j}-U_{i, j-1}}{x_{i, j}-x_{i, j-1}}\right),$ $D^{+}_{*}U_{i, j} = \frac{-U_{i, j+2}+4U_{i, j+1}-3U_{i, j}}{2h_{i, j+1}}$ (three-point approximation at right hand side) and $D^{-}_{*}U_{i, j} = \frac{U_{i, j-2}-4U_{i, j-1}+3U_{i, j}}{2h_{i, j}}$ (three point approximation at left hand side). Note that; the order of convergence for Kirchhoff's condition at the junction point; is also almost second-order by using the derivative bound from Remark 1, i.e.,

Consider that $\mathcal{T}^{N_{i}}_{\varepsilon_i}$ is the Fréchet derivative ^[43] of $T^{N_{i}}_{\varepsilon_i}$. Hence, it satisfies the following linear problem:

where $a_i(x_{i, j}) = \frac{\partial f_{i}}{\partial u_{i}}\left(x_{i, j}, \xi_j\right)$ is provided through the mean value theorem and $\mathcal{T}^{N_{i}}_{\varepsilon_i}$ is a linear operator. Note that we use the three-point discretization at the left and right boundary points and the point of discontinuity (i.e., $x_{i, N_{i}/2} = d_{i}$).

At the point $x_{i, N_{i}/2} = d_{i}$, we use the difference operator $T^{N_{i}}_{d_{i}}$ from ^[44], i.e.,

where $h_{i, N_{i}/2+1},$ and $h_{i, N_{i}/2}$ are the right and left step sizes of the point $x_{N_{i}/2}$ and $h_{i, N_{i}/2+1} = 8\tau_{i, 2}/N_{i}$ and $h_{i, N_{i}/2} = 8\tau_{i, 1}/N_{i}.$

Similarly, for the left boundary point, we employ the right-hand side three-point approximation

and at the right boundary point, we discretize Kirchhoff's condition using the following three-point approximation

Therefore, we get the following:

By Eqs (4.13) and (4.14) and the continuity condition (4.10), the matrix corresponding to Eqs (4.12)–(4.14) is not an M-matrix.} To assure the invertibility of the discrete system through the M-matrix criterion, we must transform Eqs (4.12)–(4.14) using Eq (4.7) as follows:

For simplicity, we consider that $h_{i, N_{i}/2+1} = h_{i, N_{i}/2} = h_{i}$. After using the values $U_{i, N_{i}/2+2}$ and $U_{i, N_{i}/2-2}$ in Eq (4.12), we obtain

Consider that $h_{i, N_{i}} = h_{i, 1} = H_{i}$ and from Eqs (4.18) and (4.13), we get

and from Eq (4.19) and Kirchhoff's condition (4.14), we get

Now, we define the discrete operator as follows:

and

Then, we have the following system of equations:

with the transformed Kirchoff condition

Lemma 2. Assume that $N_{i}^2/\ln^2 N_{i}\geq 128\frac{\delta^{*}_{i}}{\delta_{i}},$ where $\delta^{*}_{i} = \max_{x\in e_{i}}{a_{i}(x)}$ and $\delta_{i} = \min_{x\in e_{i}}{a_{i}(x)}.$ Then, the operator $\mathcal{T}^{N_{i}}_{M, \varepsilon_i}$ defined in Eq (4.23), with the boundary conditions (4.24)–(4.26), will lead to an M-matrix.

Proof. By Eqs (4.20)–(4.22) and (2.5), it is clear that $\frac{a_{i, j}h^2_{i}}{\varepsilon_{i}} > 0, \; \; {\mbox{for}}\; \; j = N_{i}/2-1, N_{i}/2+1$ and $\frac{a_{i, j}H^2_{i}}{\varepsilon_{i}} > 0 \; \; {\mbox{for}}\; \; j = 1, N_{i}-1$. Hence, the discretization produces an M-matrix, provided $\left(2-\frac{a_{i, j}h^2_{i}}{\varepsilon_{i}}\right)\geq 0,$ and $\left(2-\frac{a_{i, j}H^2_{i}}{\varepsilon_{i}}\right)\geq 0.$ We use a central difference scheme except for the boundary and discontinuous points. Therefore, we get

Thus, the operator $\mathcal{T}^{N_{i}}_{M, \varepsilon_i}$ satisfies the discrete maximum principle on each edge of the $k$-star graph. Therefore, by consequences of the discrete maximum principle, the discrete solution is stable in the $k$-star graph domain $\mathcal{G}$. Then, the graph Laplacian for the discrete graph produces an M-matrix and satisfies the discrete maximum principle ^[33,45]. Thus, Eqs (4.23)–(4.26) lead to a unique solution as describes an invertible M-matrix.

Remark 2. Here, we are specifying the conditions for individual edges of the $k$-star graph (Figure 1). The boundary conditions for the outer vertices will remain unchanged. However, the continuity and Kirchhoff's conditions at the junction point will not apply to single edges. For instance, we consider a Dirichlet-type boundary condition at the junction point, where $u^J$ will be treated as a known value at the boundary (junction point).

4.1.1.   Discrete stability result

Lemma 3. On each edge $e_{i}$ with $1\leq i\leq k$, let us denote $\Psi_{i, j}$ as a mesh function such that $\Psi_{i, 0} = 0,$ and $\Psi_{i, N_i} = u^{J}$. Then,

Proof. Let $\Phi_{i, j}^{\pm}$ be mesh functions defined as $\Phi_{i, j}^{\pm} = \frac{x_{i, j}}{\delta_{i} d_{i}} \max_{{1\leq j\leq N_{i}-1}}|\mathcal{T}^{N_{i}}_{M, \varepsilon_i}\Psi_{i, j}|+|u^{J}|\pm \Psi_{i, j},$ $0\le j\le N_{i}/2$ and $\Phi_{i, j}^{\pm} = \frac{(1-x_{i, j})}{\delta_{i}(1- d_{i})} \max_{{1\leq j\leq N_{i}-1}}|\mathcal{T}^{N_{i}}_{M, \varepsilon_i}\Psi_{i, j}|+|u^{J}| \pm \Psi_{i, j},$ $N_{i}/2 < j\le N_{i}.$ It is clear that $\Phi^{\pm}_{i, 0} > 0$ and $\Phi^{\pm}_{i, N_{i}} > 0$. Again, $\mathcal{T}^{N_{i}}_{M, \varepsilon_i}\Phi^{\pm}_{i, j}\geq 0,$ for $1\leq j\leq N_{i}-1.$

Hence, by using the discrete maximum principle ^[33], we get

4.2.   Convergence analysis

In this section, we provide the convergence result for each edge of the $k$-star graph domain.

Theorem 5. On each edge $e_{i}$ of the $k$-star graph domain, let $u_{i, 0}$ be the solution of the problem (2.1) at the left boundary point, and $U_{i, 0}$ be the solution of the corresponding discrete problem (4.7) at the left boundary point. Then,

where $C$ is independent of $\varepsilon_{i}$ and $N_{i}$.

Proof. At the left boundary point $x_{i, 0} = 0,$ we have

Hence, we get the following by using Lemma 3

The case $\mu_{i, 2} = \varepsilon_{i}$ leads to the same result based on a suitable choice of a suitable barrier function. Note that the derivatives are bound by $C\varepsilon_{i}^{-3/2}$ for this case (see Remark 1). However, the convergence proof will be simpler as $\mu_{i, 2} = \varepsilon_{i}$.

Theorem 6. On each edge $e_{i}$ of the $k$-star graph domain, let $u_{i}$ be the solution of the problem (2.1), and $U_{i}$ be the solution of the corresponding discrete problem (4.7). Then,

where $C$ is independent of $\varepsilon_{i}$ and $N_{i}$.

Proof. We consider the following two cases to prove this theorem.

Case Ⅰ: Consider $\mu_{i, 2}$ as a constant. Here, we decompose the discrete solution into a regular component $R_{i}(x_{i, j})$ and a singular component $S_{i}(x_{i, j})$, such that

where $R_{i}(x_{i, j})$ is the solution of the non-homogeneous problem and $S_{i}(x_{i, j})$ is the solution of the homogeneous problem ^[46]. Correspondingly, the error $(U_{i}-u_{i})(x_{i, j})$ can be decomposed as

From the original problem and the discretized problem, we get

Now, we find the error estimate separately. Let us first derive it for the regular component. Note

Here, we use the piecewise uniform mesh and its property $h_{i, j+1}+h_{i, j}\leq 2N^{-1}$, and have applied Theorem 4. Now, let us introduce the function $\Theta_{i, j} = C_iN^{-2}_{i}+C_iN^{-2}_{i}\tau_{i, 1} \left\{\begin{array}{l}\frac{x_{i, j}}{\tau_{i, 1}}, \; \; \quad 0 < x_{i, j}\leq \tau_{i, 1}, \\ 1, \; \; \; \; \; \; \; \quad \tau_{i, 1}\leq x_{i, j}\leq d_{i}-\tau_{i, 1}, \\ \frac{d_{i}-x_{i, j}}{\tau_{i, 1}}, \quad d_{i}-\tau_{i, 1}\leq x_{i, j} < d_{i}, \end{array}\right.$on the interval $(0, d_{i})$ and introduce the function $\Theta_{i, j} = C_iN^{-2}_{i}+C_iN^{-2}_{i}\tau_{i, 2} \left\{\begin{array}{l}\frac{x_{i, j}-d_{i}}{\tau_{i, 2}}, \; \; \quad d_{i} < x_{i, j}\leq d_{i}+\tau_{i, 2}, \\ 1, \; \; \; \; \; \; \; \; \quad d_{i}+\tau_{i, 2}\leq x_{i, j}\leq 1-\tau_{i, 2}, \\ \frac{1-x_{i, j}}{\tau_{i, 2}}, \quad 1-\tau_{i, 2}\leq x_{i, j}\leq 1, \end{array}\right.$ on the interval $(d_{i}, 1].$ Note that $\Theta_{i, j}\leq C_i\sqrt{\varepsilon_{i}}N^{-2}_{i}\ln N_{i}$ and

Now, consider the barrier function

By the discrete maximum principle ^[45], we obtain

Now, we decompose the singular component $S_{i}$ along the interval $(0, d)$ and $(d, 1];$ as follows

where $S_{L_i}$ is defined on $(0, d);$ and $S_{R_i}$ is defined on $(d, 1];$ and are given by

Note that either $\tau_{i, 1} = \frac{d_{i}}{4},$ and $\tau_{i, 2} = \frac{1-d_{i}}{4}$, or $\tau_{i, 1} = \tau_{i, 2} = 2\sqrt{\frac{\varepsilon_{i}}{\delta_{i}}} \ln N_{i} \le \frac{d_{i}}{4}$. Thur, if $x_{i, j} = d_{i} = \frac{1}{2}$, we have $\tau_{i, 1} = \tau_{i, 2}$. For $\tau_{i, 1} = \tau_{i, 2} = \frac{1}{8},$ the mesh will be uniform throughout the domain $\Omega_{i*}^{N_{i}}$.

Hence, by using the argument in ^[47] for the singular component $S_{{L_i}}$ on $\Omega^{N_{i}}_{i, 1},$ we get

Similarly, for the singular component $S_{R_{i}}$ on $\Omega^{N_{i}}_{i, 2}\cup \{x_{i, N_{i}}\},$ we get

Combining Eqs (4.28) and (4.29), we get

Hence, by utilizing Lemma 3, from Eqs (4.27) and (4.30), we obtain

Case Ⅱ: Here, we assume $\mu_{i, 2} = \varepsilon_{i}$, we determine the error estimate in a similar way as given in ^[33].

Theorem 7. On each edge $e_{i}$ of the $k$-star graph, let $u_{i, N_{i}/2}$ be the solution of the problem (2.1) at the discontinuity point, and $U_{i, N_{i}/2}$ be the solution of the corresponding discrete problem (4.7) at the discontinuity point. Then,

where $C$ is independent of $\varepsilon_{i}$ and $N_{i}$.

Proof. Let us take the case that $\mu_{i, 2}$'s are positive constants. Then, we have the following consistency estimate at the point of discontinuity $x_{i, N_i/2} = d_{i}$:

Hence, by applying Theorem 4 with $h_{i} = 8\sqrt{\frac{\varepsilon_{i}}{\delta_{i}}} N^{-1}_{i} \ln N_{i}$, we get

For case $\mu_{i, 2} = \varepsilon_{i},$ let us consider the barrier function $\eta_{d_{i}},$ defined by

From the discrete maximum principle on each intervals $[0, {d_i}]$ and $[{d_i}, 1]$, one can easily get

and

Define the ancillary continuous functions $w_{i, 1}$, $w_{i, 2}$ by

In addition, define

Hence, from ^[45], we get the following estimate

Therefore, for the intervals $[0, {d_i}]$ and $[{d_i}, 1]$ respectively, we have the following estimates

For $j = N_{i}/2$,

For sufficiently large $N_{i}$, let us consider the function

Note that, for $i = N_{i}/2$

Hence, for sufficiently large $N_{i}$, we have

Hence, we get the required result from Eqs (4.31) and (4.33).

Theorem 8. On each edge $e_{i}$ of the $k$-star graph domain, let $u_{i}$ be the solution of the problems (2.1) and (2.2), with the Dirichlet boundary condition (2.3) where $u^J$ is known at the junction point, and $U_{i}$ be the solution of the corresponding discrete problems (4.7)–(4.9) with the Dirichlet boundary condition (4.10) where $u^J$ is known at the junction point. Then

where $C$ is independent of $\varepsilon_{i}$ and $N_{i}$.

Proof. Utilizing the Theorems 5–7, the desired result follows.

Remark 3. The present numerical scheme provides an almost second-order parameter uniform convergence on piecewise uniform meshes; as depicted in Theorems 5–8. This order of accuracy is sharper than the numerical approximations, where we use upwind discretizations for mixed-type boundary conditions and the Kirchhoff law at the junction point. In addition, the logarithmic term inside the order of accuracy $O(N^{-2}_{i}\ln^2 N_{i})$ cannot be removed if the continuous domain is replaced by a piecewise uniform mesh ^[48]. The simple structure of this mesh assumes the boundary and interior locations and their widths in advance. If this information is unavailable, one can use other kinds of adaptive moving meshes or equidistributed meshes, like in ^[17,20,21] to enhance the rate of accuracy. However, the adaptive mesh generation algorithms are complicated and these algorithms also increase the computational costs ^[22] for coupled systems, which can be greatly reduced by splitting of the coupled matrices ^[23], appearing in the original model.

Now, we present the main result of this paper on the $k$-star graph, which gives the almost second-order $\varepsilon_{i}$-uniform estimate.

Theorem 9. On the $k$-star graph domain, let $u_{\varepsilon}$ be the solution of the problems (2.1)–(2.4), and $U_{\varepsilon}$ be the solution of the corresponding discrete problems (4.7)–(4.11). Then,

where $N_i = N$, over all edges; and $C$ is independent of $\varepsilon$ and $N.$

Proof. We have proved in Theorem 8 that an almost second-order error estimate holds throughout every edge of the $k$-star graph when $u^J$ is known. In the present $k$-star graph domain, all the edges are connected at a single point on the boundary, known as the junction point, where the $u^J$ is unknown. Thus, we need to solve the solution components over all edges simultaneously, instead of edge-wise, as $u^J$ is unknown. Note that the continuity condition (2.3) for the continuous problem; is equivalent to the continuity condition (4.10) at the discrete problem. In addition, to match the number of equations with the number of discrete variables, we need another equation, which will be given by the Kirchhoff equation (4.11). Hence, the number of variables of discrete problems matches the number of equations for the discrete problem along the $k$ star graph domain. We have proved that the rate of accuracy along every edge is preserved up to $N_i^{-2}\ln^2 N_i$ in Theorem 8. This rate is also preserved for the present discretization of the Kirchhoff equation, see Eq (4.6). In addition, the discrete problems (4.23)–(4.26) is uniformly stable, as provided in Lemma 3. Hence, we can write the following:

Now, let $U_{\varepsilon}$ and $u_{\varepsilon}$ be the discrete and analytical solutions on the $k$-star graph, respectively. Then, the error estimate on the $k$-star graph is given by:

where $N^{-2} \ln^2 N = \max(N^{-2}_{i} \ln^2 N_{i})$ for all $i = 1, 2, \ldots, k$, and $C$ is a constant independent of $\varepsilon$ and $N$, and $\bar{\mathbb{D}}^{N}_h$ represents the $k$-star graph domain.

Hence, the maximum global error on the entire $k$-star graph domain is independent of $\varepsilon$ and bounded by $C N^{-2} \ln^2 N$.

5.   Numerical experiments

In this section, we discuss the numerical experiments to verify the theoretical results, obtained in the earlier sections. Here, we consider the perturbation parameters from the set $\varepsilon = (\varepsilon_1, \varepsilon_2, \varepsilon_3)$ and use the same number of mesh intervals $N$ for $k = 3$ star graph in Example 1. We take the singular perturbation parameters from the set $\mathbf{S} = \{2^{-(p-1)}; \; p = 1, 2, \cdots, 40\},$ and discretize the problem having the transition parameters $\tau_{1}, \tau_2, \; \text{and}\; \tau_{3}$ on the $3$-star graph. We compute the nodal errors of the approximate solution using the double mesh principle ^[17,21] as follows:

where $U^{N, \varepsilon_i}_{i, j}$ represents the approximate solution components at $x_{i, j}$ across the $i$th edge of the $3$-star graph with $N+1$ mesh points for each edge, and $U^{2N, \varepsilon_i}_{i, j}$ represents the solution with $2N+1$ mesh points, obtained by bisecting the previous mesh points on each edge. The corresponding order of convergence is calculated as follows:

Example 1. Let us consider the following nonlinear reaction-diffusion problem with discontinuous source terms on the $3$-star graph over $[0, 1]$ (see Figure 2):

with any one of the boundary conditions at the free vertices

The continuity condition at the junction point is

where $u^{J}$ is unknown (here, $J$ stands for junction point), and Kirchhoff's condition at the junction point is

We use Newton's quasi-linearization technique to linearize the given nonlinear problem. Let us first approximate Kirchhoff's junction law with the standard upwind scheme. At discontinuity points, we employ the scheme (4.8), and for Neumann conditions, a first-order upwind approximation of the derivative is utilized. Tables 1–3 depict the uniform errors of overall solution components and the corresponding linear accuracy of the numerical solution. Hence, it illustrates that the upwind scheme yields an almost first-order uniformly convergent approximate solution of the current Example 1.

Moreover, we apply our proposed scheme, which utilizes a three-point scheme (4.9) for the first-order derivatives in Kirchhoff's junction law (4.11) and mixed boundary conditions. At the points of discontinuities, we employ a five-point scheme (4.8) for higher-order accuracy. We consider three different types of boundary conditions for Example 1: Dirichlet boundary conditions (5.4a), mixed boundary conditions (5.4b), and scaled mixed boundary conditions (5.4c). By applying this scheme, we obtain a higher-order accurate solution for each boundary condition. To justify this, we find the uniform maximum error on each edge of the $3$-star graph domain. Then, we take the maximum of these errors over all edges of the $3$-star graph domain to compute the uniform errors. Let us first discuss the Example 1 with the boundary conditions (5.4a): Table 4 shows that the maximum error estimate on the single edge $e_1$ for problem (5.1) with the boundary condition given in (5.4a), where $\varepsilon \in \mathbf{S}$. Similarly, Table 5 presents the maximum error for problem (5.2) on the edge $e_2$, with the boundary conditions (5.4a), where we take the perturbation parameters from the set $\varepsilon \in \mathbf{S}$. Analogously, Table 6 shows the maximum error for problem (5.3) on the edge $e_3$, with the boundary conditions (5.4a).

Note that the nonlinear problems (5.1)–(5.3) are coupled at the junction point. Hence, we impose the continuity condition (5.5) and Kirchhoff's condition (5.6) to preserve the continuity of the solution, even though the source function is non-smooth. Now, we produce the effect of these conditions and their discretizations for the coupled system (1) with the boundary conditions (5.4a), continuity condition (5.5), and Kirchhoff's condition (5.6). This is shown in Table 7, where we also observe an almost second-order convergence overall perturbation parameters from the range $\varepsilon \in \mathbf{S}$.

Next, we discuss the conditions at the free vertices of the $3$-star graph. The maximum error and corresponding order of convergence on each edge of the $3$-star graph are provided for this example, with the boundary conditions (5.4b), continuity condition (5.5), and Kirchhoff's condition (5.6). The results shown in Tables 8–10, clearly indicate almost second-order uniform convergence. The coupled system in this example with boundary conditions (5.4b) with (5.5) and (5.6), also demonstrates the almost second-order convergence, as presented in Table 11, for $\varepsilon \in \mathbf{S}$. Similarly, The results shown in Tables 12–14, clearly indicate almost second-order uniform convergence for each solution component $u_1,$ $u_2,$ and $u_3$ respectively for (5.4c) along with the (5.5) and (5.6). Finally, considering the mixed boundary conditions (5.4c), the numerical approximation for the coupled system of the $3$-star graph (see Figure 1) shows the almost uniform second-order convergence, as seen in Table 15, for $\varepsilon \in \mathbf{S}$. Note that the obtained accuracy is higher than the discretizations, where the stiffness matrices are generated through the simple upwind two-point schemes for first-order derivatives. We want to note that the round-off errors can affect the order of accuracy when the tolerance is minimal, like $10^{-6}$, and dominate the errors.

The boundary and interior layers, along with their sharpness, are visible across all edges of the solution components, as shown in Figures 3a–5a for very small values of perturbation parameters. Additionally, Figures 3b–5b illustrate that the experimental rate of convergences match with the established theoretical rate of convergences over all edges of the 3-star graphs. Moreover, Figure 6 demonstrates that the solution components along the corresponding edges agree with the continuity condition (4.10).

Hence, we conclude that a suitable transformation of row entries of stiffness matrices (corresponding to the three-point scheme at interior and boundary points and the five-point scheme at the point of discontinuity) can lead to higher-order accurate solutions (compared to first order accurate solutions as available in ^[49]) for nonlinear reaction-dominated problems that are connected by the $3$-star graph.

6.   Conclusions

We examine a semi-linear singularly perturbed reaction-diffusion problem with non-smooth data on a $k$-star graph. Here, we consider all possible boundary conditions at the free boundary, specifically at the tail of the edge. At the junction point of the graph, we impose the natural and relevant requirements of continuity and Kirchhoff's junction law. Kirchhoff's junction law involves the first derivative at the junction point. Hence this makes the original problem coupled. First, we apply the flux condition at the free boundary, which is approximated by a first-order upwind approximation. To enhance this accuracy, we employ a three-point scheme for first-order derivatives and a central difference scheme for second-order derivatives. Generally, these approximations do not lead to an M-matrix or diagonal dominant structure of the stiffness matrix. Hence, we transform this scheme and establish sufficient conditions under which the stiffness matrix will be diagonally dominant. Additionally, we conduct a rigorous convergence analysis of this scheme, demonstrating almost second-order accuracy for each problem defined on the edges of the graph and for the coupled system defined on the star graph as well. Our numerical experiments further confirm the practical significance of our scheme, supported by experimental evidence. In the future, we aim to extend this work to a system of weakly and strongly coupled differential equations connected by more general graphs, such as cycles and complete graphs. In these cases, the computational costs will be high and will be considered future goals.

Author contributions

All authors have contributed equally.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

Higinio Ramos and Pratibhamoy Das are guest editors for Networks and Heterogeneous Media and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.

Acknowledgments

The author, Pratibhamoy Das, acknowledges the support of the Science and Engineering Research Board, Government of India, under Project No. MTR/2021/000797, which facilitated the completion of this work. This work was also supported by the University Grants Commission (UGC), India, through the Junior Research Fellowship (JRF) awarded to Dilip Sarkar, who gratefully acknowledges this support. The author, Shridhar Kumar, gratefully acknowledges the Indian Institute of Technology Patna for providing the institute fellowship.