Testing the limits of cardiac electrophysiology models through systematic variation of current

1.   Introduction

High-order pantograph Volterra integro-differential equations provide a vital framework for describing a number of natural phenomena in physics and engineering ^[1,2,3]. With the integral form, the pantograph Volterra integro-differential equation has the ability to model systems with memory influence, and the pantograph term adds to the equation's capacity to study scaling properties. These attributes ensure the appearance of pantograph Volterra integro-differential equations in diverse fields, including wave propagation and quantum mechanics ^[4,5,6]. It has an essential value in studying systems where the scaled spatial or temporal variables and past states help in obtaining the behavior of the wave, for instance, quantum transitions, scattering processes, wave propagation in non-classical structures, and optical interference and diffraction ^[7,8].

On the other hand, the pantograph term as well as the differential and integral terms make it uneasy to find the exact solution of pantograph Volterra integro-differential equations, so many authors are interested in finding valid and accurate analytical and numerical solutions for different kinds of such equations. Abbaszadeh et al. ^[9] developed a virtual element framework for nonlinear partial integro-differential equations, proposing two temporal strategies (uniform and graded meshes) to address integral term singularity establishing a fully discrete scheme with rigorous error estimates, unconditional stability, convergence proofs, and validation through numerical experiments on varied physical domains. Zaky et al. ^[10] proposed a spectral tau approach to approximate the solution of high-order pantograph Volterra-Fredholm integro-differential equations in both one- and two-dimensions. Abdelhakem ^[11] utilized the Pseudo-spectral integration matrices for fractional Volterra integro-differential equations and Abel's integral equations. Ghoreyshi et al. ^[12] investigated a nonlinear time-fractional partial integro-differential equation by combining the weighted and shifted Grünwald–Letnikov formula for temporal discretization of the Caputo fractional derivative, the fractional trapezoidal rule for the Volterra integral, and Chebyshev spectral-collocation methods for spatial approximation. Ghosh and Mohapatra ^[13] applied an iterative technique for nonlinear delay Volterra integro-differential equations, where the integral term is approximated using the composite trapezoidal rule, and the Daftardar-Gejji and Jafari technique is used to solve the resulting implicit algebraic equation. Alsuyuti et al. ^[14] extended the application of the Galerkin spectral method for pantograph integro-differential and systems of pantograph differential equations of arbitrary order. Behera and Ray ^[15] carried out an operational scheme based on Bernoulli and müntz–Legendre wavelets (BMLW) for the solution of pantograph Volterra delay-integro-differential equations. Elkot et al. ^[16] presented the multi-variate Legendre-collocation spectral scheme for multi-dimensional nonlinear Volterra–Fredholm integral equations. Zhao et al. ^[17] utilized the Lagrange interpolation (LIT) and Bernstein tau (BT) methods, respectively, to deal with linear pantograph Volterra delay-integro-differential equations. The author in ^[18] applied the Jacobi spectral tau approach for systems of multi-pantograph equations. Bica and Satmari ^[19] employed an iterative approach based on the Bernstein composite quadrature and fuzzy Bernstein spline interpolation methods for the approximate solution of pantograph fuzzy Volterra integral equations. In ^[20], the spectral collocation technique is implemented, where the integral terms are approximated via the Gauss quadrature rule, for systems of multi-dimensional integral equations. In ^[21], the multistep collocation method is carried out for the pantograph second-kind Volterra integral equation. Other numerical methods include ^[22,23].

This manuscript aims to provide a numerical solution for the high-order multi-pantograph Volterra integro-differential equation

where $\xi(x), \ g(x), \ \mu_i(x), \ \nu_i(x)$ and $\gamma_i(x)\ (0 < i < s)$ are known functions and $\theta_i, \ c_i$ and $d_i\ (0 < i < s)$ are real numbers with $0 < c_i, d_i < 1$. In this regard, we use the operational matrix of differentiation and derive another one of pantograph, based on shifted Jacobi polynomials, that are utilized in conjunction with the spectral collocation approach to simplify the problem to an equivalent one of solving a system of algebraic equations. In addition, we study the application of the considered numerical technique to the system of high-order multi-pantograph Volterra integro-differential equations. On the other hand, we aim to extend the constructed numerical technique to solve the two-dimensional high-order multi-pantograph Volterra integro-differential equation

where $\mu_{i, j}(x, y), \ \xi_{i, j}(x, y), \ \nu(x, y), \ \theta_i(y)$ and $\vartheta_j(x)\ (0 < i < s, \ 0 < j < q)$ are known functions with $0 < c_i, d_j < 1.$

The paper is outlined as follows: Section 2 investigates the existence and uniqueness of the high-order multi-pantograph Volterra integro-differential equation (1.1). Section 3 provides some important properties of shifted Jacobi polynomials. Section 4 presents the application of operational matrices of pantograph and differentiation in conjunction with the Jacobi spectral collocation method to solve high-order multi-pantograph Volterra integro-differential equations. Section 5 studies the application of the considered numerical approach to the system of high-order multi-pantograph Volterra integro-differential equations. Section 6 discusses the extension of the numerical approach discussed in the previous section to solve the two-dimensional high-order multi-pantograph Volterra integro-differential equation. Section 7 introduces several test problems with their numerical solutions and comparisons with other spectral methods in the literature. Concluding remarks are displayed in Section 8.

2.   Existence and uniqueness

In this section we study the existence and uniqueness of the solution of the high-order multi-pantograph Volterra integro-differential equation (1.1).

Theorem 1. (^[24] Gronwall inequality) Let $f(x)$ be a non-negative integrable function over $(0, \alpha]$, and let $g(x)$ and $E(x)$ be continuous functions on $[0, \alpha]$. If $E(x)$ satisfies

then we have

If $g(x)$ is non-decreasing, then

Theorem 2. Assume $f(x)$ is a continuous function in $[0, \alpha];$ then the pantograph Volterra integro-differential equation (1.1) has a unique solution.

Proof. Firstly, we transform high-order multi-pantograph Volterra integro-differential equation (1.1) into its equivalent integral form as follows:

where

Suppose that for any continuous functions $f(x), \ h(x), \ \xi(x), \ \mu_i(x), \ \gamma_i(x)$ and $\nu_i(x)\ (0\leq i \leq s),$ in $[0, \alpha],$ and $c_i\ (0\leq i \leq s),$ with $0 < c_i, d_i < 1,$ the following inequalities hold:

Now, we prove that $\mathcal{K}\left(x, f(x)\right)$ satisfies the Lipschitz condition. For any two continuous functions $f(x), \ h(x),$ we have

where $\mathcal{L} = ABD+ AE\alpha+ C\alpha+AF$ is a positive constants. Consequently, Eq (2.4) can be rewritten as

Assume there are two solutions $f(x)$ and $h(x)$ of the pantograph Volterra integro-differential equation (1.1) satisfying $f^{(i)}(0) = h^{(i)}(0) = \theta_i, \ i = 0, 1, \cdots, s.$ This, along with Eq (2.6), yields

and then, it follows from Theorem 1 that $\left\| f(x)-h(x)\right\|_\infty = 0.$ Now, we consider

and define $\mathfrak{E}_{l+1} = f_{l+1}(x)-f_{l}(x).$ Following the proof in (^[25], Section 3), we get

Thus, for $n, m \rightarrow 0,$ we have

which implies that there is only one solution of (1.1). □

3.   Shifted Jacobi polynomials

Denote $S_{ \alpha, l}^{(a, b)}(x),$ for $l = 0, 1, \cdots,$ by the shifted Jacobi polynomials defined in $\Lambda_x = [0, \alpha]$, then it can be given analytically by

and satisfy

Shifted Jacobi polynomials satisfy the orthogonality relation

where $w_ \alpha^{(a, b)} = x^b (\alpha-x)^a$, $h_{ \alpha, l}^{(a, b)} = \dfrac{\alpha^{a+b+1} \Gamma(l+a+1) \Gamma(l+b+1)}{(2l+a+b+1) l! \Gamma(l+a+b+1)}$ and $\delta_{lm}$ is the well-known Kronecker delta function.

Any function $f\in L^2(\Lambda_x)$ can be expanded based on shifted Jacobi polynomials by:

The shifted Jacobi-Gauss quadrature rule satisfies

where $\left\{x^{(a, b)}_{\mathcal{L}, j}\right\}_{j = 0}^{\mathcal{L}}$ are the zeros of $S_{\mathcal{L}+1}^{(a, b)}(x)$, and $\left\{w^{(a, b)}_{\mathcal{L}, j}\right\}_{j = 0}^{\mathcal{L}}$ are their corresponding weights.

If we denote $P^{ \alpha}_{\mathcal{L}}$ by the orthogonal projection:

then, we have

where

Similarly, for $f\in L^2(\Lambda_x\times \Lambda_y), \ \Lambda_x = [0, \alpha], \ \Lambda_y = [0, \beta],$ we have

where $\mathcal{F}^T_{\mathcal{L}, \mathcal{M}}$ and $\mathfrak{S}^{ \alpha, \beta}_{\mathcal{L}, \mathcal{M}}(x, y)$ are defined by:

and

4.   Multi-pantograph Volterra integro-differential equations

In this section, we apply the Jacobi spectral collocation method to solve the following high-order multi-pantograph Volterra integro-differential equation (1.1).

The Jacobi spectral collocation approach for (1.1) is to find $f_{\mathcal{L}}\in \mathcal{S}_{\mathcal{L}}^ \alpha$, such that

Denote $f_{\mathcal{L}}(x)$ and $g_{\mathcal{L}}(x)$ by

with

The following two theorems will be of great use later.

Theorem 3. ^[26] The $i^{th}$ derivative of the vector $\mathfrak{S}^{ \alpha}_{\mathcal{L}}(x)$ is given by:

with

and

Theorem 4. For $0 \leq c \leq 1,$ the pantograph operational matrix $\boldsymbol{P}_{c}$ is given by

where

Proof. We start by expressing $S_{ \alpha, l}^{a, b}(cx)$ by:

Expanding $x^i$ in terms of $S_{ \alpha, j}^{a, b}(x), \ j = 0, 1, \cdots, \mathcal{L},$ by

A combination of (4.6) and (4.7) then yields

which completes the proof. □

Application of (4.2) with Theorems 4.1 and 4.2, we have

Substitution from (4.2), (4.3), (4.9), and (4.10) into (4.1); one gets

Thanks to (4.11), the residual $\mathcal{R}_{\mathcal{L}}(x)$ can be given by

Finally, the spectral solution of (1.1) is transformed to a problem of solving the following algebraic equations system:

5.   System of multi-pantograph Volterra integro-differential equations

Here, the operational approach discussed in the previous section is carried out to get a numerical solution for the system of high-order multi-pantograph Volterra integro-differential equations:

with

where $r, \ s, \ 0\leq c_i, d_i\leq 1\ (0\leq i \leq s)$ are real numbers and $g_l(x),$ $\gamma_{k, i}(x),$ $\mu_{k, i}(x),$ $\nu_{k, i}(x),$ $\xi_{l}(x)$ $(0\leq k \leq r, \ 1\leq i \leq s)$ are known functions defined in $0\leq x\leq \alpha$.

First, the system (5.1) may be written in the following way:

where

Now, we aim to find the vector $\textbf{F}(x)$ such that:

in this regard, we suppose

where $\textbf{f}_{l, j}\ (1\leq l\leq r, 0\leq j \leq \mathcal{L})$ are the unknowns need to be determined later, while $\textbf{g}_{l, j}\ (1\leq l\leq r, 0\leq j \leq \mathcal{L})$ can be given using the shifted Jacobi Gauss quadrature rule as follows:

Then one gets from (5.5) that

A combination of (5.5)–(5.7), Eq (5.4) can be approximated as:

Application of the spectral collocation method leads to the following system of $r (\mathcal{L}+1)$ algebraic equations:

6.   Two-dimensional case

We utilize the Jacobi spectral collocation technique to solve the two-dimensional high-order multi-pantograph Volterra integro-differential equation (1.2).

As a spectral collocation scheme, we have to find $f_{\mathcal{L}, \mathcal{M}}\in \mathcal{S}_{\mathcal{L}}^ \alpha\times\mathcal{S}_{\mathcal{M}}^ \beta$, such that

In this regard, we suppose

where $\mathcal{F}_{\mathcal{L}, \mathcal{M}}$ is an unknown vector, while

Definition 1. ^[27] The kronecker product of any two $l \times m$ and $n \times o$ dimensional matrices $P$ and $Q,$ respectively, is denoted by the $ln\times mo$ matrix $P \otimes Q$ and is given by:

Theorem 5. Assume $\mathcal{I}_{\mathcal{L}}$ and $\mathcal{I}_{\mathcal{M}}$ are the identity matrices of orders $\mathcal{L}+1$ and $\mathcal{M}+1,$ respectively; then

where $\mathcal{D}^{(i)}_x = \boldsymbol{D}^{(i)}\otimes \mathcal{I}_{\mathcal{L}}, \ \mathcal{D}^{(j)}_y = \mathcal{I}_{\mathcal{M}} \otimes \boldsymbol{D}^{(j)},$ $\mathcal{P}_{x, c} = \boldsymbol{P}_{c}\otimes \mathcal{I}_{\mathcal{M}},$ and $\mathcal{P}_{y, d} = \mathcal{I}_{\mathcal{L}}\otimes \boldsymbol{P}_{d}.$

In virtue of (6.2)–(6.6), we have

Then, the residual of (6.1) can be given by

Finally, we generate a system of $(\mathcal{L}+1)(\mathcal{M}+1)$ algebraic equations as follows:

7.   Test problems

The current section provides some test problems to ensure the validity of the numerical approaches presented in Sections 3–5. The programs used in this paper are performed using the PC machine, with CPU Intel(R) Core(TM) i3-2350M 2 Duo CPU 2.30 GHz and 6.00 GB of RAM. We also used the arithmetic symbolic program known as (Mathematica 12) to perform intermediate arithmetic operations and arithmetic tables, as well as illustrations in the paper as a whole.

7.1.   Stability test

First, we consider the pantograph Volterra integro-differential equation ^[28,29]:

with an exact solution $f(x) = \sin(x).$

To evaluate the stability of our proposed scheme, we introduce controlled perturbations into both the source term and the initial condition of Eq (7.1). This allows us to analyze how sensitive the method is to variations in input parameters. Consider the perturbed version of (7.1):

where $\epsilon_s$ and $\epsilon_i$ are perturbations to the source term and the initial conditions, respectively. The exact solution remains $g(x) = h(x) = \sin(x)$.

The numerical stability of the proposed scheme is analyzed by introducing perturbations to the source term ($\epsilon_s$) in (7.2) and initial conditions ($\epsilon_i$) in (7.3). Tables 1 and 2 show that both perturbations exhibit linear error propagation:

Table 1.  MAEs for the perturbed Problem (7.2) with source term perturbations.

$\mathcal{L}$	$\epsilon_s=0.1$	$\epsilon_s=0.01$	$\epsilon_s=0.001$
6	0.05255861	0.00525315	0.00052261
8	0.05256080	0.00525608	0.00052561
10	0.05256079	0.00525607	0.00052560

---

 | Show Table

DownLoad: CSV

Table 2.  MAEs for the perturbed Problem (7.3) with initial condition perturbations.

$\mathcal{L}$	$\epsilon_i=0.1$	$\epsilon_i=0.01$	$\epsilon_i =0.001$
6	0.19715309	0.01971260	0.00196855
8	0.19715418	0.01971542	0.00197154
10	0.19715416	0.01971541	0.00197154

---

 | Show Table

DownLoad: CSV

● For $\epsilon_s$, reducing the perturbation magnitude by a factor of $10$ (e.g., $0.1 \to 0.01 \to 0.001$) decreases the maximum absolute errors (MAEs), $\lvert f_{\mathcal{L}}(x) - g_{\mathcal{L}}(x) \rvert$, by the same factor.

● Similarly, for $\epsilon_i$, the MAEs, $\lvert f_{\mathcal{L}}(x) - h_{\mathcal{L}}(x) \rvert$, scale linearly with the perturbation size.

This linear behavior indicates the method's robustness against small perturbations in the source term or the initial condition.

Previous studies by Yüzbasi ^[28] (Laguerre operational matrix, LOM) and Gülsu & Sezer ^[29] (Taylor collocation method, TCM) provided approximate solutions. Table 3 compares the numerical results of the absolute error function of $f(x)$ using the numerical method presented in Section 3 against those given using the LOM ^[28] and TCM ^[29]. Figure 1 plots the absolute errors of $f(x)$ with $\mathcal{L} = 10$ and the logarithmic function of MAEs for $f(x)$ at $(a, b) = (0, 0)$.

Table 3.  Absolute errors of $f(x)$ with $(a, b) = (1, 0)$ for Problem (7.1).

	LOM [28]			TCM [29]			Our scheme
$x$	$N=7$	$N=10$		$N=7$	$N=10$		$\mathcal{L}=7$	$\mathcal{L}=10$
0.0	$0.00$	$0.00$		$0.00$	$0.00$		$5.55\times10^{-17}$	$5.55\times10^{-17}$
0.2	$2.34\times10^{-9}$	$3.04\times10^{-13}$		$2.35\times10^{-7}$	$7.00\times10^{-10}$		$4.05\times10^{-10}$	$9.99\times10^{-15}$
0.4	$5.13\times10^{-9}$	$6.63\times10^{-13}$		$5.27\times10^{-7}$	$1.50\times10^{-9}$		$6.29\times10^{-10}$	$1.90\times10^{-14}$
0.6	$8.33\times10^{-9}$	$1.03\times10^{-12}$		$9.32\times10^{-7}$	$2.40\times10^{-9}$		$6.07\times10^{-10}$	$7.66\times10^{-15}$
0.8	$1.15\times10^{-8}$	$1.43\times10^{-12}$		$5.60\times10^{-7}$	$3.00\times10^{-9}$		$5.86\times10^{-10}$	$2.14\times10^{-14}$
1.0	$3.64\times10^{-9}$	$3.48\times10^{-12}$		$3.37\times10^{-5}$	$1.06\times10^{-7}$		$9.27\times10^{-12}$	$3.33\times10^{-16}$

---

 | Show Table

DownLoad: CSV

Figure 1.  The logarithmic function of MAEs and absolute errors of $f(x)$ with $(a, b) = (0, 0)$ and $\mathcal{L} = 10$ for Problem (7.1).

DownLoad: Full-Size Img PowerPoint

7.2.   Convergence test

To rigorously evaluate the convergence properties of the proposed method, we analyze the Volterra integro-differential equation:

with $f(0) = 0$ and exact solution $f(x) = \sin(x)e^{2x-1}$. This benchmark problem, previously studied in ^[15,17], allows direct comparison of convergence rates.

Table 4 demonstrates the exponential convergence of our method compared to those achieved using the BMLW presented in ^[15], while Table 5 compares the absolute errors of $f(x)$ with those given using the BT and LIT methods introduced in ^[17]. Figure 2 plots the absolute errors of $f(x)$ with $\mathcal{L} = 16$ and the logarithmic function of MAEs for $f(x)$ at $(a, b) = (1, 1)$.

Table 4.  MAEs of $f(x)$ for Problem (7.4).

	$l=3$	$l=4$	$l=5$	$l=6$
BMLW ($\mathcal{M}=3$) [15]	$6.563\times10^{-4}$	$8.508\times10^{-5}$	$1.080\times10^{-5}$	$1.359\times10^{-6}$
	$\mathcal{L}=6$	$\mathcal{L}=8$	$\mathcal{L}=10$	$\mathcal{L}=12$
Our method $(a, b)=(2, 2)$	$4.507\times10^{-4}$	$4.831\times10^{-6}$	$1.137\times10^{-8}$	$9.835\times10^{-12}$

---

 | Show Table

DownLoad: CSV

Table 5.  Absolute errors of $f(x)$ at $x = 0.5$ and 1 with $(a, b) = (2, 0)$ for Problem (7.4).

	BT [17]			LIT [17]			Our method
$\mathcal{L}$	$x=0.5$	$x=1$		$x=0.5$	$x=1$		$x=0.5$	$x=1$
6	$2.75\times10^{-2}$	$3.17\times10^{-1}$		$1.85\times10^{-5}$	$4.20\times10^{-3}$		$2.47\times10^{-6}$	$6.79\times10^{-5}$
8	$1.99\times10^{-2}$	$2.46\times10^{-1}$		$2.00\times10^{-7}$	$2.30\times10^{4}$		$2.08\times10^{-8}$	$7.09\times10^{-7}$
10	$1.56\times10^{-2}$	$2.06\times10^{-1}$		$2.37\times10^{-9}$	$6.68\times10^{-6}$		$3.77\times10^{-11}$	$1.64\times10^{-9}$
12	$1.29\times10^{-2}$	$1.80\times10^{-1}$		$5.66\times10^{-11}$	$7.75\times10^{-7}$		$2.47\times10^{-14}$	$1.82\times10^{-12}$

---

 | Show Table

DownLoad: CSV

Figure 2.  The logarithmic function of MAEs and absolute errors of $f(x)$ with $(a, b) = (1, 1)$ and $\mathcal{L} = 16$ for Problem (7.4).

DownLoad: Full-Size Img PowerPoint

7.3.   Irregular solution test

To test the validity of the numerical approach presented in Section 3 for pantograph Volterra integro-differential equations with irregular solutions, we consider the following problem:

where $f(0) = \frac{df}{dx}(0) = 0$ and $g(x)$ is chosen such that the exact solution is $f(x) = x^{\frac{7}{2}}$.

To solve this problem, we apply the numerical scheme introduced in Section 3 with various choices of $\mathcal{L}$. Table 6 displays the MAEs of $f_{\mathcal{L}}(x)$ at $(a, b) = (0, 0)$ with $\mathcal{L} = 4, 8, 12, 16,$ and $20$.

Table 6.  MAEs of $f_{\mathcal{L}}(x)$ for Problem (7.5).

$\mathcal{L}$	4	8	12	16	20
MAE	$1.1425\times10^{-3}$	$7.3101\times10^{-6}$	$1.0095\times10^{-6}$	$2.7518\times10^{-7}$	$9.8240\times10^{-8}$

---

 | Show Table

DownLoad: CSV

7.4.   System of equations

Consider the system of pantograph equations:

The exact solutions are given by $f_1(x) = e^{x}$ and $f_2(x) = x^3$.

In ^[30], shifted Chebyshev polynomials were used as basis functions for the collocation spectral approach (CCS) to solve (7.6) with $c = 0.8$. Table 7 compares the MAEs of our method (using the numerical approach from Section 4 with $(a, b) = (0, 0)$ against the CCS results. Figure 3 shows the absolute error distributions for both components with $\mathcal{L} = 12$, along with the logarithmic MAE plots for $(a, b) = (1, 0)$ with $c = 0.1$.

Table 7.  MAEs of $f_1(x)$ and $f_2(x)$ for Problem (7.6).

		CCS [30]			Our scheme
$\mathcal{L}$		$f_1(x)$	$f_2(x)$		$f_1(x)$	$f_2(x)$
2		$6.0777\times10^{-2}$	$1.1843\times10^{-1}$		$6.5075\times10^{-2}$	$4.8397\times10^{-2}$
4		$9.2227\times10^{-4}$	$8.2963\times10^{-4}$		$3.9735\times10^{-5}$	$6.6931\times10^{-5}$
6		$7.3973\times10^{-7}$	$1.8431\times10^{-6}$		$5.4490\times10^{-8}$	$2.8398\times10^{-9}$
8		$9.5045\times10^{-10}$	$6.6164\times10^{-10}$		$4.7068\times10^{-11}$	$1.9821\times10^{-12}$

---

 | Show Table

DownLoad: CSV

Figure 3.  The logarithmic functions of MAEs and absolute errors of $f_1(x)$ and $f_2(x)$ with $(a, b) = (1, 0)$ and $\mathcal{L} = 12$ for Problem (7.6).

DownLoad: Full-Size Img PowerPoint

7.5.   Two-dimensional case

Consider the two-dimensional pantograph equation:

defined on the domain $0 \leq x, y \leq 1$, with exact solution $f(x, y) = (x+1)\sin(y)$.

To validate the numerical scheme from Section 5, Table 8 presents the absolute errors of $f(x, y)$ for different values of $\mathcal{L}$ and $\mathcal{M}$ with $(a, b) = (2, 2)$. Figure 4 shows the error distribution for $\mathcal{L} = \mathcal{M} = 10$ at $(a, b) = (3, 3)$, while Figure 5 displays corresponding contour plots for various values of $\mathcal{L}$ and $\mathcal{M}$.

Table 8.  Absolute errors of $f(x, y)$ with $(a, b) = (2, 2)$ for Problem (7.7).

$(x_1, x_2)$	$(\mathcal{L}, \mathcal{M})=(2, 2)$	$(\mathcal{L}, \mathcal{M})=(4, 4)$	$(\mathcal{L}, \mathcal{M})=(6, 6)$	$(\mathcal{L}, \mathcal{M})=(8, 8)$	$(\mathcal{L}, \mathcal{M})=(10, 10)$
(0.1, 0.1)	$4.7780\times10^{-3 }$	$7.5011\times10^{-6 }$	$8.7689\times10^{-8 }$	$1.8082\times10^{-10}$	$1.8707\times10^{-13}$
(0.2, 0.2)	$1.0733\times10^{-3 }$	$4.6338\times10^{-5 }$	$1.7930\times10^{-7 }$	$3.0610\times10^{-10}$	$2.7519\times10^{-13}$
(0.3, 0.3)	$4.8745\times10^{-3 }$	$6.9301\times10^{-5 }$	$2.4872\times10^{-7 }$	$3.6300\times10^{-10}$	$2.8521\times10^{-13}$
(0.4, 0.4)	$7.1270\times10^{-3 }$	$9.2293\times10^{-5 }$	$2.8985\times10^{-7 }$	$3.7009\times10^{-10}$	$3.0797\times10^{-13}$
(0.5, 0.5)	$8.6790\times10^{-3 }$	$1.1553\times10^{-4 }$	$2.9407\times10^{-7 }$	$3.9560\times10^{-10}$	$3.2218\times10^{-13}$
(0.6, 0.6)	$1.0695\times10^{-2 }$	$1.2945\times10^{-4 }$	$2.9661\times10^{-7 }$	$4.1656\times10^{-10}$	$3.3228\times10^{-13}$
(0.7, 0.7)	$1.4623\times10^{-2 }$	$1.2567\times10^{-4 }$	$3.2109\times10^{-7 }$	$4.1231\times10^{-10}$	$3.4572\times10^{-13}$
(0.8, 0.8)	$2.2154\times10^{-2 }$	$1.1240\times10^{-4 }$	$3.2223\times10^{-7 }$	$4.3780\times10^{-10}$	$3.4061\times10^{-13}$
(0.9, 0.9)	$3.5178\times10^{-2 }$	$1.3403\times10^{-4 }$	$2.8893\times10^{-7 }$	$4.0347\times10^{-10}$	$3.5038\times10^{-13}$

---

 | Show Table

DownLoad: CSV

Figure 4.  Absolute error function of $f(x, y)$ with $(a, b) = (3, 3)$ and $\mathcal{L} = \mathcal{M} = 10$ for Problem (7.7).

DownLoad: Full-Size Img PowerPoint

Figure 5.  The contour plot of absolute errors of $f(x, y)$ with $(a, b) = (3, 3)$ for Problem (7.7).

DownLoad: Full-Size Img PowerPoint

8.   Conclusions

This paper dealt with a problem of great importance in physics and engineering, namely, high-order multi-pantograph Volterra integro-differential equations with variable coefficients. The pantograph and integral terms enable the equation to study complex systems with memory effects and scaling. We investigated the existence and uniqueness, for the first time, of the solution to the high-order multi-pantograph Volterra integro-differential equation. We introduced the pantograph operational matrix, for the first time based on shifted Jacobi polynomials, and employed it with the operational matrix of differentiation to convert the studied problem into a system of algebraic equations with the aid of the spectral collocation method. The studied numerical approach is also applied for the system of high-order multi-pantograph Volterra integro-differential equations and for high-order two-dimensional multi-pantograph Volterra integro-differential equations with variable coefficients. Comparing the numerical results of five test problems with the exact solution and with other spectral methods applied to the same problems confirms the superiority of the new one. The studied numerical approach is shown to be more accurate than the Laguerre operational matrix, Taylor collocation, Chebyshev collocation spectral, Bernoulli müntz–Legendre wavelet, Lagrange interpolation, and Bernstein tau methods. At the end, we note that the presented approach can be extended to more complex integro-differential equations. This will be the subject of our future research.

Author contributions

A.H. Tedjani: Validation, Methodology, Writing the original draft; S.E. Alhazmi: Validation, Methodology, Writing the original draft; S.S. Ezz-Eldien: Writing–review & editing, Writing–original draft, Validation, Supervision, Software, Investigation. All authors have read and agreed to the published version of the manuscript.

Use of Generative-AI tools declaration

The authors declare 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-DDRSP2502).

Conflicts of interest

The authors declare that they have no conflicts of interest.

Funding

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