A Legendre tau approach for high-order pantograph Volterra-Fredholm integro-differential equations

Mahmoud A. Zaky; Weam G. Alharbi; Marwa M. Alzubaidi; R.T. Matoog; Mahmoud A. Zaky; Weam G. Alharbi; Marwa M. Alzubaidi; R.T. Matoog

doi:10.3934/math.2025322

AIMS Mathematics

2025, Volume 10, Issue 3: 7067-7085. doi: 10.3934/math.2025322

Previous Article Next Article

Research article Special Issues

A Legendre tau approach for high-order pantograph Volterra-Fredholm integro-differential equations

1.
Department of Mathematics and Statistics, College of Science, Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, Saudi Arabia
2.
Department of Mathematics, Faculty of Science, University of Tabuk, Tabuk 71491, Saudi Arabia
3.
Department of Mathematics, College of Duba, University of Tabuk, Tabuk, Saudi Arabia
4.
Mathematics Department, Faculty of Sciences, Umm Al-Qura University, Makkah, Saudi Arabia

Received: 16 January 2025 Revised: 17 March 2025 Accepted: 21 March 2025 Published: 27 March 2025
MSC : 65L05, 65R20, 65N35, 65L03

High-order pantograph Volterra-Fredholm integro-differential equations (P-VF-IDEs) arise in various scientific and engineering applications, where systems exhibit delay or scaling in their dynamics. This paper investigates a class of high-order P-VF-IDEs characterized by the presence of both Volterra and Fredholm integral terms as well as pantograph delay elements. We propose a spectral tau approach to approximate the solution P-VF-IDEs in both one- and two-dimensions. In particular, we employ the operational differentiation and integration matrices to approximate the integro-differential operator, transforming the continuous problem into a system of algebraic equations, and providing high accuracy with fewer computational modes. Numerical experiments illustrate the accuracy and convergence properties of the spectral Legendre tau method in solving high-order P-VF-IDEs and demonstrate its efficacy compared to other spectral approaches.

Keywords:

pantograph Volterra-Fredholm,
High-order integro-differential equations,
spectral method,
Legendre tau

Citation: Mahmoud A. Zaky, Weam G. Alharbi, Marwa M. Alzubaidi, R.T. Matoog. A Legendre tau approach for high-order pantograph Volterra-Fredholm integro-differential equations[J]. AIMS Mathematics, 2025, 10(3): 7067-7085. doi: 10.3934/math.2025322

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]	J. R. Ockendon, A. B. Tayler, The dynamics of a current collection system for an electric locomotive, Proc. Royal Soc. London A. Math. Phys. Sci., 322 (1971), 447–468. https://doi.org/10.1098/rspa.1971.0078 doi: 10.1098/rspa.1971.0078
[2]	M. Buhmann, A. Iserles, Stability of the discretized pantograph differential equation, Math. Comput., 60 (1993), 575–589. https://doi.org/10.1090/S0025-5718-1993-1176707-2 doi: 10.1090/S0025-5718-1993-1176707-2
[3]	P. Rahimkhani, Y. Ordokhani, E. Babolian, Müntz-legendre wavelet operational matrix of fractional-order integration and its applications for solving the fractional pantograph differential equations, Numer. Algor., 77 (2018), 1283–1305. https://doi.org/10.1007/s11075-017-0363-4 doi: 10.1007/s11075-017-0363-4
[4]	S. Behera, S. S. Ray, An efficient numerical method based on euler wavelets for solving fractional order pantograph volterra delay-integro-differential equations, J. Comput. Appl. Math., 406 (2022), 113825. https://doi.org/10.1016/j.cam.2021.113825 doi: 10.1016/j.cam.2021.113825
[5]	S. S. Ezz-Eldien, A. Alalyani, Legendre spectral collocation method for one- and two-dimensional nonlinear pantograph Volterra-Fredholm integro-differential equations, Int. J. Mod. Phys. C, 2025. https://doi.org/10.1142/S0129183125500615
[6]	N. A. Elkot, M. A. Zaky, E. H. Doha, I. G. Ameen, On the rate of convergence of the legendre spectral collocation method for multi-dimensional nonlinear volterra–fredholm integral equations, Commun. Theory Phys., 73 (2021), 025002. https://doi.org/10.1088/1572-9494/abcfb3 doi: 10.1088/1572-9494/abcfb3
[7]	D. Trif, Direct operatorial tau method for pantograph-type equations, Appl. Math. Comput., 219 (2012), 2194–2203. https://doi.org/10.1016/j.amc.2012.08.065 doi: 10.1016/j.amc.2012.08.065
[8]	S. S. Ezz-Eldien, On solving systems of multi-pantograph equations via spectral tau method, Appl. Math. Comput., 321 (2018), 63–73. https://doi.org/10.1016/j.amc.2017.10.014 doi: 10.1016/j.amc.2017.10.014
[9]	M. Sezer, S. Yalçinbaş, M. Gülsu, A taylor polynomial approach for solving generalized pantograph equations with nonhomogenous term, Int. J. Comput. Math., 85 (2008), 1055–1063. https://doi.org/10.1080/00207160701466784 doi: 10.1080/00207160701466784
[10]	F. Shakeri, M. Dehghan, Application of the decomposition method of adomian for solving the pantograph equation of order m, Zeitschrift für Naturforschung A, 65 (2010), 453–460. https://doi.org/10.1515/zna-2010-0510 doi: 10.1515/zna-2010-0510
[11]	Ş. Yüzbaşı, E. Gök, M. Sezer, Laguerre matrix method with the residual error estimation for solutions of a class of delay differential equations, Math. Meth. Appl. Sci., 37 (2014), 453–463. https://doi.org/10.1002/mma.2801 doi: 10.1002/mma.2801
[12]	S. S. Ezz-Eldien, E. H. Doha, Fast and precise spectral method for solving pantograph type volterra integro-differential equations, Numer. Algor., 21 (2019), 57–77. https://doi.org/10.1007/s11075-018-0535-x doi: 10.1007/s11075-018-0535-x
[13]	A. Ghoreyshi, M. Abbaszadeh, M. A. Zaky, M. Dehghan, Finite block method for nonlinear time-fractional partial integro-differential equations: stability, convergence, and numerical analysis, Appl. Numer. Math., 2025. https://doi.org/10.1016/j.apnum.2025.03.002
[14]	C. Yang, J. Hou, Jacobi spectral approximation for boundary value problems of nonlinear fractional pantograph differential equations, Numer. Algor., 86 (2021), 1089–1108. https://doi.org/10.1007/s11075-020-00924-7 doi: 10.1007/s11075-020-00924-7
[15]	F. Ghomanjani, M. H. Farahi, N. Pariz, A new approach for numerical solution of a linear system with distributed delays, volterra delay-integro-differential equations, and nonlinear volterra-fredholm integral equation by bezier curves, Comput. Appl. Math., 36 (2017), 1349–1365. https://doi.org/10.1007/s40314-015-0296-2 doi: 10.1007/s40314-015-0296-2
[16]	Q. Huang, H. Xie, H. Brunner, Superconvergence of discontinuous galerkin solutions for delay differential equations of pantograph type, SIAM J. Sci. Comput., 33 (2011), 2664–2684. https://doi.org/10.1137/110824632 doi: 10.1137/110824632
[17]	Y. Qin, C. Huang, An hp-version error estimate of spectral collocation methods for weakly singular volterra integro-differential equations with vanishing delays, Comput. Appl. Math., 43 (2024), 301. https://doi.org/10.1007/s40314-024-02818-z doi: 10.1007/s40314-024-02818-z
[18]	M. A. Zaky, An accurate spectral collocation method for nonlinear systems of fractional differential equations and related integral equations with nonsmooth solutions, Appl. Numer. Math., 154 (2020), 205–222. https://doi.org/10.1016/j.apnum.2020.04.002 doi: 10.1016/j.apnum.2020.04.002
[19]	M. A. Zaky, I. G. Ameen, A novel jacobi spectral method for multi-dimensional weakly singular nonlinear volterra integral equations with nonsmooth solutions, Eng. Comput., 37 (2021), 2623–2631. https://doi.org/10.1007/s00366-020-00953-9 doi: 10.1007/s00366-020-00953-9
[20]	M. A. Zaky, I. G. Ameen, A priori error estimates of a Jacobi spectral method for nonlinear systems of fractional boundary value problems and related volterra-fredholm integral equations with smooth solutions, Numer. Algor., 84 (2020), 63–89. https://doi.org/10.1007/s11075-019-00743-5 doi: 10.1007/s11075-019-00743-5
[21]	H. Moussa, M. Saker, M. A. Zaky, M. A. Babatin, S. S. Ezz-Eldien, Mapped Legendre-spectral method for high-dimensional multi-term time-fractional diffusion-wave equation with non-smooth solution, Comput. Appl. Math., 44 (2025), 167. https://doi.org/10.1007/s40314-025-03123-z doi: 10.1007/s40314-025-03123-z
[22]	P. Rahimkhani, Y. Ordokhani, Hahn wavelets collocation method combined with Laplace transform method for solving fractional integro-differential equations, Math. Sci., 18 (2024), 463–477. http://doi.org/10.1007/s40096-023-00514-3 doi: 10.1007/s40096-023-00514-3
[23]	E. Doha, A. Bhrawy, S. Ezz-Eldien, An efficient Legendre spectral tau matrix formulation for solving fractional subdiffusion and reaction subdiffusion equations, J. Comput. Nonlin. Dyn., 10 (2015), 021019. http://doi.org/10.1115/1.4027944 doi: 10.1115/1.4027944
[24]	S. S. Ezz-Eldien, R. M. Hafez, A. H. Bhrawy, D. Baleanu, A. A. El-Kalaawy, New numerical approach for fractional variational problems using shifted legendre orthonormal polynomials, J. Opt. Theory Appl., 174 (2017), 295–320. http://doi.org/10.1007/s10957-016-0886-1 doi: 10.1007/s10957-016-0886-1
[25]	J. Zhao, Y. Cao, Y. Xu, Sinc numerical solution for pantograph volterra delay-integro-differential equation, Int. J. Comput. Math., 94 (2017), 853–865. http://doi.10.1080/00207160.2016.1149577 doi: 10.1080/00207160.2016.1149577

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(341) PDF downloads(21) Cited by(3)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

AIMS Mathematics

A Legendre tau approach for high-order pantograph Volterra-Fredholm integro-differential equations

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

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