Soliton solutions and stability analysis of the stochastic nonlinear reaction-diffusion equation with multiplicative white noise in soliton dynamics and optical physics

Nafissa T. Trouba; Huiying Xu; Mohamed E. M. Alngar; Reham M. A. Shohib; Haitham A. Mahmoud; Xinzhong Zhu; Nafissa T. Trouba; Huiying Xu; Mohamed E. M. Alngar; Reham M. A. Shohib; Haitham A. Mahmoud; Xinzhong Zhu

doi:10.3934/math.2025086

AIMS Mathematics

2025, Volume 10, Issue 1: 1859-1881. doi: 10.3934/math.2025086

Previous Article Next Article

Research article

Soliton solutions and stability analysis of the stochastic nonlinear reaction-diffusion equation with multiplicative white noise in soliton dynamics and optical physics

1.
School of Computer Science and Technology, Zhejiang Normal University, Jinhua 321004, China
2.
Zhejiang Institute of Photo-electronics, Zhejiang, Jinhua 321004, China
3.
Department of Mathematics Education, Faculty of Education & Arts, Sohar University, Sohar, Oman
4.
Basic Science Department, Higher Institute of Management Sciences & Foreign Trade, Cairo, 379, Egypt
5.
Industrial Engineering Department, College of Engineering, King Saud University, Riyadh 11421, Saudi Arabia

Received: 20 November 2024 Revised: 16 January 2025 Accepted: 22 January 2025 Published: 24 January 2025
MSC : 34L30, 35A24, 35B35, 35C08, 74J35

In this article, we explored the stochastic nonlinear reaction-diffusion (RD) equation under the influence of multiplicative white noise. To obtain novel soliton solutions, we employed two powerful analytical techniques: the unified Riccati equation expansion method and the modified Kudryashov method. These methods yield a diverse set of soliton solutions, including combo-dark solitons, dark solitons, singular solitons, combo-bright-singular solitons, and periodic wave solutions. We also performed a comprehensive stability analysis of the stochastic nonlinear RD equation with multiplicative white noise. The findings provide valuable insights into the behavior of solitons in stochastic nonlinear systems, with significant implications for fields such as mathematical physics, nonlinear science, and applied mathematics. These results hold particular relevance for soliton dynamics in optical physics, where they can be applied to improve understanding of wave propagation in noisy environments, signal transmission, and the design of robust optical communication systems.

Keywords:

Citation: Nafissa T. Trouba, Huiying Xu, Mohamed E. M. Alngar, Reham M. A. Shohib, Haitham A. Mahmoud, Xinzhong Zhu. Soliton solutions and stability analysis of the stochastic nonlinear reaction-diffusion equation with multiplicative white noise in soliton dynamics and optical physics[J]. AIMS Mathematics, 2025, 10(1): 1859-1881. doi: 10.3934/math.2025086

Related Papers:

[1]	Ali H. Tedjani, Sharifah E. Alhazmi, Samer S. Ezz-Eldien . An operational approach for one- and two-dimension high-order multi-pantograph Volterra integro-differential equation. AIMS Mathematics, 2025, 10(4): 9274-9294. doi: 10.3934/math.2025426
[2]	Ahmed M. Rajab, Saeed Pishbin, Javad Shokri . Analyzing the structure of solutions for weakly singular integro-differential equations with partial derivatives. AIMS Mathematics, 2024, 9(9): 23182-23196. doi: 10.3934/math.20241127
[3]	Obaid Algahtani, M. A. Abdelkawy, António M. Lopes . A pseudo-spectral scheme for variable order fractional stochastic Volterra integro-differential equations. AIMS Mathematics, 2022, 7(8): 15453-15470. doi: 10.3934/math.2022846
[4]	A. H. Tedjani, A. Z. Amin, Abdel-Haleem Abdel-Aty, M. A. Abdelkawy, Mona Mahmoud . Legendre spectral collocation method for solving nonlinear fractional Fredholm integro-differential equations with convergence analysis. AIMS Mathematics, 2024, 9(4): 7973-8000. doi: 10.3934/math.2024388
[5]	Chuanli Wang, Biyun Chen . An $hp$ -version spectral collocation method for fractional Volterra integro-differential equations with weakly singular kernels. AIMS Mathematics, 2023, 8(8): 19816-19841. doi: 10.3934/math.20231010
[6]	Mohammed A. Almalahi, Satish K. Panchal, Fahd Jarad, Mohammed S. Abdo, Kamal Shah, Thabet Abdeljawad . Qualitative analysis of a fuzzy Volterra-Fredholm integrodifferential equation with an Atangana-Baleanu fractional derivative. AIMS Mathematics, 2022, 7(9): 15994-16016. doi: 10.3934/math.2022876
[7]	Chuanhua Wu, Ziqiang Wang . The spectral collocation method for solving a fractional integro-differential equation. AIMS Mathematics, 2022, 7(6): 9577-9587. doi: 10.3934/math.2022532
[8]	Khadijeh Sadri, David Amilo, Kamyar Hosseini, Evren Hinçal, Aly R. Seadawy . A tau-Gegenbauer spectral approach for systems of fractional integro-differential equations with the error analysis. AIMS Mathematics, 2024, 9(2): 3850-3880. doi: 10.3934/math.2024190
[9]	Mian Bahadur Zada, Muhammad Sarwar, Reny George, Zoran D. Mitrović . Darbo-Type $\mathcal{Z}_{\rm{m}}$ and $\mathcal{L}_{\rm{m}}$ contractions and its applications to Caputo fractional integro-differential equations. AIMS Mathematics, 2021, 6(6): 6340-6355. doi: 10.3934/math.2021372
[10]	Hawsar Ali Hama Rashid, Mudhafar Fattah Hama . Approximate solutions for a class of nonlinear Volterra-Fredholm integro-differential equations under Dirichlet boundary conditions. AIMS Mathematics, 2023, 8(1): 463-483. doi: 10.3934/math.2023022

Abstract

1. Introduction

The goal of this paper is to develop Legendre spectral tau approaches for the following high-order Pantograph Volterra-Fredholm integro-differential equations (P-VF-IDEs):

● One-dimensional case:

$\begin{equation} \begin{split} \sum\limits_{k = 0}^{n}b_k y^{(k)}(x) = g(x)+\sum\limits_{k = 0}^{n}c_k y^{(k)}(p_k x)+\sum\limits_{k = 0}^{n}\int_{0}^{p_k x}d_k y^{(k)}(s)ds+\int_{0}^{x}e y(s)ds+\int_{0}^{\sigma}f y(s)ds, \end{split} \end{equation}$

(1.1)

subject to

$\begin{split} \frac{d^{k} y(0) }{d x^{k}} = \varrho_k, \qquad \quad k = 0, 1, \cdots, n-1. \end{split}$

● Two-dimensional case:

$\begin{equation} \begin{split} &\sum\limits_{k_1 = 0}^{n_1}\sum\limits_{k_2 = 0}^{n_2}b_{k_1, k_2}\frac{\partial^{k_1+k_2} y(x_1, x_2) }{\partial x_1^{k_1}\partial x_2^{k_2}} = g(x_1, x_2)+\int_{0}^{\sigma_2}\int_{0}^{\sigma_1}c y(s_1, s_2) ds_1ds_2 + \int_{0}^{x_2}\int_{0}^{x_1}d y(s_1, s_2) ds_1ds_2 \\ & + \sum\limits_{k_1 = 0}^{n_1}\sum\limits_{k_2 = 0}^{n_2}\int_{0}^{p_{1} x_1}\int_{0}^{p_{2} x_2} e_{k_1, k_2} \frac{\partial^{k_1+k_2} y(s_1, s_2)}{\partial x_1^{k_1}\partial x_2^{k_2}} ds_1ds_2+ \sum\limits_{k_1 = 0}^{n_1}\sum\limits_{k_2 = 0}^{n_2}f_{k_1, k_2}\frac{\partial^{k_1+k_2} y(p_1x_1, p_2x_2)}{\partial x_1^{k_1}\partial x_2^{k_2}}, \end{split} \end{equation}$

(1.2)

subject to

$\begin{split} \frac{\partial^{k_1} y(0, x_2) }{\partial x_1^{k_1}} = \varrho_{1, k_1}(x_2), \qquad \quad k_1 = 0, 1, \cdots, n_1-1, \quad 0\leq x_2 \leq \sigma_2, \\ \frac{\partial^{k_2} y(x_1, 0) }{\partial x_2^{k_2}} = \varrho_{2, k_2}(x_1), \qquad \quad k_2 = 0, 1, \cdots, n_2-1, \quad 0\leq x_1 \leq \sigma_1.\\ \end{split}$

The problem under consideration involves several key parameters: the coefficients $b_k, c_k, d_k, e, f$ (1D case) and $b_{k_1, k_2}, e_{k_1, k_2}, f_{k_1, k_2}, c, d$ (2D case), which weight the contributions of derivatives, delay terms, and integral operators; the pantograph-type proportional delays $p_k \in (0, 1)$ (1D) and $p_1, p_2 \in (0, 1)$ (2D), introducing proportional delays such as $y(p_k x)$ ; the fixed integration limits $\sigma, \sigma_1, \sigma_2 > 0$ , defining the domains of Fredholm integrals; and the initial conditions $\varrho_k$ (1D) and $\varrho_{1, k_1}(x_2), \varrho_{2, k_2}(x_1)$ (2D), which prescribe initial values for $y$ and its derivatives. Special cases demonstrating the problem's generality include lower-order systems where $n = 0$ reduces the 1D equation to a delay Volterra-Fredholm integral equation; non-delayed systems with $p_k = 1$ , recovering classical integro-differential equations; pure Volterra dynamics achieved by omitting $f$ (1D) or the Fredholm term (2D) to isolate Volterra behavior, or removing $e$ (1D) or the Volterra term (2D) to focus on Fredholm interactions; single-integral equations formed by discarding mixed integrals; and pantograph-free cases with $c_k = 0$ and $f_{k_1, k_2} = 0$ . These examples illustrate the framework's flexibility in unifying several problems while addressing novel combinations of high-order derivatives, multi-delays, and mixed integrals.

High-order P-VF-IDEs represent a class of mathematical models that integrate effects of time delays proportional to the independent variable, known as pantograph terms, with both Volterra and Fredholm-type integro-differential equations. These equations arise in various fields, including control theory, population dynamics, and systems biology, where they model processes influenced by historical states and time-dependent interactions ^[1,2]. A pantograph is a mechanical device that collects electric current from overhead wires to supply power to electric trains and trams. Its name originates from a historical mechanical device used for copying and scaling drawings based on similar principles of motion transfer and scaling. In electric transportation systems, the pantograph ensures a consistent connection with overhead conductors, facilitating a continuous and stable flow of electricity to the vehicle. The term also originates from early drafting instruments used for geometric replication and proportional scaling ^[3,4]. Numerical approximations of P-VF-IDEs face significant challenges due to the integration of proportional delays and mixed integral terms. These complexities necessitate the use of adaptive numerical methods to achieve high accuracy. Common numerical techniques include the collocation approach ^[5,6], the direct operational tau approach ^[7], the spectral tau approach ^[8], the Taylor polynomial scheme ^[9], the Adomian decomposition method ^[10], and the Laguerre matrix scheme ^[11].

Over the years, spectral methods have proven highly effective for approximating solutions to differential and integral equations, particularly when the data and solutions are smooth. These methods use global polynomial approximations, making them especially suitable for problems requiring high accuracy. Since the solution of a delay differential equation (DDE) globally depends on its history due to the delay variable, a global spectral method presents a good approach for solving DDEs by capturing the entire solution behavior with high precision across the domain. For instance, Rahimkhani et al. ^[3] proposed a numerical method for solving fractional pantograph differential equations using the Müntz-Legendre wavelet approximation. This approach emphasized the operational matrix of fractional-order integration to reduce the differential equation to a system of algebraic equations, facilitating efficient numerical solutions. Ezz-Eldien and Doha ^[12] proposed a numerical method for solving pantograph Volterra integro-differential equations using a Chebyshev collocation technique. Yang and Hou ^[14] extended the Jacobi spectral method to nonlinear fractional pantograph differential equations by transforming them into Volterra-Fredholm integral equations. Ghomanjani et al. ^[15] introduced a Bezier curve-based approach for solving Volterra delay-integro-differential equations and linear systems with distributed input delay and saturation. Their method involved a two-step approximation strategy, where the time interval was divided into subintervals, and piecewise Bezier polynomials were applied. Huang et al. ^[16] utilized the discontinuous Galerkin method for the numerical treatment of pantograph-type differential equations. Qin et al. ^[17] developed an $hp$ -version fractional spectral collocation method for weakly singular Volterra integro-differential equations with vanishing delay, demonstrating that approximations in a fractional polynomial space can achieve exponentially decreasing errors even in the presence of weak singularities. Zaky et al. ^{[18,19,20,21]} formulated and analyzed spectral collocation schemes for solving broad classes of linear/nonlinear integral/differential equations. Rahimkhani et al. ^[22] introduced the Hahn wavelet method for fractional-order integro-differential equations, transforming them into integer-order forms using the Laplace transform. This approach allows for continuous and differentiable solutions without requiring operational matrices.

This study enhances the numerical treatment of high-order P-VF-IDEs by addressing a literature gap through the following contributions. The Legendre spectral tau method is introduced to simultaneously manage high-order derivatives, multiple pantograph delays, and mixed Volterra-Fredholm integrals within a unified framework. While existing methods, such as Chebyshev collocation ^[12,13], Jacobi spectral techniques ^[14], or Müntz-Legendre wavelets ^[3], focus on specific aspects (e.g., Volterra-only integrals or one-dimensional cases), The proposed framework unifies these elements, tackling equations previously considered intractable due to their combined complexity. By integrating high-order derivatives, multiple delays, and mixed integrals, this approach goes beyond the disjointed methods found in current literature, establishing a basis for future advancements in fractional or partial P-VF-IDEs.

This paper is organized as follows: In the next section, we provide some necessary properties of Legendre polynomial approximations. In Section 3, we construct the Legendre tau scheme for solving the one-dimensional high-order P-VF-IDEs. In Section 4, we extend the numerical scheme for the two-dimensional case. Section 5 provides some numerical results to clarify the schemes. The conclusion is given in the last section.

2. Legendre polynomials

The shifted Legendre polynomials $P^\sigma_i(x)$ , defined on the interval $\Omega = [0, \sigma]$ , are the eigenfunctions of the Sturm-Liouville problem:

$\frac{d}{dx} \left( x(\sigma - x) \frac{d}{dx} P^\sigma_i(x) \right) + i(i + 1) P^\sigma_i(x) = 0, \quad x \in \Omega,$

and can be expressed as:

$P^\sigma_i(x) = \sum\limits_{p = 0}^i C_{i, p} \frac{x^p}{\sigma^p} = \sum\limits_{p = 0}^i (-1)^{i-p} \frac{(i + p)!}{(i - p)! (p!)^2 \sigma^p} x^p,$

where $C_{i, p} = (-1)^{i-p} \frac{(i + p)!}{(i - p)! (p!)^2 \sigma^p}$ . These polynomials satisfy the derivative conditions:

$\frac{d^s}{dx^s} P^\sigma_i(0) = (-1)^{i-s} \frac{(i + s)!}{(i - s)! s! \sigma^s}, \quad \text{and} \quad \frac{d^s}{dx^s} P^\sigma_i(\sigma) = \frac{(i + s)!}{(i - s)! s! \sigma^s}.$

The shifted Legendre polynomials satisfy the orthogonality relation:

$\begin{equation} \int_0^\sigma P^\sigma_i(x) P^\sigma_p(x) \, dx = \frac{\sigma}{2i + 1} \delta_{ip}, \end{equation}$

(2.1)

where $\delta_{ip}$ is the Kronecker delta. Any function $y \in L^2(\Omega)$ can be expanded as:

$\begin{equation} y(x) = \sum\limits_{i = 0}^\infty a_i P^\sigma_i(x), \quad \text{where} \quad a_i = \frac{2i + 1}{\sigma} \int_0^\sigma y(x) P^\sigma_i(x) \, dx. \end{equation}$

(2.2)

Let $\text{Proj}^\sigma_{\mathcal{N}}$ denote the orthogonal projection operator:

$\text{Proj}^\sigma_{\mathcal{N}}: L^2(\Omega) \rightarrow S_{\mathcal{N}}, \quad S_{\mathcal{N}} = \text{Span}\{ P^\sigma_i(x) : 0 \leq i \leq \mathcal{N} \}.$

This projection is given by:

$\begin{equation} \text{Proj}^\sigma_{\mathcal{N}} y = y_{\mathcal{N}}(x) = \sum\limits_{i = 0}^{\mathcal{N}} a_i P^\sigma_i(x) = A_{\mathcal{N}}^T \Phi^\sigma_{\mathcal{N}}(x), \end{equation}$

(2.3)

where:

$\begin{equation} A_{\mathcal{N}} = \left[ a_i \quad \text{for} \quad 0 \leq i \leq \mathcal{N} \right]^T, \quad \Phi^\sigma_{\mathcal{N}}(x) = \left[ P^\sigma_i(x) \quad \text{for} \quad 0 \leq i \leq \mathcal{N} \right]^T. \end{equation}$

(2.4)

For functions $y \in L^2(\Omega_1 \times \Omega_2)$ with $\Omega_1 = [0, \sigma_1]$ and $\Omega_2 = [0, \sigma_2]$ , the approximation becomes:

$\begin{equation} y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) = \sum\limits_{i = 0}^{\mathcal{N}_1} \sum\limits_{j = 0}^{\mathcal{N}_2} a_{i, j} P^{\sigma_1}_i(x_1) P^{\sigma_2}_j(x_2) = A_{\mathcal{N}_1, \mathcal{N}_2}^T \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \end{equation}$

(2.5)

where:

$\begin{equation} \begin{split} A_{\mathcal{N}_1, \mathcal{N}_2} & = \left[ a_{i, j} \quad \text{for} \quad 0 \leq i \leq \mathcal{N}_1, \, 0 \leq j \leq \mathcal{N}_2 \right]^T, \\ \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) & = \left[ P^{\sigma_1}_i(x_1) P^{\sigma_2}_j(x_2) \quad \text{for} \quad 0 \leq i \leq \mathcal{N}_1, \, 0 \leq j \leq \mathcal{N}_2 \right]^T. \end{split} \end{equation}$

(2.6)

The coefficients $a_{i, j}$ are computed via:

$\begin{equation} a_{i, j} = \frac{(2i + 1)(2j + 1)}{\sigma_1 \sigma_2} \int_0^{\sigma_1} \int_0^{\sigma_2} y(x_1, x_2) P^{\sigma_1}_i(x_1) P^{\sigma_2}_j(x_2) \, dx_1 dx_2. \end{equation}$

(2.7)

3. One-dimensional case

The spectral tau approach for (1.1) is to find $y_{\mathcal{N}}\in P_{\mathcal{N}}$ , such that

$\begin{equation} \begin{split} \sum\limits_{k = 0}^{n}b_k y_{\mathcal{N}}^{(k)}(x) = g_{\mathcal{N}}(x)+\sum\limits_{k = 0}^{n}c_k y_{\mathcal{N}}^{(k)}(p_k x)+\sum\limits_{k = 0}^{n}\int_{0}^{p_k x}d_k y_{\mathcal{N}}^{(k)}(s)ds+\int_{0}^{x}e y_{\mathcal{N}}(s)ds+\int_{0}^{\sigma}f y_{\mathcal{N}}(s)ds. \end{split} \end{equation}$

(3.1)

Now, we denote

$\begin{equation} \begin{split} y_{\mathcal{N}}(x)& = A^T_\mathcal{N} \Phi^{\sigma}_{\mathcal{N}}(x), \\ g_{\mathcal{N}}(x)& = \mathcal{G}^T_\mathcal{N} \Phi^{\sigma}_{\mathcal{N}}(x), \\ \end{split} \end{equation}$

(3.2)

where $\mathcal{G}_\mathcal{N}$ is a known vector given by

$\mathcal{G}_\mathcal{N} = \left[g_0, g_1, \cdots, g_\mathcal{N}\right]^T; \qquad g_j = \frac{2j+1}{\sigma} \int_{0}^{\sigma} g(x) P^{\sigma}_{j}(x) dx.$

If we denote $x_{\mathcal{N}, j, \sigma}$ by the Legendre Gauss quadrature nodes on the $(0, \sigma)$ and $\varpi_{\mathcal{N}, j, \sigma}$ by its corresponding Christoffel numbers, then we have

$g_j\simeq \frac{2j+1}{\sigma} \sum\limits_{k = 0}^{\mathcal{N}} \varpi_{\mathcal{N}, k, \sigma}g(x_{\mathcal{N}, k, \sigma}) P^{\sigma}_{j}(x_{\mathcal{N}, k, \sigma}).$

The following theorems are of great benefit later.

Theorem 1. ^[23] The derivative of order $k$ for the vector $\Phi^{\sigma}_{\mathcal{N}}(x)$ is given by:

$\begin{equation} \frac{d^k}{d x^k}\Phi^{\sigma}_{\mathcal{N}}(x) = \boldsymbol{D}_{(\mathcal{N}, k)}\Phi^{\sigma}_{\mathcal{N}}(x); \qquad \boldsymbol{D}_{(\mathcal{N}, k)} = \left(\boldsymbol{D}^{(1)}\right)^k, \end{equation}$

(3.3)

where $\boldsymbol{D}^{(1)} = (d_{r, j})_{0\leq r, j \leq \mathcal{N}},$ and

$d_{r, j} = \begin{cases}\dfrac{2}{\sigma}(2j+1), \ & r = j+s, \ \begin{cases}s = 1, 3, \cdots, \mathcal{N}-1, \ &\mathcal{N} \;\mathit{\text{is even}}, \\ s = 1, 3, \cdots, \mathcal{N}, \ &\mathcal{N} \;\mathit{\text{is odd}}, \end{cases}\\ 0, &otherwise. \end{cases}$

Theorem 2. ^[24] The integration of the vector $\Phi^{\sigma}_{\mathcal{N}}(x)$ is given by:

$\begin{equation} \int_0^{x} \Phi^{\sigma}_{\mathcal{N}}(x)dx = \boldsymbol{I}_{(\mathcal{N}, 1)} \Phi^{\sigma}_{\mathcal{N}}(x), \end{equation}$

(3.4)

where

$\boldsymbol{I}_{(\mathcal{N}, 1)} = (\beta_{i, j})_{0\leq i, j\leq \mathcal{N}}; \qquad \quad \beta_{i, j} = \sum\limits_{k = 0}^i\sum\limits_{l = 0}^{j}\frac{(-1)^{i+j+k+l} (i+k)!(j+l)!\sigma}{(i-k)!k!(k+1)!(j-l)!(l!)^2(k+l+2)}.$

Theorem 3. For $0 \leq p \leq 1,$ the pantograph operational matrix $\boldsymbol{Q}_{\mathcal{N}, p}$ can be defined by

$\begin{equation} \Phi^{\sigma}_{\mathcal{N}}(px) = \boldsymbol{Q}_{\mathcal{N}, p} \Phi^{\sigma}_{\mathcal{N}}(x), \end{equation}$

(3.5)

where $\boldsymbol{Q}_{\mathcal{N}, p} = (q^p_{k, j})_{0\leq k, j\leq \mathcal{N}},$ and $q^p_{k, j} = \sum\limits_{r = 0}^{k} \dfrac{C_{k, r}p^rf_{r, j}}{\sigma^r}$ .

Proof. We start by expressing $P^{\sigma}_{k}(px)$ by:

$\begin{equation} P^{\sigma}_{k}(px) = \sum\limits_{r = 0}^{k} C_{k, r} \frac{p^rx^r}{\sigma^r}. \end{equation}$

(3.6)

Expanding $x^r$ in terms of $P^{\sigma}_{j}(x), \ j = 0, 1, \cdots, \mathcal{N},$ by

$\begin{equation} x^r = \sum\limits_{j = 0}^{\mathcal{N}}f_{r, j}P^{\sigma}_{j}; \qquad \qquad f_{r, j} = \frac{2j+1}{\sigma} \int_{0}^{\sigma} x^r P^{\sigma}_{j}(x) dx, \qquad 0\leq j. \end{equation}$

(3.7)

A combination of (3.6) and (3.7) then yields

$\begin{equation} \begin{split} P^{\sigma}_{k}(px)& = \sum\limits_{r = 0}^{k} \frac{C_{k, r}p^r}{\sigma^r}\left(\sum\limits_{j = 0}^{\mathcal{N}}f_{r, j}P^{\sigma}_{j}\right) = \sum\limits_{j = 0}^{\mathcal{N}}P^{\sigma}_{j}\left(\sum\limits_{r = 0}^{k} \frac{C_{k, r}p^rf_{r, j}}{\sigma^r}\right)\\ & = \left[\sum\limits_{r = 0}^{k} \frac{C_{k, r}p^rf_{r, 0}}{\sigma^r}, \sum\limits_{r = 0}^{k} \frac{C_{k, r}p^rf_{r, 1}}{\sigma^r}, \cdots, \sum\limits_{r = 0}^{k} \frac{C_{k, r}p^rf_{r, \mathcal{N}}}{\sigma^r} \right]^T \Phi^{\sigma}_{\mathcal{N}}(x), \end{split} \end{equation}$

(3.8)

which completes the proof. □

Application of Theorems 1–3, we have

$\begin{equation} \begin{split} \frac{d^k y_{\mathcal{N}}(x)}{dx^k} & = A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \Phi^{\sigma}_{\mathcal{N}}(x), \\ \int_{0}^{x} y_{\mathcal{N}}(s)ds & = A^T_\mathcal{N} \textbf{I}_{(\mathcal{N}, 1)} \Phi^{\sigma}_{\mathcal{N}}(x), \\ y_{\mathcal{N}}(px)& = A^T_{\mathcal{N}}\textbf{Q}_{\mathcal{N}, p} \Phi^{\sigma}_{\mathcal{N}}(x), \\ \frac{d^k y_{\mathcal{N}}(p x)}{dx^k}& = A^T_{\mathcal{N}}\textbf{D}_{(\mathcal{N}, k)}\textbf{Q}_{\mathcal{N}, p} \Phi^{\sigma}_{\mathcal{N}}(x), \\ \int_{0}^{x} y_{\mathcal{N}}(p s)ds & = A^T_\mathcal{N} \textbf{I}_{(\mathcal{N}, 1)} \textbf{Q}_{\mathcal{N}, p} \Phi^{\sigma}_{\mathcal{N}}(x). \end{split} \end{equation}$

(3.9)

Using Theorems 1–3, and with the help of (3.9), the residual $\mathcal{R}_{\mathcal{N}}(x)$ for (3.1) is given by:

$\begin{equation} \begin{split} \mathcal{R}_{\mathcal{N}}(x)& = \sum\limits_{k = 0}^{n}b_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \Phi^{\sigma}_{\mathcal{N}}(x)-\sum\limits_{k = 0}^{n}c_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \textbf{Q}_{\mathcal{N}, p_k} \Phi^{\sigma}_{\mathcal{N}}(x)-e A^T_\mathcal{N} \textbf{I}_{(\mathcal{N}, 1)} \Phi^{\sigma}_{\mathcal{N}}(x)\\ &-\sum\limits_{k = 0}^{n}d_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \textbf{I}_{(\mathcal{N}, 1)} \textbf{Q}_{\mathcal{N}, p_k} \Phi^{\sigma}_{\mathcal{N}}(x)-f a_0\sigma-\mathcal{G}^T_\mathcal{N} \Phi^{\sigma}_{\mathcal{N}}(x), \end{split} \end{equation}$

(3.10)

According to the spectral tau method, the numerical solution for (1.1) is obtained by generating $\mathcal{N}+1$ linear algebraic equations as follows:

$\begin{equation} \begin{split} \int_{0}^{\sigma} \mathcal{R}_{\mathcal{N}}(x) P^{\sigma}_{k}(x) dx& = 0, \qquad \qquad k = 0, 1, \cdots, \mathcal{N}-n, \\ \end{split} \end{equation}$

(3.11)

$\begin{equation} \begin{split} A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \Phi^{\sigma}_{\mathcal{N}}(0)& = \varrho_k, \qquad \qquad k = 0, 1, \cdots, n-1. \end{split} \end{equation}$

(3.12)

Define the vector ${}_{\sigma}\textbf{M}_{\mathcal{N}}^i, \ i = 0, 1, \cdots, \mathcal{N}-n,$ as follows:

${}_{\sigma}\textbf{M}_{\mathcal{N}}^i = \frac{\sigma}{2i+1}\mathbf{e}_i,$

where $\mathbf{e}_i$ is the $\mathcal{N}+1$ standard basis vector. Then the system of equations (3.11) simplifies to:

$\begin{equation} \begin{split} \sum\limits_{k = 0}^{n}b_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} & {}_{\sigma}\textbf{M}_{\mathcal{N}}^i-\sum\limits_{k = 0}^{n}c_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \textbf{Q}_{\mathcal{N}, p_k} {}_{\sigma}\textbf{M}_{\mathcal{N}}^i -\sum\limits_{k = 0}^{n}d_k A^T_\mathcal{N} \textbf{D}_{(\mathcal{N}, k)} \textbf{I}_{(\mathcal{N}, 1)} \textbf{Q}_{\mathcal{N}, p_k} {}_{\sigma}\textbf{M}_{\mathcal{N}}^i\\ & -e A^T_\mathcal{N} \textbf{I}_{(\mathcal{N}, 1)} {}_{\sigma}\textbf{M}_{\mathcal{N}}^i-f \sigma^2 \delta_{i0} A^T_\mathcal{N} \mathbf{e}_0 -\mathcal{G}^T_\mathcal{N} {}_{\sigma}\textbf{M}_{\mathcal{N}}^i = 0, \qquad i = 0, 1, \cdots, \mathcal{N}-n, \end{split} \end{equation}$

(3.13)

If we denote $E_i, \ i = 0, 1, \cdots, \mathcal{N}-n,$ by

$\begin{split} E_i = \sum\limits_{k = 0}^{n}b_k \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^i\right)^T \textbf{D}_{(\mathcal{N}, k)}^T & -\sum\limits_{k = 0}^{n} c_k \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^i\right)^T \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{D}_{(\mathcal{N}, k)}^T -\sum\limits_{k = 0}^{n}d_k \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^i\right)^T \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{I}_{(\mathcal{N}, 1)}^T \textbf{D}_{(\mathcal{N}, k)}^T \\ & -e \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^i\right)^T \textbf{I}_{(\mathcal{N}, 1)}^T -f \delta_{i0} \sigma^2 \mathbf{e}^T_0, \qquad i = 0, 1, \cdots, \mathcal{N}-n, \end{split}$

and

$C_i = \left(\Phi^{\sigma}_{\mathcal{N}}(0)\right)^T\textbf{D}_{(\mathcal{N}, i)}^T, \qquad i = 0, 1, \cdots, n-1,$

then the solution of the main problem is reduced to the system

$\textbf{S} A_{\mathcal{N}} = \textbf{B}.$

The coefficient matrix S is expressed as:

$\textbf{S} = \begin{bmatrix} \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^0\right)^T \left( \sum\nolimits_{k = 0}^{n} b_k \textbf{D}_{(\mathcal{N}, k)}^T - c_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{D}_{(\mathcal{N}, k)}^T - d_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{I}_{(\mathcal{N}, 1)}^T \textbf{D}_{(\mathcal{N}, k)}^T - e \textbf{I}_{(\mathcal{N}, 1)}^T \right) -f \sigma^2 \mathbf{e}_0^T \\ \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^1\right)^T \left( \sum\nolimits_{k = 0}^{n} b_k \textbf{D}_{(\mathcal{N}, k)}^T - c_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{D}_{(\mathcal{N}, k)}^T - d_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{I}_{(\mathcal{N}, 1)}^T \textbf{D}_{(\mathcal{N}, k)}^T - e \textbf{I}_{(\mathcal{N}, 1)}^T \right) \\ \vdots \\ \left({}_{\sigma}\textbf{M}_{\mathcal{N}}^{\mathcal{N}-n}\right)^T \left( \sum\nolimits_{k = 0}^{n} b_k \textbf{D}_{(\mathcal{N}, k)}^T - c_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{D}_{(\mathcal{N}, k)}^T - d_k \textbf{Q}_{\mathcal{N}, p_k}^T \textbf{I}_{(\mathcal{N}, 1)}^T \textbf{D}_{(\mathcal{N}, k)}^T - e \textbf{I}_{(\mathcal{N}, 1)}^T \right) \\ \left( \Phi^{\sigma}_{\mathcal{N}}(0) \right)^T \textbf{D}_{(\mathcal{N}, 0)}^T \\ \vdots \\ \left( \Phi^{\sigma}_{\mathcal{N}}(0) \right)^T \textbf{D}_{(\mathcal{N}, n-1)}^T \end{bmatrix}.$

Compactly, we can write S as:

$\textbf{S} = \begin{bmatrix} E_0 \\ E_1 \\ \vdots \\ E_{\mathcal{N}-n} \\ C_0 \\ C_1 \\ \vdots \\ C_{n-1} \end{bmatrix},$

where each $E_i$ is a row vector of length $\mathcal{N}+1$ , and the total system is of size $(\mathcal{N}+1)\times (\mathcal{N}+1)$ . The right-hand side vector B is given by:

$\textbf{B} = \begin{bmatrix} \frac{g_0 \sigma}{1} \\ \frac{g_1 \sigma}{3} \\ \frac{g_2 \sigma}{5} \\ \vdots \\ \frac{g_{\mathcal{N}-n} \sigma}{2\mathcal{N}-2n+1} \\ \varrho_0 \\ \varrho_1 \\ \vdots \\ \varrho_{n-1} \end{bmatrix}.$

4. Two-dimensional case

In this section, we employ the Legendre spectral tau method to solve the two-dimensional high-order pantograph-type Volterra-Fredholm integro-differential equation:

$\begin{equation} \begin{split} &\sum\limits_{k_1 = 0}^{n_1}\sum\limits_{k_2 = 0}^{n_2} b_{k_1, k_2} \frac{\partial^{k_1+k_2} y(x_1, x_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} = g(x_1, x_2) + \int_{0}^{\sigma_2} \int_{0}^{\sigma_1} c \, y(s_1, s_2) \, ds_1 ds_2 \\ &+ \int_{0}^{x_2} \int_{0}^{x_1} d \, y(s_1, s_2) \, ds_1 ds_2 + \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} \int_{0}^{p_1 x_1} \int_{0}^{p_2 x_2} e_{k_1, k_2} \frac{\partial^{k_1+k_2} y(s_1, s_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} \, ds_1 ds_2 \\ &+ \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \frac{\partial^{k_1+k_2} y(p_1 x_1, p_2 x_2)}{\partial x_1^{k_1} \partial x_2^{k_2}}, \end{split} \end{equation}$

(4.1)

subject to the initial conditions:

$\begin{equation} \begin{split} \frac{\partial^{k_1} y(x_1, x_2)}{\partial x_1^{k_1}} & = \varrho_{1, k_1}(0, x_2), \quad k_1 = 0, 1, \dots, n_1-1, \quad 0 \leq x_2 \leq \sigma_2, \\ \frac{\partial^{k_2} y(x_1, x_2)}{\partial x_2^{k_2}} & = \varrho_{2, k_2}(x_1, 0), \quad k_2 = 0, 1, \dots, n_2-1, \quad 0 \leq x_1 \leq \sigma_1, \end{split} \end{equation}$

(4.2)

where $\sigma_1, \sigma_2, c, d, b_{k_1, k_2}, e_{k_1, k_2}$ , and $f_{k_1, k_2}$ (for $0 \leq k_1 \leq n_1$ , $0 \leq k_2 \leq n_2$ ) are known real constants.

As a spectral approach, we seek an approximate solution $y_{\mathcal{N}_1, \mathcal{N}_2} \in P_{\mathcal{N}_1} \times P_{\mathcal{N}_2}$ such that:

$\begin{equation} \begin{split} &\sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} b_{k_1, k_2} \frac{\partial^{k_1+k_2} y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} = g_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) \\ &+ \int_{0}^{x_2} \int_{0}^{x_1} d \, y_{\mathcal{N}_1, \mathcal{N}_2}(s_1, s_2) \, ds_1 ds_2 + \int_{0}^{\sigma_2} \int_{0}^{\sigma_1} c \, y_{\mathcal{N}_1, \mathcal{N}_2}(s_1, s_2) \, ds_1 ds_2 \\ &+ \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \frac{\partial^{k_1+k_2} y_{\mathcal{N}_1, \mathcal{N}_2}(p_1 x_1, p_2 x_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} \\ &+ \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} \int_{0}^{p_1 x_1} \int_{0}^{p_2 x_2} e_{k_1, k_2} \frac{\partial^{k_1+k_2} y_{\mathcal{N}_1, \mathcal{N}_2}(s_1, s_2)}{\partial x_1^{k_1} \partial x_2^{k_2}} \, ds_1 ds_2. \end{split} \end{equation}$

(4.3)

We represent the approximate solution $y_{\mathcal{N}_1, \mathcal{N}_2}$ and the right-hand side function $g_{\mathcal{N}_1, \mathcal{N}_2}$ in terms of Legendre basis functions:

$\begin{equation} \begin{split} y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ g_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) & = \mathcal{G}^T_{\mathcal{N}_1, \mathcal{N}_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \end{split} \end{equation}$

(4.4)

where $\mathcal{G}_{\mathcal{N}_1, \mathcal{N}_2}$ is a known vector given by:

$\mathcal{G}_{\mathcal{N}_1, \mathcal{N}_2} = \left[ g_{0, 0}, g_{1, 0}, \dots, g_{\mathcal{N}_1, 0}, g_{0, 1}, \dots, g_{\mathcal{N}_1, \mathcal{N}_2} \right]^T,$

and the coefficients $g_{i, j}$ are computed as:

$g_{i, j} = \frac{(2i+1)(2j+1)}{\sigma_1 \sigma_2} \int_{0}^{\sigma_2} \int_{0}^{\sigma_1} g(x_1, x_2) \, P^{\sigma_1}_i(x_1) \, P^{\sigma_2}_j(x_2) \, dx_1 dx_2.$

The following matrix operators are defined for differentiation, integration, and scaling:

Theorem 4. Let $\mathbb{I}_{\mathcal{N}_1}$ and $\mathbb{I}_{\mathcal{N}_2}$ denote the identity matrices of orders ( $\mathcal{N}_1$ +1) and $(\mathcal{N}_2+1)$ , respectively. Then:

$\begin{equation} \frac{\partial^k}{\partial x_1^k} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) = \mathcal{D}_{(1, k)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{D}_{(1, k)} = \mathbf{D}_{(\mathcal{N}_1, 1)} \otimes \mathbb{I}_{\mathcal{N}_2}, \end{equation}$

(4.5)

$\begin{equation} \frac{\partial^k}{\partial x_2^k} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) = \mathcal{D}_{(2, k)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{D}_{(2, k)} = \mathbb{I}_{\mathcal{N}_1} \otimes \mathbf{D}_{(\mathcal{N}_2, 2)}, \end{equation}$

(4.6)

$\begin{equation} \int_{0}^{x_1} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(s_1, x_2) \, ds_1 = \mathcal{I}_{(1, 1)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{I}_{(1, 1)} = \mathbf{I}_{(\mathcal{N}_1, 1)} \otimes \mathbb{I}_{\mathcal{N}_2}, \end{equation}$

(4.7)

$\begin{equation} \int_{0}^{x_2} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, s_2) \, ds_2 = \mathcal{I}_{(2, 1)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{I}_{(2, 1)} = \mathbb{I}_{\mathcal{N}_1} \otimes \mathbf{I}_{(\mathcal{N}_2, 2)}, \end{equation}$

(4.8)

$\begin{equation} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(p_1 x_1, x_2) = \mathcal{Q}_{1, p_1} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{Q}_{1, p_1} = \mathbf{Q}_{\mathcal{N}_1, p_1} \otimes \mathbb{I}_{\mathcal{N}_2}, \end{equation}$

(4.9)

$\begin{equation} \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, p_2 x_2) = \mathcal{Q}_{2, p_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \quad \mathcal{Q}_{2, p_2} = \mathbb{I}_{\mathcal{N}_1} \otimes \mathbf{Q}_{\mathcal{N}_2, p_2}. \end{equation}$

(4.10)

In virtue of (4.4)–(4.10), we have

$\begin{equation} \begin{split} \frac{\partial^k}{\partial x_1^k} y_{\mathcal{N}_1, \mathcal{N}_2} (x_1, x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{D}_{(1, k)}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ \frac{\partial^k}{\partial x_2^k} y_{\mathcal{N}_1, \mathcal{N}_2} (x_1, x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{D}_{(2, k)}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ \int_{0}^{x_1} y_{\mathcal{N}_1, \mathcal{N}_2} (s_1, x_2) ds_1 & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{I}_{(1, 1)}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ \int_{0}^{x_2} y_{\mathcal{N}_1, \mathcal{N}_2} (x_1, s_2) ds_2 & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{I}_{(2, 1)}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ y_{\mathcal{N}_1, \mathcal{N}_2} (p_1 x_1, x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{Q}_{1, p_1}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ y_{\mathcal{N}_1, \mathcal{N}_2} (x_1, p_2x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{Q}_{2, p_2}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2), \\ y_{\mathcal{N}_1, \mathcal{N}_2} (p_1x_1, p_2x_2) & = A^T_{\mathcal{N}_1, \mathcal{N}_2} \mathcal{Q}_{1, p_2}\mathcal{Q}_{2, p_2}\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2). \\ \end{split} \end{equation}$

(4.11)

Using the basis representation and matrix operators, the residual $\mathcal{R}_{\mathcal{N}_1, \mathcal{N}_2}$ of (4.3) is given by:

$\begin{split} \mathcal{R}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) & = \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} b_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) - \mathcal{G}^T_{\mathcal{N}_1, \mathcal{N}_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2)\\ &- d \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{I}_{(1, 1)} \, \mathcal{I}_{(2, 1)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) - c \sigma_1\sigma_2a_{0, 0} \\ &- \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} e_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \mathcal{I}_{(1, 1)} \, \mathcal{I}_{(2, 1)} \, \mathcal{Q}_{1, p_1} \, \mathcal{Q}_{2, p_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2) \\ &- \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \mathcal{Q}_{1, p_1} \, \mathcal{Q}_{2, p_2} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2). \end{split}$

The system of $(\mathcal{N}_1+1)(\mathcal{N}_2+1)$ algebraic equations is obtained by enforcing the residual to be orthogonal to the basis functions:

$\begin{equation} \begin{split} &\int_{0}^{\sigma_2} \int_{0}^{\sigma_1} \mathcal{R}_{\mathcal{N}_1, \mathcal{N}_2}(x_{1, i}, x_{2, j}) \, P^{\sigma_1}_i(x_1) \, P^{\sigma_2}_j(x_2) \, dx_1 dx_2 = 0, \\ &\quad 0 \leq i \leq \mathcal{N}_1 - n_1, \quad 0 \leq j \leq \mathcal{N}_2 - n_2, \\ \end{split} \end{equation}$

(4.12)

$\begin{equation} \begin{split} &A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(0, x_{2, j}) = \varrho_{1, k_1}(x_{2, j}), \quad 0 \leq k_1 \leq n_1-1, \quad 0 \leq j \leq \mathcal{N}_2, \\ \end{split} \end{equation}$

(4.13)

$\begin{equation} \begin{split} &A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(2, k_2)} \, \Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_{1, j}, 0) = \varrho_{2, k_2}(x_{1, j}), \quad 0 \leq k_2 \leq n_2-1, \quad n_1 \leq j \leq \mathcal{N}_1, \end{split} \end{equation}$

(4.14)

where $x_{1, i}$ and $x_{2, j}$ are the roots of $P^{\sigma_1}_{\mathcal{N}_1+1}(x_1)$ and $P^{\sigma_2}_{\mathcal{N}_2+1}(x_2)$ , respectively.

Define the vector $\mathbf{e}_{00} = \mathbf{e}_0\otimes\mathbf{e}_0$ and ${}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}, \ i = 0, 1, \cdots, \mathcal{N}_1-n_1, \ j = 0, 1, \cdots, \mathcal{N}_2-n_2,$ as follows:

${}_{\sigma_1, \sigma_2} \textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j} = {}_{\sigma_1} \textbf{M}^i_{\mathcal{N}_1}\otimes {}_{\sigma_2} \textbf{M}^j_{\mathcal{N}_2}, \quad i = 0, 1, \cdots, \mathcal{N}_1-n_1, \ j = 0, 1, \cdots, \mathcal{N}_2-n_2,$

then the system (4.12) is simplified to:

$\begin{split} \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} & b_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j} - \mathcal{G}^T_{\mathcal{N}_1, \mathcal{N}_2} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}- d \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{I}_{(1, 1)} \, \mathcal{I}_{(2, 1)} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\\ &- c \sigma_1^2 \sigma_2^2 \delta_{i0} \delta_{j0} A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathbf{e}_{00} - \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \mathcal{Q}_{1, p_1} \, \mathcal{Q}_{2, p_2} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j} \\ &- \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} e_{k_1, k_2} \, A^T_{\mathcal{N}_1, \mathcal{N}_2} \, \mathcal{D}_{(1, k_1)} \, \mathcal{D}_{(2, k_2)} \, \mathcal{I}_{(1, 1)} \, \mathcal{I}_{(2, 1)} \, \mathcal{Q}_{1, p_1} \, \mathcal{Q}_{2, p_2} \, {}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}. \end{split}$

Denoting $E_{i, j}, \ i = 0, 1, \cdots, \mathcal{N}_1-n_1, \ j = 0, 1, \cdots, \mathcal{N}_2-n_2,$ by

$\begin{split} E_{i, j}& = \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} b_{k_1, k_2} \, \left({}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\right)^T \, \mathcal{D}^T_{(2, k_2)} \, \mathcal{D}^T_{(1, k_1)} - d \, \left({}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\right)^T \mathcal{I}^T_{(2, 1)} \mathcal{I}^T_{(1, 1)}\\ &- c \, \sigma_1^2 \sigma_2^2 \delta_{i0} \delta_{j0} \mathbf{e}_{00}^T - \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} f_{k_1, k_2} \, \left({}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\right)^T \mathcal{Q}^T_{2, p_2} \mathcal{Q}^T_{1, p_1} \mathcal{D}^T_{(2, k_2)}\mathcal{D}^T_{(1, k_1)} \\ &- \sum\limits_{k_1 = 0}^{n_1} \sum\limits_{k_2 = 0}^{n_2} e_{k_1, k_2} \, \left({}_{\sigma_1, \sigma_2}\textbf{M}_{\mathcal{N}_1, \mathcal{N}_2}^{i, j}\right)^T \mathcal{Q}^T_{2, p_2} \mathcal{Q}^T_{1, p_1} \mathcal{I}^T_{(2, 1)} \mathcal{I}^T_{(1, 1)} \mathcal{D}^T_{(2, k_2)} \mathcal{D}^T_{(1, k_1)} , \end{split}$

and

$\begin{split} &C_{1, i, j} = \left(\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(0, x_{2, j}) \right)^T\mathcal{D}^T_{(1, i)}, \qquad i = 0, 1, \cdots, n_1-1, \quad j = 0, 1, \cdots, \mathcal{N}_2, \\ &C_{2, i, j} = \left(\Phi^{\sigma_1, \sigma_2}_{\mathcal{N}_1, \mathcal{N}_2}(x_{1, j}, 0)\right)^T\mathcal{D}^T_{(2, i)}, \qquad i = 0, 1, \cdots, n_2-1, \quad j = n_1, n_1+1, \cdots, \mathcal{N}_1. \end{split}$

The solution of the main problem is reduced to the system

$\textbf{S} A_{\mathcal{N}_1, \mathcal{N}_2} = \textbf{B},$

where

$\begin{split} \textbf{S} = \Big[E_{0, 0};\ E_{0, 1}&;\ \cdots;\ E_{\mathcal{N}_1-n_1;\mathcal{N}_2-n_2};\ C_{1, 0, 0};\ C_{1, 0, \mathcal{N}_2};\ \cdots;\ C_{1, n_1-1, \mathcal{N}_2};\ C_{2, 0, n_1};\ C_{2, 0, n_1+1};\ \cdots;\ C_{2, n_2-1, \mathcal{N}_1}\Big]^T, \\ \textbf{B} = \Big[\frac{g_{0, 0}\sigma_1\sigma_2}{(1)(1)}&, \ \frac{g_{0, 1}\sigma_1\sigma_2}{(1)(3)}, \frac{g_{0, 2}\sigma_1\sigma_2}{(1)(5)}, \ \cdots, \frac{g_{\mathcal{N}_1-n_1, \mathcal{N}_2-n_2}\sigma_1\gamma_2}{(2\mathcal{N}_1-2n_1+1)(2\mathcal{N}_2-2n_2+1)} , \varrho_{1, 0}(x_{2, 0}), \\ & \varrho_{1, 0}(x_{2, 1}), \ \cdots, \ \varrho_{1, n_1-1}(x_{2, \mathcal{N}_2}), \ \varrho_{2, 0}(x_{1, n_1}), \varrho_{2, 0}(x_{1, n_1+1}), \cdots, \ \varrho_{2, n_2-1}(x_{1, \mathcal{N}_1]}) \Big]^T. \end{split}$

5. Numerical results

5.1. Convergence test

In this subsection, we analyze the convergence behavior of the proposed numerical method. We consider a test problem with an irregular solution and compute the maximum absolute errors (MAEs) for different values of $\mathcal{N}$ .

We consider the following P-VF-IDE:

$\begin{equation} \frac{d^2 y(x)}{dx^2} + \frac{dy(x)}{dx} = \frac{1}{2} y\left(\frac{1}{3}x\right) + \int_{0}^{x} y(t) \, dt + \int_{0}^{\frac{1}{3}x} y(t) \, dt + \int_{0}^{1} y(t) \, dt + g(x), \quad 0 \leq x \leq 1, \end{equation}$

(5.1)

with the initial conditions:

$y(0) = 0, \quad \frac{dy(0)}{dx} = 0.$

The function $g(x)$ is chosen such that the exact solution is:

$y(x) = x^{\frac{11}{3}}.$

We apply the proposed numerical scheme to this problem for various values of $\mathcal{N}$ . The MAEs of the numerical solution $y_{\mathcal{N}}(x)$ are computed and summarized in . The MAEs decrease as $\mathcal{N}$ increases, demonstrating the method's high accuracy. The convergence behavior is further visualized in . The results demonstrate that the proposed numerical method achieves excellent convergence for problems with irregular solutions. The decay of errors with increasing $\mathcal{N}$ highlights the efficiency and robustness of the method.

Table 1. The MAEs of

$y_{\mathcal{N}}(x)$ for problem 5.1.

$\mathcal{N}$	4	8	12	16	20
MAE	$8.2705\times10^{-4}$	$4.9139\times10^{-6}$	$6.0563\times10^{-7}$	$1.5001\times10^{-7}$	$4.9759\times10^{-8}$

| Show Table

DownLoad: CSV

Figure 1. Convergence of the approximate solution for problem 5.1.

DownLoad: Full-Size Img PowerPoint

5.2. The numerical stability test

To numerically investigate the stability of the spectral tau method, we consider the following Volterra integro-differential equation ^[25]:

$\begin{equation} \begin{split} \frac{d y(x)}{dx} = y(0.5x) + \int_{0}^{x} y(t) \, dt + \int_{0}^{0.5x} y(t) \, dt + 1 - 1.5x, \quad 0 \leq x \leq 1, \end{split} \end{equation}$

(5.2)

with the initial condition $y(0) = 0$ and the exact solution $y(x) = 1 - e^{-x}$ .

To assess the stability of the method, we introduce perturbations to the right-hand side and the initial condition. Specifically, we consider the following perturbed problems:

$(i)$ Perturbed right-hand side: The perturbed right-hand side problem is given by:

$\begin{equation} \frac{d z(x)}{dx} = z(0.5x) + \int_{0}^{x} z(t) \, dt + \int_{0}^{0.5x} z(t) \, dt + 1 - 1.5x + \epsilon_r, \end{equation}$

(5.3)

with the initial condition $z(0) = 0$ and the exact solution $z(x) = 1 - e^{-x}$ , where $\epsilon_r$ is a small perturbation parameter. This perturbation tests the sensitivity of the method to changes in the forcing term of the equation.

The maximum absolute errors, $\lvert y_{\mathcal{N}}(x) - z_{\mathcal{N}}(x) \rvert$ , for the perturbed right-hand side problem (5.3) are computed for several values of $\epsilon_r$ . The results are summarized in Table 2.

Table 2. Maximum absolute errors for the perturbed right-hand side problem (5.3).

$\mathcal{N}$	$\epsilon_r=0.1$	$\epsilon_r=0.01$	$\epsilon_r=0.001$
6	0.15237708	0.01523770	0.00152376
8	0.15237709	0.01523770	0.00152377
10	0.15237709	0.01523770	0.00152377

| Show Table

DownLoad: CSV

$(ii)$ Perturbed initial condition: The perturbed initial condition problem is given by:

$\begin{equation} \frac{d z(x)}{dx} = z(0.5x) + \int_{0}^{x} z(t) \, dt + \int_{0}^{0.5x} z(t) \, dt + 1 - 1.5x, \end{equation}$

(5.4)

with the perturbed initial condition:

$z(0) = \epsilon_i,$

where $\epsilon_i$ is a small perturbation parameter. This perturbation tests the sensitivity of the method to changes in the initial condition.

The maximum absolute errors, $\lvert y_{\mathcal{N}}(x) - z_{\mathcal{N}}(x) \rvert$ , for the perturbed initial condition problem (5.4) are computed for several values of $\epsilon_i$ . The results are summarized in Table 3.

Table 3. Maximum absolute errors for the perturbed initial condition problem (5.4).

$\mathcal{N}$	$\epsilon_i=0.1$	$\epsilon_i=0.01$	$\epsilon_i =0.001$
6	0.34154212	0.03415420	0.00341541
8	0.34154212	0.03415421	0.00341542
10	0.34154212	0.03415421	0.00341542

| Show Table

DownLoad: CSV

The numerical results demonstrate that the spectral tau method remains stable under small perturbations to both the right-hand side and the initial condition. The errors introduced by the perturbations are proportional to the perturbation parameters $\epsilon_r$ and $\epsilon_i$ , indicating that the method is robust and well-conditioned. Specifically:

● For the perturbed right-hand side, the errors decrease linearly with $\epsilon_r$ .

● For the perturbed initial condition, the errors decrease linearly with $\epsilon_i$ .

These findings confirm the stability of the spectral tau method and its ability to handle small perturbations without significant loss of accuracy.

5.3. Comparison test

In this subsection, we evaluate the performance of the proposed numerical method by comparing it with existing methods. Specifically, we consider the following Volterra integro-differential equation ^[25]:

$\begin{equation} \begin{split} \frac{d y(x)}{dx} = y(0.5x) + \int_{0}^{x} y(t) \, dt + \int_{0}^{0.5x} y(t) \, dt + 1 - 1.5x, \quad 0 \leq x \leq 1, \end{split} \end{equation}$

(5.5)

with the initial condition $y(0) = 0$ and the exact solution $y(x) = 1 - e^{-x}$ .

For the solution of this problem, Zhao et al. ^[25] used the spectral collocation approach based on the Sinc function (SCSFM). This method reduces the problem to solving a system of algebraic equations. In , we compare the MAEs of $y_{\mathcal{N}}(x)$ obtained using our method with those obtained using the SCSFM ^[25].

Table 4. Comparing MAEs of

$y_{\mathcal{N}}(x)$ given using our method against those given using the SCSFM ^[25] for 5.5.

SCSFM ^[25]		Our scheme
$N$	MAE	$M$	MAE
10	$1.2608\times10^{-4}$	2	$7.1205\times10^{-3}$
20	$2.2598\times10^{-6}$	4	$1.6006\times10^{-5}$
30	$9.5984\times10^{-8}$	6	$2.2379\times10^{-8}$
40	$7.1028\times10^{-9}$	8	$1.9964\times10^{-11}$
50	$7.4603\times10^{-10}$	10	$1.0957\times10^{-14}$
60	$9.7609\times10^{-11}$	12	$4.0571\times10^{-16}$

| Show Table

DownLoad: CSV

The results in demonstrate that our method achieves significantly higher accuracy compared to the SCSFM, especially for larger values of $\mathcal{N}$ . For instance, with $\mathcal{N} = 12$ , our method achieves an MAE of $4.0571 \times 10^{-16}$ , whereas the SCSFM achieves $9.7609 \times 10^{-11}$ with $\mathcal{N} = 60$ . This highlights the superior convergence properties of our method.

To further illustrate the accuracy of our method, we plot the absolute error function $|y(x) - y_{\mathcal{N}}(x)|$ for $\mathcal{N} = 6, 10$ and $14$ in Figure 2.

Figure 2. Absolute errors with

$\mathcal{N} = 6, 10,$ and

$14$ for problem 5.5.

DownLoad: Full-Size Img PowerPoint

5.4. Two-dimensional case

In this subsection, we extend the application of the proposed numerical method to a two-dimensional problem. The goal is to demonstrate the effectiveness of the method in solving higher-dimensional integro-differential equations. We consider the following problem:

$\begin{equation} \begin{split} \frac{\partial^{5} y(x_1, x_2)}{\partial x_1^{3} \partial x_2^{2}} & = -y(x_1, x_2) + y(0.5x_1, 0.25x_2) - \int_{0}^{0.25x_2} \int_{0}^{0.5x_1} \frac{\partial y(s_1, s_2)}{\partial x_1} \, ds_1 ds_2 \\ &+ \int_{0}^{x_2} \int_{0}^{x_1} y(s_1, s_2) \, ds_1 ds_2 + \int_{0}^{1} \int_{0}^{1} y(s_1, s_2) \, ds_1 ds_2 + g(x_1, x_2), \end{split} \end{equation}$

(5.6)

where $0 \leq x_1, x_2 \leq 1$ , and the initial conditions are given by:

$y(0, x_2) = \log(x_2 + 1), \quad \frac{\partial y(0, x_2)}{\partial x_1} = \frac{\partial^2 y(0, x_2)}{\partial x_1^2} = y(x_1, 0) = 0, \quad \frac{\partial y(x_1, 0)}{\partial x_2} = x_1^4 + 1.$

The function $g(x_1, x_2)$ is chosen such that the exact solution is:

$y(x_1, x_2) = (x_1^4 + 1) \log(1 + x_2).$

We apply the numerical scheme presented in Section 4 to solve this problem. The discretization is performed with $\mathcal{N}_1 = 4$ and varying values of $\mathcal{N}_2$ . The absolute errors of the numerical solution $y_{4, \mathcal{N}_2}(x_1, x_2)$ are computed for $\mathcal{N}_2 = 4, 8, 12, 16,$ and $20$ . The results are summarized in Table 5.

Table 5. Absolute errors of

$y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2)$ at

$\mathcal{N}_1 = 4$ and different choices of

$\mathcal{N}_2$ for problem 5.6.

$(x_1, x_2)$	$\mathcal{N}_2=4$	$\mathcal{N}_2=8$	$\mathcal{N}_2=12$	$\mathcal{N}_2=16$	$\mathcal{N}_2=20$
(0.1, 0.1)	$5.695\times10^{-5}$	$2.231\times10^{-8}$	$7.399\times10^{-12}$	$3.955\times10^{-15}$	$4.163\times10^{-17}$
(0.2, 0.2)	$1.763\times10^{-5}$	$7.956\times10^{-8}$	$1.358\times10^{-10}$	$1.518\times10^{-13}$	$3.330\times10^{-16}$
(0.3, 0.3)	$2.306\times10^{-4}$	$3.891\times10^{-7}$	$5.919\times10^{-10}$	$7.792\times10^{-13}$	$1.111\times10^{-15}$
(0.4, 0.4)	$6.320\times10^{-4}$	$1.224\times10^{-6}$	$1.783\times10^{-9}$	$2.277\times10^{-12}$	$2.997\times10^{-15}$
(0.5, 0.5)	$1.367\times10^{-3}$	$2.739\times10^{-6}$	$4.033\times10^{-9}$	$5.164\times10^{-12}$	$6.772\times10^{-15}$
(0.6, 0.6)	$2.603\times10^{-3}$	$5.225\times10^{-6}$	$7.698\times10^{-9}$	$9.859\times10^{-12}$	$1.276\times10^{-14}$
(0.7, 0.7)	$4.467\times10^{-3}$	$8.894\times10^{-6}$	$1.303\times10^{-8}$	$1.664\times10^{-11}$	$2.164\times10^{-14}$
(0.8, 0.8)	$7.000\times10^{-3}$	$1.369\times10^{-5}$	$2.004\times10^{-8}$	$2.563\times10^{-11}$	$3.330\times10^{-14}$
(0.9, 0.9)	$1.007\times10^{-2}$	$1.953\times10^{-5}$	$2.860\times10^{-9}$	$3.654\times10^{-11}$	$4.707\times10^{-14}$

| Show Table

DownLoad: CSV

The results in demonstrate the high accuracy of the proposed method for two-dimensional problems. As $\mathcal{N}_2$ increases, the absolute errors decrease significantly, indicating the convergence of the numerical solution to the exact solution. To further illustrate the accuracy of the method, we plot the absolute error function $|y(x_1, x_2) - y_{\mathcal{N}_1, \mathcal{N}_2}(x_1, x_2)|$ for $(\mathcal{N}_1, \mathcal{N}_2) = \{(4, 10), (4, 14), (4, 18)\}$ in . To analyze the convergence behavior, we plot the logarithmic function of the MAEs for $\mathcal{N}_1 = 4$ and varying $\mathcal{N}_2$ in Figure 4.

Figure 3. Absolute error function of

$y(x_1, x_2)$ with

$(\mathcal{N}_1, \mathcal{N}_2) = \{(4, 10), (4, 14), (4, 18)\}$ for problem 5.6.

DownLoad: Full-Size Img PowerPoint

Figure 4. The logarithmic function of MAEs of

$y(x_1, x_2)$ with

$\mathcal{N}_1 = 4$ for problem 5.6.

DownLoad: Full-Size Img PowerPoint

The numerical results confirm the effectiveness of the proposed method for solving two-dimensional integro-differential equations. The errors decrease exponentially as $\mathcal{N}_2$ increases, demonstrating the spectral accuracy of the method. The logarithmic error plot further highlights the rapid convergence of the numerical solution.

6. Conclusions

This paper investigated a class of high-order P-VF-IDEs, which incorporate both Volterra and Fredholm integral components along with pantograph delay elements. We introduced a spectral tau approach for approximating solutions to P-VF-IDEs in one and two dimensions, utilizing operational differentiation and integration matrices to transform the continuous problem into a manageable system of algebraic equations. This method demonstrated high accuracy with a few numbers of computational modes. Through comprehensive numerical experiments, we highlighted the accuracy and convergence properties of the spectral Legendre tau method, affirming its effectiveness in solving high-order P-VF-IDEs in comparison to other spectral techniques. The results indicate that the proposed approach is a powerful tool for addressing intricate integro-differential equations with integral and proportional delay features.

Author contributions

M. A. Zaky: Writing–review & editing, Writing–original draft, Validation, Supervision, Software, Investigation; W. G. Alharbi: Validation, Formal analysis; M. M. Alzubaidi: Validation, Methodology, Writing the original draft. R. T. Matoog: Validation, Methodology, Writing the original draft. All authors have read and agreed to the published version of the manuscript.

Use of Generative-AI tools declaration

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

Acknowledgments

This work was supported and funded by the Deanship of Scientific Research at Imam Mohammad Ibn Saud Islamic University (IMSIU) (grant number IMSIU-DDRSP2503).

Funding

This work was supported and funded by the Deanship of Scientific Research at Imam Mohammad Ibn Saud Islamic University (IMSIU) (grant number IMSIU-DDRSP2503).

Conflict of interest

The authors declare that they do not have a conflict of interest.

References

[1]	M. Iqbal, A. Ganie, M. Miah, M. Osman, Extracting the ultimate new soliton solutions of some nonlinear time fractional PDEs via the conformable fractional derivative, Fractal Fract., 8 (2024), 210. https://doi.org/10.3390/fractalfract8040210 doi: 10.3390/fractalfract8040210
[2]	S. Ahmad, N. Becheikh, L. Kolsi, T. Muhammad, Z. Ahmad, M. Nasrat, Uncovering the stochastic dynamics of solitons of the Chaffee-Infante equation, Sci. Rep., 14 (2024), 19485. https://doi.org/10.1038/s41598-024-67116-4 doi: 10.1038/s41598-024-67116-4
[3]	I. Alraddadi, M. Chowdhury, M. Abbas, K. El-Rashidy, J. Borhan, M. Miah, et al., Dynamical behaviors and abundant new soliton solutions of two nonlinear PDEs via an efficient expansion method in industrial engineering, Mathematics, 12 (2024), 2053. https://doi.org/10.3390/math12132053 doi: 10.3390/math12132053
[4]	N. Mhadhbi, S. Gana, M. Alsaeedi, Exact solutions for nonlinear partial differential equations via a fusion of classical methods and innovative approaches, Sci. Rep., 14 (2024), 6443. https://doi.org/10.1038/s41598-024-57005-1 doi: 10.1038/s41598-024-57005-1
[5]	S. Zhang, G. Zhu, W. Huang, H. Wang, C. Yang, Y. Lin, Symbolic computation of analytical solutions for nonlinear partial differential equations based on bilinear neural network method, Nonlinear Dyn., in press. https://doi.org/10.1007/s11071-024-10715-7
[6]	Y. Ye, H. Fan, Y. Li, X. Liu, H. Zhang, Conformable bilinear neural network method: a novel method for time-fractional nonlinear partial differential equations in the sense of conformable derivative, Nonlinear Dyn., in press. https://doi.org/10.1007/s11071-024-10495-0
[7]	H. Hu, L. Qi, X. Chao, Physics-informed neural networks (PINN) for computational solid mechanics: numerical frameworks and applications, Thin-Wall. Struct., 205 (2024), 112495. https://doi.org/10.1016/j.tws.2024.112495 doi: 10.1016/j.tws.2024.112495
[8]	J. Linghu, W. Gao, H. Dong, Y. Nie, Higher-order multi-scale physics-informed neural network (HOMS-PINN) method and its convergence analysis for solving elastic problems of authentic composite materials, J. Comput. Appl. Math., 456 (2025), 116223. https://doi.org/10.1016/j.cam.2024.116223 doi: 10.1016/j.cam.2024.116223
[9]	S. Sarker, R. Karim, M. Akbar, M. Osman, P. Dey, Soliton solutions to a nonlinear wave equation via modern methods, J. Umm Al-Qura Univ. Appll. Sci., 10 (2024), 785–792. https://doi.org/10.1007/s43994-024-00137-x doi: 10.1007/s43994-024-00137-x
[10]	X. Gao, Symbolic computation on a (2+1)-dimensional generalized nonlinear evolution system in fluid dynamics, plasma physics, nonlinear optics and quantum mechanics, Qual. Theory Dyn. Syst., 23 (2024), 202. https://doi.org/10.1007/s12346-024-01045-5 doi: 10.1007/s12346-024-01045-5
[11]	X. Gao, In plasma physics and fluid dynamics: symbolic computation on a (2+1)-dimensional variable-coefficient Sawada-Kotera system, Appl. Math. Lett., 159 (2025), 109262. https://doi.org/10.1016/j.aml.2024.109262 doi: 10.1016/j.aml.2024.109262
[12]	P. Singh, K. Senthilnathan, Evolution of a solitary wave: optical soliton, soliton molecule and soliton crystal, Discov. Appl. Sci., 6 (2024), 464. https://doi.org/10.1007/s42452-024-06152-1 doi: 10.1007/s42452-024-06152-1
[13]	X, Gao, Hetero-Bäcklund transformation, bilinear forms and multi-solitons for a (2+1)-dimensional generalized modified dispersive water-wave system for the shallow water, Chinese J. Phys., 92 (2024), 1233–1239. https://doi.org/10.1016/j.cjph.2024.10.004 doi: 10.1016/j.cjph.2024.10.004
[14]	H. Wilhelmsson, Simultaneous diffusion and reaction processes in plasma dynamics, Phys. Rev. A, 38 (1988), 1482. https://doi.org/10.1103/PhysRevA.38.1482 doi: 10.1103/PhysRevA.38.1482
[15]	D. Anderson, R. Jancel, H. Wilhelmsson, Similarity solution of the evolution equation describing the combined effects of diffusion and recombination in plasmas, Phys. Rev. A, 30 (1984), 2113. https://doi.org/10.1103/PhysRevA.30.2113 doi: 10.1103/PhysRevA.30.2113
[16]	H. Wilhelmsson, Similarity solution of two coupled reaction-diffusion rate equations, Phys. Rev. A, 35 (1987), 1957. https://doi.org/10.1103/PhysRevA.35.1957 doi: 10.1103/PhysRevA.35.1957
[17]	R. Kumar, R. Kaushal, A. Prasad, Soliton-like solutions of certain types of nonlinear diffusion-reaction equations with variable coefficient, Phys. Lett. A, 372 (2008), 1862–1866. https://doi.org/10.1016/j.physleta.2007.10.061 doi: 10.1016/j.physleta.2007.10.061
[18]	R. Kumar, R. Kaushal, A. Prasad, Some new solitary and travelling wave solutions of certain nonlinear diffusion-reaction equations using auxiliary equation method, Phys. Lett. A, 372 (2008), 3395–3399. https://doi.org/10.1016/j.physleta.2008.01.062 doi: 10.1016/j.physleta.2008.01.062
[19]	M. Morgado, M. Rebelo, Numerical approximation of distributed order reaction-diffusion equations, J. Comput. Appl. Math., 275 (2015), 216–227. https://doi.org/10.1016/j.cam.2014.07.029 doi: 10.1016/j.cam.2014.07.029
[20]	J. Li, Barycentric rational collocation method for fractional reaction-diffusion equation, AIMS Mathematics, 8 (2023), 9009–9026. https://doi.org/10.3934/math.2023451 doi: 10.3934/math.2023451
[21]	M. Abdelrahman, W. Mohammed, M. Alesemi, S. Albosaily, The effect of multiplicative noise on the exact solutions of nonlinear Schrödinger equation, AIMS Mathematics, 6 (2021), 2970–2980. https://doi.org/10.3934/math.2021180 doi: 10.3934/math.2021180
[22]	S. Albosaily, W. Mohammed, M. Aiyashi, M. Abdelrahman, Exact solutions of the (2+1)-dimensional stochastic chiral nonlinear Schrödinger equation, Symmetry, 12 (2020), 1874. https://doi.org/10.3390/sym12111874 doi: 10.3390/sym12111874
[23]	W. Mohammed, H. Ahmad, H. Boulares, F. Khelifi, M. El-Morshedy, Exact solutions of Hirota-Maccari system forced by multiplicative noise in the Itô sense, J. Low Freq. Noise V. A., 41 (2022), 74–84. https://doi.org/10.1177/14613484211028100 doi: 10.1177/14613484211028100
[24]	W. Mohammed, H. Ahmad, A. Hamza, E. Aly, M. El-Morshedy, E. Elabbasy, The exact solutions of the stochastic Ginzburg-Landau equation, Results Phys., 23 (2021), 103988. https://doi.org/10.1016/j.rinp.2021.103988 doi: 10.1016/j.rinp.2021.103988
[25]	W. Mohammed, N. Iqbal, A. Ali, M. El-Morshedy, Exact solutions of the stochastic new coupled Konno-Oono equation, Results Phys., 21 (2021), 103830. https://doi.org/10.1016/j.rinp.2021.103830 doi: 10.1016/j.rinp.2021.103830
[26]	W. Mohammed, M. El-Morshedy, The influence of multiplicative noise on the stochastic exact solutions of the Nizhnik-Novikov-Veselov system, Math. Comput. Simulat., 190 (2021), 192–202. https://doi.org/10.1016/j.matcom.2021.05.022 doi: 10.1016/j.matcom.2021.05.022
[27]	W. Mohammed, S. Albosaily, N. Iqbal, M. El-Morshedy, The effect of multiplicative noise on the exact solutions of the stochastic Burgers' equation, Wave. Random Complex, 34 (2024), 274–286. https://doi.org/10.1080/17455030.2021.1905914 doi: 10.1080/17455030.2021.1905914
[28]	T. Shaikh, M. Baber, N. Ahmed, N. Shahid, A. Akgül, M. De la Sen, On the soliton solutions for the stochastic Konno-Oono system in magnetic field with the presence of noise, Mathematics, 11 (2023), 1472. https://doi.org/10.3390/math11061472 doi: 10.3390/math11061472
[29]	M. Baber, N. Ahmed, M. Yasin, S. Ali, M. Ali, A. Akgül, et al., Abundant soliton solution for the time-fractional stochastic Gray-Scot model under the influence of noise and M-truncated derivative, Discov. Appl. Sci., 6 (2024), 119. https://doi.org/10.1007/s42452-024-05759-8 doi: 10.1007/s42452-024-05759-8
[30]	E. Zayed, M. El-Horbaty, B. Saad, A. Arnous, Y. Yildirim, Novel solitary wave solutions for stochastic nonlinear reaction-diffusion equation with multiplicative noise, Nonlinear Dyn., 112 (2024), 20199–20213. https://doi.org/10.1007/s11071-024-10085-0 doi: 10.1007/s11071-024-10085-0
[31]	Sirendaoreji, Unified Riccati equation expansion method and its application to two new classes of Benjamin-Bona-Mahony equations, Nonlinear Dyn., 89 (2017), 333–344. https://doi.org/10.1007/s11071-017-3457-6 doi: 10.1007/s11071-017-3457-6
[32]	K. Hosseini, A. Bekir, R. Ansari, New exact solutions of the conformable time-fractional Cahn-Allen and Cahn-Hilliard equations using the modified Kudryashov method, Optik, 132 (2017), 203–209. https://doi.org/10.1016/j.ijleo.2016.12.032 doi: 10.1016/j.ijleo.2016.12.032
[33]	M. Alosaimi, M. Al-Malki, K. Gepreel, Optical soliton solutions in optical metamaterials with full nonlinearity, J. Nonlinear Opt. Phys., in press. https://doi.org/10.1142/S0218863524500218

This article has been cited by:

1.	Amin Ghoreyshi, Mostafa Abbaszadeh, Mahmoud A. Zaky, Mehdi Dehghan, Finite block method for nonlinear time-fractional partial integro-differential equations: Stability, convergence, and numerical analysis, 2025, 214, 01689274, 82, 10.1016/j.apnum.2025.03.002
2.	Mostafa Abbaszadeh, Mahmoud A. Zaky, Mehdi Dehghan, Virtual element approximation and BDF2 time-discrete scheme for a partial integro-differential equation with a singular Abel's kernel, 2025, 501, 00963003, 129451, 10.1016/j.amc.2025.129451
3.	Ali H. Tedjani, Sharifah E. Alhazmi, Samer S. Ezz-Eldien, An operational approach for one- and two-dimension high-order multi-pantograph Volterra integro-differential equation, 2025, 10, 2473-6988, 9274, 10.3934/math.2025426

Reader Comments

Your name:*

Email:*
© 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

AIMS Mathematics

1.8 3.4

Metrics

Article views(819) PDF downloads(176) Cited by(3)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(6)

AIMS Mathematics

Soliton solutions and stability analysis of the stochastic nonlinear reaction-diffusion equation with multiplicative white noise in soliton dynamics and optical physics

Related Papers:

Abstract

1. Introduction

2. Legendre polynomials

3. One-dimensional case

4. Two-dimensional case

5. Numerical results

5.1. Convergence test

5.2. The numerical stability test

5.3. Comparison test

5.4. Two-dimensional case

6. Conclusions

Author contributions

Use of Generative-AI tools declaration

Acknowledgments

Funding

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

AIMS Mathematics

Soliton solutions and stability analysis of the stochastic nonlinear reaction-diffusion equation with multiplicative white noise in soliton dynamics and optical physics

Related Papers:

Abstract

1. Introduction

2. Legendre polynomials

3. One-dimensional case

4. Two-dimensional case

5. Numerical results

5.1. Convergence test

5.2. The numerical stability test

5.3. Comparison test

5.4. Two-dimensional case

6. Conclusions

Author contributions

Use of Generative-AI tools declaration

Acknowledgments

Funding

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog