A mathematical model for assessing the impact of dual-level toxicity on aquatic biospecies and its optimal control analysis

K. O. Achema; D. Okuonghae; C. J. Alhassan; K. O. Achema; D. Okuonghae; C. J. Alhassan

doi:10.3934/mmc.2022012

Mathematical Modelling and Control

2022, Volume 2, Issue 3: 100-121. doi: 10.3934/mmc.2022012

Previous Article Next Article

Research article Special Issues

A mathematical model for assessing the impact of dual-level toxicity on aquatic biospecies and its optimal control analysis

1.
Department of Mathematics and Statistics, Confluence University of Science and Technology, Osara, Kogi State, Nigeria
2.
Department of Mathematics, University of Benin, Benin City, Nigeria
3.
Department of Mathematics and Computer Science, Edo State University, Uzairue, Edo State, Nigeria

Received: 28 January 2022 Revised: 26 May 2022 Accepted: 19 July 2022 Published: 27 September 2022

Ecological models have become paramount for assessing the pesticides effect on the function and structure of aquatic ecosystems. The most paramount concerns are assessments of pesticides/toxicants that have the potential to change from one form to another when they are released into the aquatic ecosystem. Optimal control model is formulated from the nonlinear mathematical model for assessing dual-level toxicity of pesticides effect on aquatic species with the goal to minimizing the pesticides concentration in the aquatic species environment and maximizing the aquatic species population. Two control functions were introduced to represent a policy of not allowing pesticides concentration into the aquatic species environment and the removal of those pesticides that are already in the aquatic environment. The resulting optimal controls are characterized in terms of the optimality system and it was solved quantitatively for different scenarios using both forward and backward sweep iterative method with Runge-Kutta fourth order scheme. The result of the system showed different levels of the aquatic species population stability due to the different levels of the pesticides influx. It was also observed that the degradation of pesticides concentration causes pesticides concentration to vary significantly between the water body and the sediment region with significant level effect on the aquatic species.

Keywords:

Citation: K. O. Achema, D. Okuonghae, C. J. Alhassan. A mathematical model for assessing the impact of dual-level toxicity on aquatic biospecies and its optimal control analysis[J]. Mathematical Modelling and Control, 2022, 2(3): 100-121. doi: 10.3934/mmc.2022012

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]	V. Aihie, D. Okuonghae, Optimal control measures for tuberculosis mathematical model including immigration and isolation of infective, J. Biol. Syst., 18 (2010), 17–54. https://doi.org/10.1142/S0218339010003160 doi: 10.1142/S0218339010003160
[2]	T. C. M. Brock, R. P. A. Wijngarden, P. J. Van Den Brink, Threshold levels for effects of insecticides in freshwater ecoystems: a review, Ecotoxicology, 14 (2005), 355–380. https://doi.org/10.1007/s10646-004-6371-x doi: 10.1007/s10646-004-6371-x
[3]	L. Feiyu, X. Mengshi, L. Derong, L. Xraomei, W. Zhijun, H. Yichen, et al., Microbial degradation of pesticide residue and an emphasis on the degradation of cypermethrin and 3-phenoxy benzoic acid: a review, Molecules, 11 (2018), 2313. https://doi.org/10.3390/molecules23092313 doi: 10.3390/molecules23092313
[4]	S. A. Patin, Pollution and the Biological Resources of the Ocean, Butter worth scientific, London, 1982.
[5]	A. Rescigno, The struggle for life-v. one species living in a limited environment, B. Math. Biol., 39 (1997), 479–485. https://doi.org/10.1016/S0092-8240(77)90008-8 doi: 10.1016/S0092-8240(77)90008-8
[6]	A. Hasan, A. K. Misra, A. Kumar, A. K. Agrawal, Modelling the effect of toxicant on the deformity in a subclass of a biological species, Model. Earth Syst. Env., 2 (2016), 40. https://doi.org/10.1007/s40808-016-0086-x doi: 10.1007/s40808-016-0086-x
[7]	B. Ghosh, T. Kar, Sustainability and optimal control of exploited prey predator system through provision of alternative food predator, Biosystems, 109 (2012), 220–232. https://doi.org/10.1016/j.biosystems.2012.02.003 doi: 10.1016/j.biosystems.2012.02.003
[8]	N. Othax, J. G. Castelain, F. Peluso, S. Dubny, Environmental risk of pesticides applying the del azul pest risk model to freshwater of an agricultural area of Argentina, Human and Ecological Risk Assessment: An International Journal, 20 (2014), 1177–1199. https://doi.org/10.1080/10807039.2014.883800 doi: 10.1080/10807039.2014.883800
[9]	W. Liu, C. Xu, W. Tu, L. Niu, Embryonic exposure to butachlor in zebra fish (danio renio) endocrine disruption, developmental toxicity and immunotoxicity, Ecology and Environmental Safety, 89 (2013), 189–195. https://doi.org/10.1016/j.ecoenv.2012.11.031 doi: 10.1016/j.ecoenv.2012.11.031
[10]	M. Chowdhury, M. Rahman, M. A. Uddin, M. Saha, Pesticide residues in some selected pond water samples of meherpur region of Bangladesh, Journal of the Asiative Society of Bangladesh, Science, 39 (2013), 77–82. https://doi.org/10.3329/jasbs.v39i1.16036 doi: 10.3329/jasbs.v39i1.16036
[11]	United State Geological Survey (USGS) Fact Sheet, Pesticides in stream sediment and aquatic biota, 2019. Available from: http://water.usgs.gov/nawqa.
[12]	H. D. Miller, P. D. Noyees, M. K. McElwee, The toxicology of climate change: environmental contaminants in a worming world, Environmental International, 35 (2009), 971–986. https://doi.org/10.1016/j.envint.2009.02.006 doi: 10.1016/j.envint.2009.02.006
[13]	A. K. Agrawal, J. B. Shukla, Some mathematical models in ecotoxicology: Effects of toxicants on biological species, Sadahana, 24 (1999), 25–40. https://doi.org/10.1007/BF02747550 doi: 10.1007/BF02747550
[14]	S. M. Chowdhury, J. T. Chowdhury, S. F. Ahmed, P. Agarwal, I. A. Badruddin, S. Kamangar, Mathematical modelling of COVID-19 disease dynamics: Interaction between immune system and SARS-CoV-2 within host, AIMS Mathematics, 7 (2021), 2618–2633. https://doi.org/10.3934/math.2022147 doi: 10.3934/math.2022147
[15]	E. Lawrence, O. Ozekeke, I. Tongo, Risk assessment of agricultural pesticides in water, sediments, and fish from Owan River, Edo State, Nigeria, Environ. Monit. Assess., 187 (2015), 654. https://doi.org/10.1007/s10661-015-4840-8 doi: 10.1007/s10661-015-4840-8
[16]	R. B. Schäfer, P. J. van den Brink, M. Liess, Impacts of pesticides on freshwater ecosystems, Ecological impacts of toxic chemicals, 2011 (2011), 111–137. https://doi.org/10.2174/978160805121211101010111 doi: 10.2174/978160805121211101010111
[17]	C. Baines, A. Lerebours, F. Thomas, J. Fort, R. Kreitsberg, S. Gentes, et al., Linking pollution and cancer in aquatic environments: A review, Environ. Int., 149 (2021), 106391. https://doi.org/10.1016/j.envint.2021.106391 doi: 10.1016/j.envint.2021.106391
[18]	United Nations Environment Programme, From Pollution to Solution, A global assessment of marine litter and plastic pollution, Nairobi, 2021.
[19]	I. Tongo, K. O. Achema, D. Okuonghae, Dual-level toxicity assessment of biodegradable pesticides to aquatic species, Ecol. Complex., 45 (2021), 1–15. https://doi.org/10.1016/j.ecocom.2021.100911 doi: 10.1016/j.ecocom.2021.100911
[20]	Environmental Protection Agency (EPA), Office of water office of science and technology engineering and analysis division(4303T), Washington, 2007.
[21]	Y. Wan, B. Hassard, N. Kazarinoff, Theory and Application of Hopf bifurcation, London Mathematical Society Lecture Note Series, Cambridge University Press, Cambridge, 1981.
[22]	P. Kumar, R. Cochard, S. Maneepifak, Aquatic final abundance and diversity in relation to synthetic and natural pesticides applications in rice fields of central Thailand, International Journal of Biodiversity Science, Ecosystem Services and Management, 10 (2014), 157–173. https://doi.org/10.1080/21513732.2014.892029 doi: 10.1080/21513732.2014.892029
[23]	A. K. Agrawal, K. A. Anuj, A. W. Khan, The effect of an external toxicant on a biological species in case of deformity: a model, Model. Earth Syst. Env., 2 (2016), 1–8. https://doi.org/10.1007/s40808-016-0203-x doi: 10.1007/s40808-016-0203-x
[24]	C. E. Clark, T. G. Hallam, Non-autonomous logistic equation as models of population in a deteriorating environment, J. Theor. Biol., 93 (1982), 303–311. https://doi.org/10.1016/0022-5193(81)90106-5 doi: 10.1016/0022-5193(81)90106-5
[25]	B. Thomas, Ecotoxicology and environmental toxicology, an introduction, University of Gothenburg, 2012.
[26]	A. Chaturvedi, K. Ramesh, G. A. Vatsala, A mathematical approach to study the effect of pollutants/toxicants in aquatic environment, International Journal of Research - Granthalayah, 5 (2017), 33–38. https://doi.org/10.29121/granthaalayah.v5.i4RAST.2017.3299 doi: 10.29121/granthaalayah.v5.i4RAST.2017.3299
[27]	S. K. Arun, G. P. Satish, D. M. Ashokrao, Microbial degradation of pesticides: a review, Afr. J. Microbiol. Res., 11 (2017), 992–1012. https://doi.org/10.5897/AJMR2016.8402 doi: 10.5897/AJMR2016.8402
[28]	S. Jain, A. A. Alderremy, S. Aly, P. Agarwal, S. Deniz, A new analysis of a partial differential equation arising in biology and population genetics via semi analytical techniques, Physica A, 542 (2020), 122769.
[29]	P. Agarwal, A. U. Rehman, R. Singh, Modeling, analysis and prediction of new variants of covid-19 and dengue co-infection on complex network, Chaos, Solitons and Fractals, 150 (2021), 1–19. https://doi.org/10.1016/j.chaos.2021.111008 doi: 10.1016/j.chaos.2021.111008
[30]	S. F. Ahmed, P. Agarwal, A. B. M. S. Ali, S. M. Muyeen, S. M. E. K. Chowdhury, M. Forkan, Modeling the sars-cov-2 parallel transmission dynamics: Asymptomatic and symptomatic pathways, Comput. Biol. Med., 143 (2022), 1–11. https://doi.org/10.1016/j.compbiomed.2022.105264 doi: 10.1016/j.compbiomed.2022.105264
[31]	M. Ruzhansky, D. F. M. Torres, P. Agarwal, J. J. Nieto, Analysis of infectious disease problems (covid-19) and their global impact, Springer: Singapore, 2021.
[32]	H. W. Dehne, E. C. Oerke, Safeguarding production losses in major crops and the role of crop protection, Crop Prot., 23 (2004), 275–285. https://doi.org/10.1016/j.cropro.2003.10.001 doi: 10.1016/j.cropro.2003.10.001
[33]	S. Rahman, Pesticide consumption and productivity and the potential of IPM in Bangladesh, Sci. Total Environ., 445 (2013), 48–56. https://doi.org/10.1016/j.scitotenv.2012.12.032 doi: 10.1016/j.scitotenv.2012.12.032

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:*
© 2022 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)