Non-separated inclusion problem via generalized Hilfer and Caputo operators

1.   Introduction

Non-linear partial differential equations are widely used by scientists and engineers to interpret various physical phenomena, such as shallow water wave behavior, viscoelastic behavior of fluids, and shock wave behavior in fluids. In different branches of science and engineering, various physical processes are described theoretically through ordinary or partial differential equations. The Lane-Emden equation ^[1,2], Emden Fowler equation ^[3,4,5], and Riccati equation ^[6] are examples of ordinary differential equations which are used to describe various physical, chemical, and biological phenomena, including the problems in fluid flow and engineering. Burger's equation is used to study the physical processes in aerospace engineering and fluid dynamics ^[7]. The non-local diffusion equation is used to describe neutron transport in a nuclear reactor ^[8]. The Burger-Huxley equation, Burger-Fisher equation ^[7,9,10], Hodgkin-Huxley equation ^[11,12], and Fitzough-Nagumo equation ^[13,14] have numerous applications in the field of physics, biology, fluid dynamics, engineering, optics, and plasma physics, among others.

Consider the one-dimensional non-linear reaction-diffusion equation:

where $f(u)$ is the non-linear source function, $\varepsilon(\tau)$ is the function of time whose value lies between 0 and 1, and $I = (a, b)$. The initial condition is taken to be continuous whereas boundary conditions are assumed to be of the Dirichlet type. The initial condition is defined as:

Both the homogeneous Dirichlet type and non-homogeneous Dirichlet type boundary conditions have been considered and are defined as:

where, $K_{1}$ and $K_{2}$ are continuous functions of $\tau$.

The simplified form of Eq (1.1) for $\varepsilon(\tau) = 1$ was first proposed by Hodgkin and Huxley to explain the propagation of action potentials in the squid giant axiom. The generalized Burger's, generalized Burger's-Huxley, and generalized Burger's-Fisher equations are the special cases of Eq (1.1) with non-linear advection term and source term. These equations have numerous applications in reaction mechanics and diffusion transport in fluids and semi-solid particles.

A variety of numerical techniques have been developed by different investigators to study the solution behavior of non-linear partial differential equations ^{[1,5,6,15,16,17,18]}. Lagurre polynomials were utilized in ^[1] to study the behavior of the Lane-Emden equation whereas ^[2] has used the Bessel collocation approach. The Emden-Fowler equation was solved using artificial neural networks by ^[5]. Orthogonal collocation on finite elements was used by ^[19] to discretize the non-linear heat conduction equation. The coupled Burger and coupled Burger-Huxley equations were discretized using the Bessel collocation method in ^[7]. The two phase non-linear reaction diffusion model was discretized using cubic Hermite splines in ^[20], whereas singularly perturbed reaction diffusion equations were solved by orthogonal cubic splines of order 3 and order 5 in ^[21]. The non-local Nagumo type equation was solved by using a finite volume scheme in ^[22]. The non-local fourth-order sub-diffusion equation was solved using orthogonal spline collocation in ^[18], and the non-local heat model was discretized using the sinc collocation method in ^[23].

In the present study, cubic Hermite splines with the Crank-Nicolson scheme is used to discretize the non-linear reaction-diffusion equation. The cubic Hermite splines are third-order Hermite splines that interpolate the function as well as its first-order derivative at node points. These splines have an advantage over the cubic B-splines as these splines do not involve any fictitious node points and thus are easily adapted to a computer program. The combination of cubic Hermite splines with the CN scheme in the time direction gives an advantage to discretize the equation in both the the space and, the time directions. It reduces the stiffness not only in space direction but also in time direction.

The present paper has been weaved into different sections and sub-sections, from the introduction to the conclusions. The CN scheme is discussed in Section 2 and the cubic Hermite splines are explained in Section 3. The stability and convergence of the proposed technique is discussed in Section 4, whereas Section 5 is composed of the numerical applications. The conclusions of the study are presented in Section 6.

2.   Crank-Nicolson scheme

A variety of numerical techniques for solving two point boundary value problems in the time direction have been developed. The Crank-Nicolson (CN) and Alternate Direction Implicit scheme are two such examples. The Alternate Direction Implicit scheme, commonly known as ADI scheme, was introduced by Peaceman and Rachford to study second-order parabolic boundary value problems over the rectangular domain. ADI methods are capable of solving multidimensional problems into collection of one-dimensional problems ^[24]. Other than the simple ADI scheme or the Peaceman-Rachford ADI scheme, ^[25] discussed certain other classifications of ADI schemes, such as the Douglas scheme, Craig-Sneyd scheme, modified Craig-Sneyd scheme, Hundsdorfer-Verwer scheme, etc. Crank-Nicolson scheme, or CN scheme, is also an efficient time integration scheme to solve one-dimensional reaction-diffusion problems. In the present study, the CN scheme as proposed by ^[26] operates on the mean of a function over the interval $[\tau_{j}\; \; \tau_{j+1}]$ with a uniform distribution of points ^[16,27].

Define the partition $\pi_{\tau}:0 = \tau_{0} < \tau_{1} < \ldots < \tau_{M} = T$ with $\Delta t = \tau_{j+1}-\tau_{j}$. Consider the general non-linear reaction-diffusion equation defined by Eq (1.1)

For convenience, we write $u(\zeta, \tau_{j+1}) = u^{j+1}$.

Lemma 1. ^[16,28] If $\left|\frac{\partial^{i}u}{\partial \tau^{i}}\right|\leq C, \ i = 0, 1$, and $(\zeta, \tau)\in\bar\Omega$, then the local truncation error $e$ is given by

and the global truncation error $E$ is given by

Theorem 1. ^[16,28] Let U($\zeta$, $\tau$) be a function such that ${\rm{£}}$U = 0 and $U\leq C$. Then, $U_\tau$ is also bounded by some constant C, where C is an arbitrary constant.

For the implementation of the technique on non-linear singularly perturbed equations, let $u^{j+1}$ be the approximating polynomial function at time step $\tau_{j+1}$. Also, let $\Delta t = T/M, \ M > 0$. Then, $\tau_j = j\Delta t$ for $j = 0, 1, 2, …, M$. The semi-discretized form of Eq (2.1) follows as:

with corresponding boundary conditions are given by Eqs (1.2) and (1.3), respectively. For any function F,

Throughout this paper, C denotes the generic constant.

After the implementation of the CN scheme, Eq (1.1) reduces to the semi-discretized form as follows:

The non-linear function term $f(u^{j+1})$ is quasi-linearized using the formula given by ^[29]:

where, $\frac{df}{du_j}$ represents the derivative of $f(u)$ at $\tau_j$. After substituting Eq (2.5) in Eq (2.4) and rearranging the terms, one gets the iterated scheme

where $f_u^j$ represents $\frac{df}{du_j}$ for convenience. The iterative form of Eq (2.6) can be expressed in operator form at the $(j+1)^{th}$ time step in the following manner:

where, $\mathbb{L}(u^{j+1}) = u^{j+1}-\frac{\varepsilon_j \Delta t}{2} u_{\zeta\zeta}^{j+1}-\frac{\alpha_j\Delta t f_{u^j}}{2} u^{j+1}$ and $G_{\tau_j} = u^{j}+\frac{\varepsilon_j \Delta t}{2} u_{\zeta\zeta}^{j}-\frac{\alpha_j\Delta t f_{u^j}}{2} u^{j}+\alpha_j \Delta tf(u^j)$. Since there is no singular term on the right-hand side, Eq (2.7) is convergent, which leads to the solution of Eq (1.1).

3.   Cubic Hermite splines

Orthogonal splines represent the piecewise orthogonal polynomials used to interpolate the function at node points. Hermite and Lagrange's interpolating polynomials are such orthogonal splines that are often used to interpolate functions. Hermite interpolating polynomials of order '$k$' are considered as an extension of $k^{th}$-order Lagrangian interpolating polynomials ^{[7,9,10,16,20,21,27]}. A continuity condition is imposed at node points in the Lagrangian interpolating polynomials as these polynomials can only interpolate the function at node points. Hermite interpolating polynomials of order '$2k+1$' interpolate the function as well as its $k^{th}$ order derivative at node points. This feature of Hermite interpolating polynomials makes it superior over Lagrangian interpolating polynomials. Hence, th $k^{th}$-order Hermite polynomial in $\zeta$ is a polynomial of order $2k+1$, and, therefore, cubic Hermite interpolating polynomials are a particular case of general Hermite interpolating polynomials for $k = 1$. This consists of two node points and two tangents at these node points. The explanation of cubic Hermite splines is given hereunder.

Consider an interval $I = (a, b)$, and let $\pi$ be the partition of $(a, b)$ such that:

Let $P_{3}$ be the set of all polynomials of degree less than or equal to 3 defined on $[\zeta_{i-1}, \zeta_{i}]$. Let $v$ be a continuously differentiable function defined on $\bar{I}$ such that $v$ is the cubic Hermite approximation of $u$ on $I$. Let $\mathbb{M}_{\mu}$ be the space of all continuously differentiable functions defined on $[\zeta_{i-1}, \zeta_{i}]$ such that

If boundary conditions are of the homogeneous type, then the dimension of $\mathbb{M}_{\mu} = 2n$, and if the boundary conditions are of the non-homogeneous type, then the dimension of $\mathbb{M}_{\mu} = 2n+2$ ^[24].

The cubic Hermite polynomials $\mathcal{P}_{\mu}(\zeta)$ and $\mathcal{Q}_{\mu}(\zeta)$ for $k = 1$ are defined as

These piecewise cubics are designed such that the following identities hold: $\mathcal{P}_{\mu}(\zeta_{i}) = \delta_{\mu i}$; $\mathcal{P}'_{\mu}(\zeta_{i}) = 0$; $\mathcal{Q}_{\mu}(\zeta_{i}) = 0$; $\mathcal{Q}'_{\mu}(\zeta_{i}) = \delta_{\mu i}$.

To apply collocation in the interval $[\zeta_{i-1}, \zeta_{i}]$, a new variable $\xi$ is introduced in such a way that, as $\zeta$ varies from $\zeta_{i-1}$ to $\zeta_{i}$, the variable $\xi$ varies from 0 to 1. It converts the polynomials defined in Eqs (3.2) and (3.3) into the following form:

These four cubic Hermite interpolating polynomials form the basis of $\mathbb{M}_{\mu}$ for $k = 1$. These interpolating polynomials are bounded and vary from $-0.5$ to $1$. The graphical representation of these polynomials is shown in Figure 1. In Table 1, the values of these polynomials and their respective derivatives are shown at boundary points 0 and 1.

3.1.   Collocation points

Let $\mathbb{G}$ be the space of all those points which are zeros of the shifted Legendre polynomials of order $n$. Let $\eta_{j}$ be the node points of the Gauss-Legendre quadrature formula defined on $[0, 1]$, where $j$ = $1, 2$, $\ldots$, $n-1$. The collocation points are defined by

where $h_{\mu} = \zeta_{\mu}-\zeta_{\mu-1}$. Runge's divergence formula also states that non-uniform collocation points give less error compared to uniform collocation points ^[15,19].

Theorem 2. ^[30] The necessary and sufficient condition for a simple set $\mathbb{B}_n(\xi)$ of real polynomials to be orthogonal w.r.t. the weight function $\omega(\xi)$ on $a \leq \xi \leq b$ is

where, $P_n(\xi) \in \mathbb{B}_n(\xi)$. The roots of these polynomials are always real and distinct and lie in the interval $a \leq \xi \leq b$.

The collocation points are chosen to be the zeros of the Legendre polynomials as these polynomials have the tendency to reduce the error at corners as well as on averages ^[19]. The zeros of the Legendre polynomials are calculated using the following recurrence relation given in ^[31]:

The details of these polynomials are given elsewhere ^[19,31].

An application of orthogonal collocation using cubic Hermite splines as basis polynomials is applied with consideration of the approximate solution of a semi-discrete scheme within each sub-element of the partition $\pi_{\zeta}$. A linear combination of four stencils with time dependent coefficients is proposed as follows:

where $\sigma_{i}^{j}$ are the unknown coefficients to be determined. After substituting Eq (3.4) into Eq (2.6), the following system of equations is obtained:

At the $k^{th}$ collocation point, Eq (3.5) reduces to

where $H_{ki}$ and $H''_{ki}$ are the Hermite spline and the second-order derivative of the Hermite spline at the $k^{th}$ collocation point, respectively, and $f_k^j = f\bigg(\sum_{i = 1}^4 H_{ki}\sigma^{j}\bigg)$. This system of equations reduces to the matrix form

where $\mathcal{Q}$ and $\mathbb{H}$ are almost quad-diagonal dominant matrices of order $2n\times2n$, where $n$ denotes the number of sub-elements in $\pi_{\zeta}$. $\rho^{n} = [\sigma_{1}^{n}, \sigma_{2}^{n}, \ldots, \sigma_{2n}^{n}]'$; $F^{n} = [f^{n}(\xi_{1}), f^{n}(\xi_{2}), \ldots, f^{n}(\xi_{2n})]'$. The resulting system of matrix equations can be solved by any iterative technique.

In combined form, Eq (3.7) can be written as

where $\mathcal{H} = \mathbb{H}\rho^{j}+F^{j}$.

4.   Stability and convergence analysis

The implementation of any numerical technique depends upon its stability. A method is said to be stable if its $j^{th}$ iteration is bounded ^[9,10]. The function $f(u)$ is bounded over the domain $I \times [0, T]$ and satisfies the Lipschitz condition. $\frac{df}{du}$ is bounded over the desired domain such that $\mid f(u) \mid \leq K_1$ and $\mid \frac{df}{du} \mid \leq K$.

Define the inner product $< u, v > = \int_{\Omega}{u(\zeta)v(\zeta) d\zeta}$, where $< u, u > = \mid\mid u\mid\mid^2$.

Lemma 2. Let $f(u)$ be defined on the domain $I \times [0, T]$ such that $f(u)$ satisfies Lagrange's theorem. If the derivative $\frac{df}{du}$ is bounded over the domain $I \times [0, T]$ s.t. $\mid \frac{df}{du} \mid \leq K$, then $f(u)$ satisfies the Lipschitz condition

Theorem 3. Cauchy-Schwarz inequality: Let $u$ and $v$ be any two functions defined on $\mathcal{R}$. Then,

Lemma 3. ^[28] Let $u(\zeta, 0)$ be the positive continuous initial solution of the non-linear reaction-diffusion equation with $0 \leq u(\zeta, 0)\leq 1$. Then, the bounds for solution $u(\zeta, t)$ remains constant with time.

Throughout this paper, $C$ denotes the generic constant.

Lemma 4. ^[28] Let $u(\zeta, \tau)$ be the solution of Eq (1.1) in $\bar{\Omega}$. Then, the bounds on $u(\zeta, \tau)$ is given by

Lemma 5. ^[28] Maximum Principle: Let $u(\zeta, \tau)$ be the solution of Eq (1.1) in $\bar{\Omega}$, $u(\zeta, \tau) \geq 0$ on $\partial\Omega$, and $\mathcal{L}u(\zeta, \tau) \geq 0$ on $\Omega$. Then, $u(\zeta, \tau) \geq 0$ on $\bar{\Omega}$.

Lemma 6. ^[28] Let $u(\zeta, \tau)$ be the solution of Eq (1.1) in $\bar{\Omega}$. Then, the bound on the derivative of u w.r.t. $\tau$ is given by

Lemma 7. ^[28] Let $u(\zeta, \tau)$ be the solution of Eq (1.1) in $\bar{\Omega}$. Then, the bound on the derivative of u w.r.t. $\zeta$ is given by

Lemma 8. ^[28] Let $u(\zeta, \tau)$ be the solution of Eq (1.1) in $\bar{\Omega}$. Then, the bound on the derivatives of u are given by

Theorem 4. Let $u^j$ be the approximate solution of Eq (1.1) at $\tau_j$. Then, the following inequality holds:

Proof. Now, using the bound for $\frac{df}{du}$ and redefining Eq (2.6), the following equation is obtained:

where $\eta_1 = \Big(1-\frac{\alpha_j\Delta t K}{2}\Big)$ and $\eta_2 = \frac{\varepsilon_j \Delta t}{2}$.

Taking the inner product of both sides of $u^{j+1}$

Now, using the fact that $- < u_{\zeta\zeta}^{j+1}, u^{j+1} > \ \geq \ 0$, and from Cauchy-Schwarz inequality that $\mid < u^{j+1}, u^j > \mid \ \leq \ \mid\mid u^{j+1} \mid\mid$ $\mid\mid u^j \mid\mid$, Eq (4.8) takes the following form:

Dividing both sides by $\eta_1\mid\mid u^{j+1}\mid\mid$,

Using Lemmas 2, 3, and 8 and the principle of induction, it is clear that right-hand side is bounded. This implies $\mid\mid u^{j+1}\mid\mid$ is bounded, and hence the result. □

Lemma 9. ^[32] If any Hermite interpolation polynomial based on at most two distinct points has a unique solution in $\Omega$, then any Hermite interpolation polynomial has a unique solution in $\Omega$. This implies that any Hermite interpolation polynomial has a unique solution in the polynomial space $\mathbb{P}_n$ of degree $n$ on I.

Lemma 10. ^[32] The univariate Hermite interpolating polynomials are regular for any set of node points. Moreover, Hermite interpolating polynomials are regular for any nodal set of points as well as for any choice of derivatives to be interpolated.

Lemma 11. Let $\bar{\mathbb{H}}$ be the space of all Hermite interpolating polynomials of order 3 defined on [0 1]. Then, the Hermite splines of order 3 are bounded with upper bound unity. Moreover, $\sum_{i = 1}^{4}\mid H_i(\xi)\mid \leq 4$ for all $0 \leq \xi \leq 1$.

Lemma 12. ^[27,33,34] Let H $\in H_{\Delta \xi}^{4}(\xi)$ be the piecewise Hermite spline of degree $3$ over the subinterval $[\xi_i, \xi_{i+1}]$ approximating $\boldsymbol{U}$ $\in C^{4}$[a, b]. Then, $\mid\mid \boldsymbol{H}^{(r)}-\boldsymbol{U}^{(r)}\mid\mid_{\infty} \leq C\gamma_r h^{4-r}, \; r = 0, 1, 2, 3$, where $H^{(r)}$ represents the $r^{th}$-order derivative of H. The values of $\gamma_r$ are given in ^[27,33,34].

Theorem 5. Let $u(\zeta)$ be the exact solution of Eq (1.1) and $\bar{u}(\zeta)$ be the approximate solution of Eq (1.1) in the space $\mathbb{H}^3$ of cubic Hermite interpolating polynomials of order 3 such that $\bar{u}(\zeta) \in C^4[a, b]$. Then, the uniform error estimate is given by

Proof. Let ${{\bf{H}}}(\zeta)$ be the unique cubic Hermite spline interpolate of $u(\zeta)$, and $\pi_{\zeta} = {a = \zeta_{0}\leq \zeta_{1}\leq, ..., \leq \zeta_{n} = b}$ be the equi-spaced partition of [a, b] with uniform step-size $h$.

Consider

where $\bar{\alpha}/2$ = $\big| 1-\frac{\alpha_j f_{u^j}\Delta t}{2} \big|$.

Using Lemma 12, the following inequality is obtained:

where $C$ is the generic constant and the values of $\gamma_0$ and $\gamma_2$ are given in ^[27,33,34].

Now, consider

Now, we use the fact that $\mathbb{L}(\bar{u}(\zeta_{k})) = \bar{G}_{\tau_j}(\zeta_{k})$ and $\mathcal{Q}\rho = \mathcal{H}$.

This implies that $\mathcal{Q}(\rho - \bar{\rho}) = \mathcal{H}-\bar{\mathcal{H}}$, and using Eq (4.14) we get

Also, $\mathcal{Q}$, the matrix of collocation coefficients, is bounded. Therefore, $||\mathcal{Q}||\leq C$ and

Now, $\bar{u}(\zeta)-H(\zeta) = \sum_{i = 1}^{4}H_{i}(\zeta)(\sigma_{i}-\bar{\sigma}_{i})$, and we get

According to ^[33],

Now, using the triangle inequality,

□

Theorem 6. Let $u(\zeta, \tau)$ represent the exact solution and $\bar{u}(\zeta, \tau)$ be the approximate solution. Then, the uniform error estimate is given by

Proof. The convergence of the proposed method in the time direction is $O(\Delta t^2)$ as given in Lemma 1, and from Theorem 5 the error estimate in the space direction is $Ch^2 \big(\gamma_2 \varepsilon_j\Delta t + \gamma_0(1+\bar{\alpha})h^2\big)$. Therefore, combining the results from Lemma 1 and Theorem 5, the following error estimate is obtained:

□

5.   Results and discussions

The error analysis of any technique is incomplete if it does not consider the Eucledian and supremum norms, also known as the $L_2$-norm and $L_{\infty}$-norm, respectively. These norms are defined as

Problem 1. Consider the Fitzough-Nagumo equation

The boundary as well as initial conditions are given below:

Problem 1 was solved using cubic Hermite splines with $a = -10$ and $b = 10$. The numerical findings of the $L_{2}$-norm and the $L_{\infty}$-norm have been calculated for $\mu = 0.25$ and $\mu = -0.5$, which is presented in Table 2. The norms show the stability and the efficiency of the technique. In Tables 3 and 4, the absolute error is shown for $\mu = 0.25$ and $\mu = -0.5$ for varying number of node points. It is clear from these tables that the absolute error varies from order $10^{-6}$ to $10^{-4}$. In Table 5, the numerical values obtained from present technique are compared with the values given in ^[35]. The $L_{\infty}$-norm is found to be better than the values given in the literature, and the $L_2$-norm is found to be on par with the $L_2$-norm given in the literature, and at some points even better. The graphical representation of the numerical solutions is visualized in Figures 2–5 for varying values of $\mu$. It is observed from these figures that the numerical values are bounded and lie between 0 and 1. From these tables and figures, it is observed that these quantitative results show the efficiency and accuracy of the method.

Problem 2. Consider the variable order Fitzough-Nagumo equation

Initially,

Problem 2 was solved using cubic Hermite splines as basis functions. To check the stability and the efficiency of the technique, the numerical findings of the $L_{2}$-norm and $L_{\infty}$-norm were calculated for $\mu = 0.25$ and $\mu = 0.05$. The numerical values of the $L_2$-norm and the $L_{\infty}$-norm are presented in Table 6. It is observed that the $L_2$-norm varies from $10^{-6}$ to $10^{-5}$, whereas the $L_{\infty}$-norm varies from $10^{-5}$ to $10^{-3}$. In Tables 7 and 8, the absolute error is shown for $\mu = 0.25$ and $\mu = 0.05$ for varying number of node points. It is observed that the absolute error stabilizes after 100 node points for $\mu = 0.25$ and $\mu = 0.05$. In Table 9, the $L_2$-norm and $L_{\infty}$-norm are compared to the values given in ^[35]. Both the $L_2$-norm and $L_{\infty}$-norm are found to be better than the values given in the literature. It is also observed that the absolute error is reduced for small values of $\mu$. The graphical representation of the numerical solution is presented in Figures 6–9 for different values of $\mu$. It is clear from these figures that $u$ is bounded and lies between 0 and 0.25.

Problem 3. Consider the Newell-Whitehead-Segel equation

The boundary conditions can be taken as

The initial condition is

The above problem is solved using cubic Hermite splines. The absolute error obtained from the cubic Hermite splines is compared with values given in ^[36] at different time periods in Tables 10–14. It is observed that the values obtained from the cubic Hermite splines are better than those given in the literature. In Table $15$, the absolute error is shown for $\varepsilon = 1$ and $\tau = 0.5$ for different node points. It is observed that the absolute error stabilizes after 200 node points. The graphical representation of the numerical solution is shown in Figures 10 and 11 for different values of $\tau$. The 2D and 3D graphs show that the absolute value of $u$ is less than 1.

6.   Conclusions

Orthogonal collocation with cubic Hermite splines was applied on non-linear reaction-diffusion type equations. Two equations of Fitzough-Nagumo type and one of Newell-Whitehead-Segel type were solved numerically using cubic Hermite splines with the CN scheme. It is observed that the proposed algorithm is applicable to non-linear reaction-diffusion equations with both variable and constant coefficients. The absolute error is found to be of order $10^{-4}$, which proves the efficiency of the proposed algorithm. The proposed technique can also be applied on higher-dimensional linear as well as non-linear problems. The proposed technique can also be applied to fractional-order partial differential equations.

Use of AI tools declaration

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

Acknowledgments

Dr. Shelly Arora is thankful to DST-SERB for providing grant via SPG/2022/001269 and Mr. Abdul-Majeed Ayebire is thankful to ICCR for providing financial assistance via reference number LY7930408566174.

Authors are thankful to the worthy reviewers for their valuable comments to improve this manuscript.

Conflict of interest

Authors declare that they have no conflict of interest.