# A numerical technique for solving variable order time fractional differential-integro equations

Mohammadhossein Derakhshan

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mohammadhossein Derakhshan. A numerical technique for solving variable order time fractional differential-integro equations. Communications in Mathematics, 2023, Volume 32 (2024), Issue 1, 10.46298/cm.ABCD. hal-03914959v2

HAL Id: hal-03914959

https://hal.science/hal-03914959v2

Submitted on 6 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# A numerical technique for solving variable order time fractional differential-integral equations

Mohammad Hossein Derakhshan

**Abstract.** In this manuscripts, we consider the coupled differential-integral equations including the variable-order Caputo fractional operator. To solve numerically these type of equations, we apply the shifted Jacobi-Gauss collocation scheme. Using this numerical method a system of algebraic equations is constructed. We solve this system with a recursive method in the nonlinear case and we solve it in linear case with algebraic formulas. Finally, for the high performance of the suggested method three Examples are illustrated.

## 1 Introduction

A coupled differential-integral equation including the Caputo fractional operator is introduced as follows:

$$ \underbrace{C_D^{\varrho_1(x)}[p_1(x)] + p'_1(x) + \int_0^x p_1(v)dv}_{\text{the initial condition for this types of equation is } p_1(0)=p'_1(0)=0} = \theta_1(x), \quad (1) $$

$$ \underbrace{C_D^{\varrho_2(x)}[p_2(x)] + p'_2(x) + \int_0^x p_2(v)dv}_{\text{the initial condition for this types of equation is } p_2(0)=p'_2(0)=0}, = \theta_2(x), \quad (2) $$

where $C_D^{\varrho_1(x)}, C_D^{\varrho_2(x)}$ are the Caputo fractional operators of variable orders $0 < \varrho_i(x) \le 1$, given by:

$$ C_D^{\varrho_i(x)} p(x) = \begin{cases} p'(x), & \varrho_i(x) = 1, \\ \frac{1}{\Gamma(1-\varrho_i(x))} \int_0^x (x-\nu)^{-\varrho_i(x)} p'(\nu) d\nu, & 0 < \varrho_i(x) \le 1. \end{cases} \quad (3) $$

MSC 2020: 26A33,35R11.

Keywords: Coupled differential-integral equation, Caputo fractional operator, Shifted fractional Jacobi collocation method, Variable-order.

Affiliation:

Zand institute of higher education, faculty of technology and engineering, Shiraz, Iran

Department of Industrial Engineering, Apadana Institute of Higher Education, Shiraz, Iran

E-mail: m.h.derakhshan.20@gmail.com
In the coupled differential-integral equations (1), (2), $p_1(x), p_2(x)$ and $\theta_1(x), \theta_2(x)$ are in-
distinct and definite respectively. The main reason for choosing these kind of differential-
integral equations is because of their role and applications in various fields as computational
and mathematical sciences [9], engineering sciences which are modelled using these types
of equations, physics and mathematical models in chemistry.

Several types of integral equations (IEs) are considered very important in various func-
tional analysis topics because of their essential role in engineering, physics, economics, and
natural sciences. Many real-life applications can be well described using IEs. For instance,
IEs play an important role in kinetic molecular, radiative, neutron transport, traffic, and
queuing theories [3]. Fractional calculus constitutes a powerful mathematical tool in terms
of analysis when it is applicable to the investigation of arbitrary-order integrals and deriva-
tives via the generalisations of integer-order differentiation and n-fold integration. Some
interesting research works concerning the relation between symmetry and generalised frac-
tional calculus properties have been conducted [10]. Recently, the study of the integrals
and derivatives of fractional order was limited to a purely mathematical context, but in
recent decades, they have been applied in various natural sciences and technology fields
such as theoretical physics, fluid mechanics, biology, image processing, and entropy theory
which reveal a great potential of these fractional integrals and derivatives in modelling
scientific phenomena [26]. Both of the theoretical and practical studies of the fractional
differential equations' element have received a global research interest. The main pur-
pose of this paper is to obtain the solution of the coupled differential-integral equations
are introduced in (1) and (2) by using a numerical method that this numerical method
is called a shifted Jacobi-Gauss collocation algorithm. This numerical method contains
a shifted Jacobi polynomials. Obtaining a solution by using different numerical methods
was studied by some authors in the mathematical field. In this paper a few numerical
methods are mentioned, for Example, in [11] was studied a numerical method based on
Chebyshev cardinal wavelets, the new algorithm base on the discretisation method to ob-
tain solution of the fractional functional integral equations of variable-order in [14] was
discussed, the explicit and implicit Euler methods [26], the Legendre wavelets method [4],
the numerical method base on Laplace and Sumudu transform methods [9], the Adams-
Bashforth-Moulton scheme [5], the optimisation method [23].

So, this article is divided into four sections as follows. In Section 2, Lemma and definitions of Jacobi polynomials, shifted Jacobi polynomials and their properties are introduced. In Section 3, approximation function and the numerical method to obtaining solution of the coupled differential-integral equation are described. In order to accuracy of the presented method four examples in Section 4 are described.

## 2 Some main lemma and definitions about the Jacobi polynomials and shifted Jacobi polynomials

This section deals with some definitions about the Jacobi polynomials and shifted
Jacobi polynomials and Lemma, which will be applied in next section.
**Definition 2.1.** Let $\mu, \nu > -1, x \in [-1, 1]$. Then the Jacobi polynomial of degree $n$ is defined by [7]:

$$ \mathfrak{P}_{n+1}^{\mu, \nu}(x) = (a_n^{\mu, \nu} x - b_n^{\mu, \nu}) \mathfrak{P}_n^{\mu, \nu}(x) - c_n^{\mu, \nu} \mathfrak{P}_{n-1}^{\mu, \nu}(x), \quad (4) $$

where for $n = 0, 1$ the functions $\mathfrak{P}_0^{\mu, \nu}(x)$ and $\mathfrak{P}_1^{\mu, \nu}(x)$ are defined by:

$$ \mathfrak{P}_0^{\mu, \nu}(x) = 1, \quad \mathfrak{P}_1^{\mu, \nu}(x) = \frac{(\mu + \nu + 2)x + \mu - \nu}{2}, \quad (5) $$

and the coefficients $a_n^{\mu, \nu}$, $b_n^{\mu, \nu}$, $c_n^{\mu, \nu}$ are given by:

$$ \begin{aligned} a_n^{\mu, \nu} &= \frac{(2n + \mu + \nu + 1)(2n + \mu + \nu + 2)}{(2n + 2)(n + \mu + \nu + 1)}, \\ b_n^{\mu, \nu} &= \frac{(\nu^2 - \mu^2)(2n + \mu + \nu + 1)}{(2n + 2)(n + \mu + \nu + 1)(2n + \mu + \nu)}, \\ c_n^{\mu, \nu} &= \frac{(n + \mu)(n + \nu)(2n + \mu + \nu + 2)}{(n + 1)(n + \mu + \nu + 1)(2n + \mu + \nu)}. \end{aligned} \quad (6) $$

Also, the Jacobi polynomial has a finite series as follows:

$$ \mathfrak{P}_n^{\mu, \nu}(x) = \frac{\Gamma(\mu + n + 1)}{n!\Gamma(\mu + \nu + n + 1)} \sum_{k=0}^{n} \binom{n}{k} \frac{\Gamma(\mu + \nu + k + n + 1)}{\Gamma(\mu + k + 1)} \left(\frac{x-1}{2}\right)^k, \quad (7) $$

where $\Gamma(\cdot)$ is the gamma function.

**Definition 2.2.** Let $\mu, \nu > -1$. Then the shifted Jacobi polynomial of degree $n$ is defined by [7]:

$$ \mathbb{P}_n^{\mu, \nu}(x) = \mathfrak{P}_n^{\mu, \nu}(\underbrace{2x-1}_{2x-1 \in [-1,1]}) \quad (8) $$

$$ \begin{aligned} &= \frac{\Gamma(\mu + n + 1)}{n!\Gamma(\mu + \nu + n + 1)} \sum_{k=0}^{n} \binom{n}{k} \frac{\Gamma(\mu + \nu + k + n + 1)}{\Gamma(\mu + k + 1)} (x-1)^k \\ &= \frac{\Gamma(\nu + n + 1)}{\Gamma(\mu + \nu + n + 1)} \sum_{k=0}^{n} (-1)^{n-k} \frac{\Gamma(\mu + \nu + k + n + 1)}{k!(n-k)!\Gamma(\nu + k + 1)} x^k. \end{aligned} \quad (9) $$

Also, the shifted Jacobi polynomial respect to the weight function $\omega^{\mu, \nu}(x) = (1-x)^{\mu}x^{\nu}$ is orthogonal, that the orthogonality condition is given as:

$$ \langle \mathbb{P}_n^{\mu, \nu}(x), \mathbb{P}_m^{\mu, \nu}(x) \rangle_{\omega^{\mu, \nu}(x)} = \int_0^1 \omega^{\mu, \nu}(x) \mathbb{P}_n^{\mu, \nu}(x) \mathbb{P}_m^{\mu, \nu}(x) dx = \begin{cases} 0, & n \neq m, \\ \frac{\Gamma(n+\mu+\nu+1)}{n!\Gamma(2n+\mu+\nu+1)}, & n=m, \end{cases} \quad (10) $$

where $\langle \cdot, \cdot \rangle$ is the inner product.
**Lemma 2.3.** Let $0 < \varrho_i \le 1, i = 1, 2$. Then for any $\zeta \ge -1$ the Caputo fractional operator of variable order $\varrho_i, i = 1, 2$ is given by [24]:

$$ C D_x^{\varrho_i(x)} x^{\zeta} = \left\{ \begin{array}{ll} \frac{\Gamma(\zeta+1)}{\Gamma(\zeta-\varrho_i(x)+1)} x^{\zeta-\varrho_i(x)}, & \zeta = \text{constant}, \\ \frac{d^n}{dx^n}(x^{\zeta}), & \zeta \ge -1, \\ \varrho_i(x) = m, m \in \mathbb{N}. & \end{array} \right. \quad (11) $$

Then, applying the relation (11) on (8), for $i = 1, 2$ we obtain:

$$ \begin{aligned} & C D_x^{\varrho_i(x)} \mathbb{P}_n^{\mu, \nu}(x) \\ &= \frac{\Gamma(\nu + n + 1)}{\Gamma(\mu + \nu + n + 1)} \sum_{k=0}^{n} (-1)^{n-k} \frac{\Gamma(k + 1)}{\Gamma(k - \varrho_i(x) + 1)} \frac{\Gamma(\mu + \nu + k + n + 1)}{k!(n-k)!\Gamma(\nu + k + 1)} x^{k-\varrho_i(x)} \\ &= \sum_{k=0}^{n} \Omega_{k,n}^{\mu, \nu} x^{k-\varrho_i(x)}, \end{aligned} \quad (12) $$

$$ \underbrace{\text{Caputo} D_x^{\varrho_i(x)} \mathbb{P}_n^{\mu, \nu}(x)}_{\text{for } \varrho_i(x)=1, m=1} = \frac{\Gamma(\nu + n + 1)}{\Gamma(\mu + \nu + n + 1)} \sum_{k=1}^{n} (-1)^{n-k} \frac{k\Gamma(\mu + \nu + k + n + 1)}{k!(n-k)!\Gamma(\nu + k + 1)} x^{k-1} = \sum_{k=1}^{n} \Psi_{k,n}^{\mu, \nu} x^{k-1}, \quad (13) $$

where

$$ \Omega_{k,n}^{\mu, \nu} = \frac{\Gamma(\nu + n + 1)}{\Gamma(\mu + \nu + n + 1)} \times (-1)^{n-k} \frac{\Gamma(k + 1)}{\Gamma(k - \varrho_i(x) + 1)} \times \frac{\Gamma(\mu + \nu + k + n + 1)}{k!(n-k)!\Gamma(\nu + k + 1)} $$

and

$$ \Psi_{k,n}^{\mu, \nu} = \frac{\Gamma(\nu + n + 1)}{\Gamma(\mu + \nu + n + 1)} \times (-1)^{n-k} \frac{k\Gamma(\mu + \nu + k + n + 1)}{k!(n-k)!\Gamma(\nu + k + 1)}. $$

### 3 The function approximation and proposed method algorithm

Due to the orthogonality condition is presented in (10) about the shifted Jacobi polynomial, any $p_1(x), p_2(x) \in L^2[0, 1]$ can be approximated in terms of the shifted Jacobi polynomial as follows:

$$ p_1(x) = \sum_{k=0}^{\infty} a_k \mathbb{P}_k^{\mu, \nu}(x), \quad (14) $$

$$ p_2(x) = \sum_{k=0}^{\infty} b_k \mathbb{P}_k^{\mu, \nu}(x), \quad (15) $$
considering the first $n$ sentence of the relations (14), (15), we have:

$$p_1(x) \approx p_{1,n}(x) = \sum_{k=0}^{n} a_k \mathbb{P}_k^{\mu,\nu}(x), \qquad (16)$$

$$p_2(x) \approx p_{2,n}(x) = \sum_{k=0}^{n} b_k \mathbb{P}_k^{\mu,\nu}(x), \qquad (17)$$

where the coefficients $a_j, b_j$ are calculated as:

$$a_j = \langle \mathbb{P}_n^{\mu,\nu}(x), p_1(x) \rangle_{\omega^{\mu,\nu}(x)} = \frac{\int_0^1 \omega^{\mu,\nu}(x) \mathbb{P}_n^{\mu,\nu}(x) p_1(x) dx}{\| \mathbb{P}_n^{\mu,\nu}(x) \|^2}, \qquad (18)$$

$$b_j = \langle \mathbb{P}_n^{\mu,\nu}(x), p_2(x) \rangle_{\omega^{\mu,\nu}(x)} = \frac{\int_0^1 \omega^{\mu,\nu}(x) \mathbb{P}_n^{\mu,\nu}(x) p_2(x) dx}{\| \mathbb{P}_n^{\mu,\nu}(x) \|^2}. \qquad (19)$$

**Theorem 3.1.** Let $p_{1,n}$ be as an approximation of the function $p_1(x)$, that it satisfy in the Equation (1) and the following conditions for $\lambda_1, \lambda_2 \in (0, 1)$ hold:

$$\| p_{1,n}(x) - p_{1,n-1}(x) \| \le \lambda_1 \| p_{1,n-1}(x) - p_{1,n-2}(x) \| \qquad (20)$$

$$\| p'_{1,n}(x) - p'_{1,n-1}(x) \| \le \lambda_2 \| p'_{1,n-1}(x) - p'_{1,n-2}(x) \|. \qquad (21)$$

Then the sequence $p_{1,n}$ to $p_1(x)$ converges, that $p_1(x)$ is the exact solution of Equation (1).

*Proof.* Suppose $\mathbf{R}_n(x)$ be a residual function, which is defined as:

$$\mathbf{R}_n(x) = \mathcal{C} \mathbf{D}_x^{\rho_1(x)} [p_{1,n-1}(x)] + p'_{1,n-1}(x) + \int_0^x p_{1,n-1}(v) dv - \theta_1(x),$$

$$\mathbf{R}_{n-1}(x) = \mathcal{C} \mathbf{D}_x^{\rho_1(x)} [p_{1,n-2}(x)] + p'_{1,n-2}(x) + \int_0^x p_{1,n-2}(v) dv - \theta_1(x). \qquad (22)$$
Then, we have:

$$
\begin{align*}
& \| \mathbf{R}_n(x) - \mathbf{R}_{n-1}(x) \| \\
&= \| {}^{C}\mathrm{D}_x^{\varrho_1(x)} [p_{1,n-1}(x) - p_{1,n-2}(x)] + (p'_{1,n-1}(x) - p'_{1,n-2}(x)) \\
&\quad + \int_0^x (p_{1,n-1}(v) - p_{1,n-2}(v)dv) \| \\
&= \| \frac{1}{\Gamma(1-\varrho_1(x))} \int_0^x (x-v)^{-\varrho_1(x)} [p'_{1,n-1}(v) - p'_{1,n-2}(v)] dv \\
&\quad + (p'_{1,n-1}(x) - p'_{1,n-2}(x)) \\
&\quad + \int_0^x (p_{1,n-1}(v) - p_{1,n-2}(v)dv) \| \\
&\le \frac{1}{\Gamma(1-\varrho_1(x))} \int_0^x (x-v)^{-\varrho_1(x)} \| p'_{1,n-1}(v) - p'_{1,n-2}(v) \| dv \\
&\quad + \| p'_{1,n-1}(x) - p'_{1,n-2}(x) \| \\
&\quad + \int_0^x \| p_{1,n-1}(v) - p_{1,n-2}(v) \| dv \\
&= \frac{x^{1-\varrho_1(x)}}{\Gamma(2-\varrho_1(x))} \int_0^x \| p'_{1,n-1}(v) - p'_{1,n-2}(v) \| dv + \| p'_{1,n-1}(x) - p'_{1,n-2}(x) \| \\
&\quad + \int_0^x \| p_{1,n-1}(v) - p_{1,n-2}(v) \| dv. \tag{23}
\end{align*}
$$

Using (20), (21), we get:

$$
\begin{align*}
& \| \mathbf{R}_n(x) - \mathbf{R}_{n-1}(x) \| \\
&\le \frac{\lambda_2 x^{1-\varrho_1(x)}}{\Gamma(2-\varrho_1(x))} \int_0^x \| p'_{1,n-2}(v) - p'_{1,n-3}(v) \| dv + \lambda_2 \| p'_{1,n-2}(x) - p'_{1,n-3}(x) \| \\
&\quad + \lambda_1 \int_0^x \| p_{1,n-2}(v) - p_{1,n-3}(v) \| dv \\
&\vdots \\
&\le \frac{\lambda_2^{n-2} x^{1-\varrho_1(x)}}{\Gamma(2-\varrho_1(x))} \int_0^x \| p'_{1,1}(v) - p'_{1,0}(v) \| dv + \lambda_2^{n-2} \| p'_{1,1}(x) - p'_{1,0}(x) \| \\
&\quad + \lambda_1^{n-2} \int_0^x \| p_{1,1}(v) - p_{1,0}(v) \| dv \\
&\le \frac{\lambda_2^{n-2} x^{2-\varrho_1(x)}}{\Gamma(2-\varrho_1(x))} \max_{x \in [0,1]} \| p'_{1,1}(x) - p'_{1,0}(x) \| + \lambda_2^{n-2} \max_{x \in [0,1]} \| p'_{1,1}(x) - p'_{1,0}(x) \| \\
&\quad + \lambda_1^{n-2} \max_{x \in [0,1]} \| p_{1,1}(x) - p_{1,0}(x) \|. \tag{24}
\end{align*}
$$

So, $\| \mathbf{R}_n(x) - \mathbf{R}_{n-1}(x) \| \to 0$ when $n \to \infty$, then $\mathbf{R}_n(x)$ a Cauchy sequence in $L^2[0, 1]$ and
$L^2[0, 1]$ is complete, therefore the sequence $\mathbf{R}_n(x)$ is convergent that is

$$ \lim_{n \to \infty} \mathbf{R}_n(x) = \mathbf{R}(x), \qquad (25) $$

where

$$ \begin{aligned} \mathbf{R}_n(x) &= {}^C\mathbf{D}_x^{\varrho_1(x)}[p_{1,n-1}(x)] + p'_{1,n-1}(x) + \int_0^x p_{1,n-1}(\upsilon)\,d\upsilon - \theta_1(x), \\ \mathbf{R}(x) &= {}^C\mathbf{D}_x^{\varrho_1(x)}[p_1(x)] + p'_1(x) + \int_0^x p_1(\upsilon)\,d\upsilon - \theta_1(x). \end{aligned} \qquad (26) $$

The proof is complete.

<figure><img src="image_2.png" /></figure>

Like the Theorem 3.1 can be considered a similar process for the Equation (2).

**Theorem 3.2 ([17]).** Suppose $f(x)$ belongs to the space $\mathbf{H}^{\sigma}(0, 1)$, $\sigma \ge 0$ that $\mathbf{H}^{\sigma}(0, 1)$ is defined by:

$$ \mathbf{H}^{\sigma}(0, 1) = \left\{ f \in L^2(0, 1) : \| f \|_{\mathbf{H}^{\sigma}(0, 1)} = \left( \sum_{k=0}^{\sigma} \| f^{(k)} \|_{L^2(0, 1)}^2 \right)^{\frac{1}{2}} < \infty \right\}, \quad (27) $$

and $P_{k,M}$ as the best approximation of $f$ is considered. Then the following inequality is hold:

$$ \| f - P_{k,M} \|_{L^2(0,1)} \le cM^{-\sigma}2^{-k\sigma} \| f^{(\sigma)} \|_{L^2(0,1)}^2, \qquad (28) $$

$$ \| f - P_{k,M} \|_{\mathbf{H}^q(0,1)} \le cM^{2q-\frac{1}{2}-\sigma} 2^{k(q-\sigma)} \| f^{(\sigma)} \|_{L^2(0,1)}^2, \quad q \ge 1. \qquad (29) $$

**Theorem 3.3.** Suppose $p_1(x)$ belongs to the space $\mathbf{H}^{\sigma}(0, 1)$, $\sigma \ge 0$ and $p_N$ is the best approximation of $p_1(x)$. Then we have:

$$ \| \mathbb{E} \|_{L^2(0,1)} \le \frac{cN^{-\sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2}{\Gamma(2-\varrho_1(x))} + 2cN^{-\sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2, \qquad (30) $$

$$ \| \mathbb{E} \|_{L^2(0,1)} \le \frac{cN^{2q-\frac{1}{2}-\sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2}{\Gamma(2-\varrho_1(x))} + 2cN^{2q-\frac{1}{2}-\sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2, \quad q \ge 1, \qquad (31) $$

here $\mathbb{E} = p_1(x) - p_N(x)$ is considered.
Proof. Using (1), we obtain

$$
\begin{align*}
\| \mathbb{E} \|_{L^2(0,1)} &= \| {}^{C}D_x^{\varrho_1(x)} [p_1(x) - p_N(x)] + (p'_1(x) - p'_N(x)) \\
&\qquad + \int_0^x (p_1(\upsilon) - p_N(\upsilon)\,d\upsilon) \|_{L^2(0,1)} \\
&= \| \frac{1}{\Gamma(1-\varrho_1(x))} \int_0^x (x-\upsilon)^{-\varrho_1(x)} [p'_1(\upsilon) - p'_N(\upsilon)] \,d\upsilon \\
&\qquad + (p'_1(x) - p'_N(x)) + \int_0^x (p_1(\upsilon) - p_N(\upsilon)\,d\upsilon) \|_{L^2(0,1)} \\
&\leq \frac{1}{\Gamma(1-\varrho_1(x))} \| \int_0^x (x-\upsilon)^{-\varrho_1(x)} [p'_1(\upsilon) - p'_N(\upsilon)] \,d\upsilon \|_{L^2(0,1)} \\
&\qquad + \| p'_1(x) - p'_N(x) \|_{L^2(0,1)} + \| \int_0^x (p_1(\upsilon) - p_N(\upsilon)\,d\upsilon) \|_{L^2(0,1)} . \tag{32}
\end{align*}
$$

We use $\| f * g \|_{L^2(0,1)} \le \| f \|_1 \| g \|_{L^2(0,1)}$ for Eq.(32), we have:

$$
\begin{align*}
\| \mathbb{E} \|_{L^2(0,1)} &\le \frac{1}{\Gamma(1-\varrho_1(x))} \| \int_0^x (x-\upsilon)^{-\varrho_1(x)} [p'_1(\upsilon) - p'_N(\upsilon)] \,d\upsilon \|_{L^2(0,1)} \\
&\quad + \| p'_1(x) - p'_N(x) \|_{L^2(0,1)} + \| \int_0^x (p_1(\upsilon) - p_N(\upsilon)\,d\upsilon) \|_{L^2(0,1)} \\
&\le \frac{x^{1-\varrho_1(x)}}{\Gamma(2-\varrho_1(x))} \| p'_1(x) - p'_N(x) \|_{L^2(0,1)} \\
&\quad + \| p'_1(x) - p'_N(x) \|_{L^2(0,1)} + \| \int_0^x (p_1(\upsilon) - p_N(\upsilon)\,d\upsilon) \|_{L^2(0,1)} \\
&\le \frac{x^{1-\varrho_1(x)}}{\Gamma(2-\varrho_1(x))} \| p_1(x) - p_N(x) \|_{L^2(0,1)} \\
&\quad + \| p_1(x) - p_N(x) \|_{L^2(0,1)} + x \| p_1(x) - p_N(x) \|_{L^2(0,1)} \\
&\le \frac{\| p_1(x) - p_N(x) \|_{L^2(0,1)}}{\Gamma(2-\varrho_1(x))} + 2 \| p_1(x) - p_N(x) \|_{L^2(0,1)}, \tag{33}
\end{align*}
$$

so, we employ Eqs.(28) and (29) on the relation (33), we get:

$$
\| \mathbb{E} \|_{L^2(0,1)} \le \frac{c N^{-\sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2}{\Gamma(2 - \varrho_1(x))} + 2c N^{-\sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2, \quad (34)
$$

$$
\| \mathbb{E} \|_{L^2(0,1)} \le \frac{c N^{2q - \frac{1}{2} - \sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2}{\Gamma(2 - \varrho_1(x))} + 2c N^{2q - \frac{1}{2} - \sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2, \quad q \ge 1. \quad (35)
$$

The desired result is achieved.

<figure><img src="image_3.png" /></figure>

**Theorem 3.4.** Suppose $p_2(x)$ is belongs to the space $H^\sigma(0,1)$, $\sigma \ge 0$ and $p_{N'}$ is the best
approximation of $p_1(x)$. Then we have:

$$ \| \mathbb{E} \|_{L^2(0,1)} \le \frac{c N'^{-\sigma} \| p_1^{(\sigma)} \|_{L^2(0,1)}^2}{\Gamma(2 - \varrho_2(x))} + 2c N'^{-\sigma} \| p_2^{(\sigma)} \|_{L^2(0,1)}^2, \quad (36) $$

$$ \| \mathbb{E} \|_{L^2(0,1)} \le \frac{c N'^{2q-\frac{1}{2}-\sigma} \| p_2^{(\sigma)} \|_{L^2(0,1)}^2}{\Gamma(2 - \varrho_2(x))} + 2c N'^{2q-\frac{1}{2}-\sigma} \| p_2^{(\sigma)} \|_{L^2(0,1)}^2, \quad q \ge 1. \quad (37) $$

*Proof.* The proof of this theorem is similar to the proof of the Theorem 3.3. □

## 3.1 The shifted Jacobi-Gauss collocation algorithm

This subsection describes the proposed method for solving Eqs.(1), (2). By substituting the relations (16) and (17) in Eqs.(1) and (2), we obtain:

$$ \underbrace{\mathbb{C} \mathbb{D}_x^{\varrho_1(x)} \left[ \sum_{k=0}^{n} a_k \mathbb{P}_k^{\mu,\nu}(x) \right] + \frac{d}{dx} \left( \sum_{k=0}^{n} a_k \mathbb{P}_k^{\mu,\nu}(x) \right) + \int_0^x \sum_{k=0}^{n} a_k \mathbb{P}_k^{\mu,\nu}(v) dv}_{\text{the initial condition for this types of equation is } \sum_{k=0}^{n} a_k \mathbb{P}_k^{\mu,\nu}(0) = \sum_{k=0}^{n} a_k \mathbb{P}_k^{\mu,\nu}(0)=0,} = \theta_1(x), \quad (38) $$

$$ \underbrace{\mathbb{C} \mathbb{D}_x^{\varrho_2(x)} \left[ \sum_{k=0}^{n} b_k \mathbb{P}_k^{\mu,\nu}(x) \right] + \sum_{k=0}^{n} b_k \mathbb{P}_k^{\mu,\nu}(x) + \int_0^x \sum_{k=0}^{n} b_k \mathbb{P}_k^{\mu,\nu}(v) dv}_{\text{the initial condition for this types of equation is } \sum_{k=0}^{n} b_k \mathbb{P}_k^{\mu,\nu}(0) = \sum_{k=0}^{n} b_k \mathbb{P}_k^{\mu,\nu}(0)=0,} = \theta_2(x), \quad (39) $$

by applying Eqs.(12),(13) on the equations (38) and (39), we have:

$$ \sum_{k=0}^{n} \sum_{m=0}^{k} a_k \Omega_{m,k}^{\mu,\nu} x^{m-\varrho_1(x)} + \sum_{k=0}^{n} \sum_{m=1}^{k+1} a_k \Psi_{m,k}^{\mu,\nu} x^{m-1} + \sum_{k=0}^{n} \sum_{m=0}^{k} a_k \Upsilon_{m,k}^{\mu,\nu} x^{m+1} = \theta_1(x), \quad (40) $$

$$ \sum_{k=0}^{n} \sum_{m=0}^{k} b_k \Omega_{m,k}^{\mu,\nu} x^{m-\varrho_2(x)} + \sum_{k=0}^{n} \sum_{m=1}^{k+1} b_k \Psi_{m,k}^{\mu,\nu} x^{m-1} + \sum_{k=0}^{n} \sum_{m=0}^{k} b_k \Upsilon_{m,k}^{\mu,\nu} x^{m+1} = \theta_2(x), \quad (41) $$

where $\Upsilon_{m,k}^{\mu,\nu} = (-1)^{k-m} \frac{\Gamma(\nu+k+1)\Gamma(\mu+\nu+m+k+1)}{(m+1)!(k-m)!(\nu+m+1)}$. Now we calculate the equations (40) and (41) in points $x_{\mu,\nu,i} = \frac{x_{\mu,\nu,i+1}}{2}$ which is the nodes of the standard Jacobi-Gauss interpolation in
the interval [-1, 1] and their introduced in [7]. Then for $i = 1, 2, \dots, n$ we get:

$$ \sum_{k=0}^{n} \sum_{m=0}^{k} a_k \Omega_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m-\varrho_1(x)} + \sum_{k=0}^{n} \sum_{m=1}^{k+1} a_k \Psi_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m-1} + \sum_{k=0}^{n} \sum_{m=0}^{k} a_k \Upsilon_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m+1} = \theta_1(x_{\mu,\nu,i}), \\ \sum_{k=0}^{n} a_k \mathbb{P}_k^{\mu,\nu}(0) = \sum_{k=0}^{n} a_k \mathbb{P}_k'^{\mu,\nu}(0) = 0. \qquad (42) $$

$$ \sum_{k=0}^{n} \sum_{m=0}^{k} b_k \Omega_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m-\varrho_2(x)} + \sum_{k=0}^{n} \sum_{m=1}^{k+1} b_k \Psi_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m-1} + \sum_{k=0}^{n} \sum_{m=0}^{k} b_k \Upsilon_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m+1} = \theta_2(x_{\mu,\nu,i}), \\ \sum_{k=0}^{n} b_k \mathbb{P}_k^{\mu,\nu}(0) = \sum_{k=0}^{n} b_k \mathbb{P}_k'^{\mu,\nu}(0) = 0. \qquad (43) $$

By solving these equations together with the initial conditions, the uncertain coefficients $a_k, b_k$ are obtained. When the functions $\theta_1(x), \theta_2(x)$ are nonlinear, in this case the equations (42) and (43) changes to the following equations:

$$ \sum_{k=0}^{n} \sum_{m=0}^{k} a_k \Omega_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m-\varrho_1(x)} + \sum_{k=0}^{n} \sum_{m=1}^{k+1} a_k \Psi_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m-1} + \sum_{k=0}^{n} \sum_{m=0}^{k} a_k \Upsilon_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m+1} \\ = \theta_1(x_{\mu,\nu,i}, p_1(x_{\mu,\nu,i}), p_1'(x_{\mu,\nu,i})), \\ \sum_{k=0}^{n} a_k \mathbb{P}_k^{\mu,\nu}(0) = \sum_{k=0}^{n} a_k \mathbb{P}_k'^{\mu,\nu}(0) = 0. \qquad (44) $$

$$ \sum_{k=0}^{n} \sum_{m=0}^{k} b_k \Omega_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m-\varrho_2(x)} + \sum_{k=0}^{n} \sum_{m=1}^{k+1} b_k \Psi_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m-1} + \sum_{k=0}^{n} \sum_{m=0}^{k} b_k \Upsilon_{m,k}^{\mu,\nu} x_{\mu,\nu,i}^{m+1} \\ = \theta_2(x_{\mu,\nu,i}, p_2(x_{\mu,\nu,i}), p_2'(x_{\mu,\nu,i})), \\ \sum_{k=0}^{n} b_k \mathbb{P}_k^{\mu,\nu}(0) = \sum_{k=0}^{n} b_k \mathbb{P}_k'^{\mu,\nu}(0) = 0. \qquad (45) $$

By applying a recursive method on the relations (44) and (45) can be obtained the uncertain coefficients $a_k, b_k$.

# 4 Numerical Examples

This section examines four examples by using the proposed method to show its performance and accuracy.
<figure><img src="image_4.png" /><figcaption>Figure 1: Numerical and exact solutions of the Example 1 with $\rho_1(x) = 1 - 0.001x$ and for the values of the parameters $n, \mu = \frac{1}{2}, \nu = \frac{1}{4}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 2: Numerical and exact solutions of the Example 1 with $\rho_2(x) = 1 - \sin^2(x)$ and for the values of the parameters $n, \mu = \frac{1}{2}, \nu = \frac{1}{4}$.</figcaption></figure>
<table><caption>Table 1: The absolute error for the values of the parameters $n = \mu = \frac{1}{2}$, $\nu = \frac{1}{4}$ and $\varrho_1(x) = 1 - 0.001x$ for Example 1.</caption><thead><tr><th></th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th></tr><tr><th>x</th><th>n=8</th><th>n=16</th><th>n=32</th><th>n=64</th></tr></thead><tbody><tr><td>0.1</td><td>0.96485e-15</td><td>0.99077e-16</td><td>1.01719e-17</td><td>1.04410e-18</td></tr><tr><td>0.2</td><td>0.86589e-15</td><td>0.88993e-16</td><td>0.91443e-17</td><td>0.93940e-18</td></tr><tr><td>0.3</td><td>0.77425e-15</td><td>0.79650e-16</td><td>0.81918e-17</td><td>0.84231e-18</td></tr><tr><td>0.4</td><td>0.68954e-15</td><td>0.71009e-16</td><td>0.73105e-17</td><td>0.75244e-18</td></tr><tr><td>0.5</td><td>0.61139e-15</td><td>0.63033e-16</td><td>0.64967e-17</td><td>0.66940e-18</td></tr><tr><td>0.6</td><td>0.53947e-15</td><td>0.55688e-16</td><td>0.57467e-17</td><td>0.59284e-18</td></tr><tr><td>0.7</td><td>0.47343e-15</td><td>0.48940e-16</td><td>0.50573e-17</td><td>0.52242e-18</td></tr><tr><td>0.8</td><td>0.41297e-15</td><td>0.42757e-16</td><td>0.44252e-17</td><td>0.45780e-18</td></tr><tr><td>0.9</td><td>0.35778e-15</td><td>0.37110e-16</td><td>0.38473e-17</td><td>0.39869e-18</td></tr></tbody></table>

<table><caption>Table 2: The absolute error for the values of the parameters $n = \mu = \frac{1}{2}$, $\nu = \frac{1}{4}$ and $\varrho_2(x) = 1 - \sin^2(x)$ for Example 1.</caption><thead><tr><th></th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th></tr><tr><th>x</th><th>n=8</th><th>n=16</th><th>n=32</th><th>n=64</th></tr></thead><tbody><tr><td>0.1</td><td>1.02256e-15</td><td>1.12266e-17</td><td>1.22921e-19</td><td>1.34242e-21</td></tr><tr><td>0.2</td><td>1.46252e-17</td><td>1.58973e-19</td><td>1.72428e-21</td><td>1.86639e-23</td></tr><tr><td>0.3</td><td>0.68232e-17</td><td>0.75878e-19</td><td>0.84084e-21</td><td>0.92869e-23</td></tr><tr><td>0.4</td><td>0.42820e-17</td><td>0.48438e-19</td><td>0.54532e-21</td><td>0.61123e-23</td></tr><tr><td>0.5</td><td>0.24707e-17</td><td>0.28622e-19</td><td>0.32932e-21</td><td>0.37658e-23</td></tr><tr><td>0.6</td><td>0.12620e-17</td><td>0.15145e-19</td><td>0.17989e-21</td><td>0.21170e-23</td></tr><tr><td>0.7</td><td>0.05314e-17</td><td>0.06759e-19</td><td>0.08446e-21</td><td>0.10393e-23</td></tr><tr><td>0.8</td><td>0.01572e-17</td><td>0.02239e-19</td><td>0.03073e-21</td><td>0.04091e-23</td></tr><tr><td>0.9</td><td>0.00196e-17</td><td>0.00383e-19</td><td>0.00663e-21</td><td>0.01053e-23</td></tr></tbody></table>

## 4.1 Example 1.

We consider the following coupled differential-integral equation of order $\varrho_1(x) = 1 - 0.001x$ and $\varrho_2(x) = 1 - \sin^2(x)$:

$$ C_D^{\varrho_1(x)}[p_1(x)] + p'_1(x) + \int_0^x p_2(\nu)d\nu = \theta_1(x), \\ p_1(0) = p'_1(0) = 0, \tag{46} $$

$$ C_D^{\varrho_2(x)}[p_2(x)] + p'_2(x) + \int_0^x p_1(\nu)d\nu = \theta_2(x), \\ p_2(0) = p'_2(0) = 0, \tag{47} $$

where

$$ \theta_1(x) = x^2 \cosh(x) - 2t \sinh(x) + 2 \cosh(x) + 2(e^x + 2xe^x - 1) $$
<figure><img src="image_6.png" /><figcaption>Figure 3: Numerical and exact solutions of the Example 2 with $\rho_1(x) = \sin(\frac{x}{3})$ and for the values of the parameters $n, \mu = \frac{1}{2}, \nu = \frac{1}{4}$.</figcaption></figure>

and

$$ \theta_2(x) = 4x \cosh(x) + 2x^2 \cosh(x) + e^x(x - 1). $$

The analytical solutions for these questions are $p_1(x) = xe^x - x$, $p_2(x) = x^2 \sinh(x)$. This example is solved by the proposed method and the results of the approximate and exact solutions are shown in the Figures 1 and 2. In the Tables 1 and 2, the absolute error between the approximate and analytical solutions are displayed.

<table><caption>Table 3: The absolute error for the values of the parameters $n = \mu = \frac{1}{2}$, $\nu = \frac{1}{4}$ and $\rho_1(x) = \sin(\frac{x}{3})$ for Example 2.</caption><thead><tr><th></th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th></tr><tr><th>x</th><th>n=8</th><th>n=16</th><th>n=32</th><th>n=64</th></tr></thead><tbody><tr><td>0.1</td><td>0.20113e-15</td><td>0.20639e-16</td><td>0.21173e-17</td><td>0.21717e-18</td></tr><tr><td>0.2</td><td>0.18095e-15</td><td>0.18586e-16</td><td>0.19087e-17</td><td>0.19595e-18</td></tr><tr><td>0.3</td><td>0.16210e-15</td><td>0.16669e-16</td><td>0.17136e-17</td><td>0.17611e-18</td></tr><tr><td>0.4</td><td>0.14456e-15</td><td>0.14883e-16</td><td>0.15317e-17</td><td>0.15760e-18</td></tr><tr><td>0.5</td><td>0.12827e-15</td><td>0.13223e-16</td><td>0.13626e-17</td><td>0.14037e-18</td></tr><tr><td>0.6</td><td>0.11319e-15</td><td>0.11685e-16</td><td>0.12058e-17</td><td>0.12439e-18</td></tr><tr><td>0.7</td><td>0.09928e-15</td><td>0.10265e-16</td><td>0.10609e-17</td><td>0.10960e-18</td></tr><tr><td>0.8</td><td>0.08649e-15</td><td>0.08959e-16</td><td>0.09275e-17</td><td>0.09598e-18</td></tr><tr><td>0.9</td><td>0.07479e-15</td><td>0.07762e-16</td><td>0.08051e-17</td><td>0.08347e-18</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 4: Numerical and exact solutions of the Example 2 with $\rho_2(x) = \frac{x}{3}$ and for the values of the parameters $n, \mu = \frac{1}{2}, \nu = \frac{1}{4}$.</figcaption></figure>

## 4.2 Example 2.

Consider the following coupled differential-integral equation with $\rho_1(x) = \sin(\frac{x}{3})$ and $\rho_2(x) = \frac{x}{3}$:

$$ 
\begin{aligned}
{}^C\mathbf{D}_x^{\rho_1(x)}[p_1(x)] &+ p'_1(x) + \int_0^x p_1(\nu)d\nu = 4x(x+1) + 2x^2 + \frac{x^4}{4} + \frac{x^3}{3}, \\
p_1(0) & = p'_1(0) = 0,
\end{aligned}
\quad (48) $$

$$ 
\begin{aligned}
{}^C\mathbf{D}_x^{\rho_2(x)}[p_2(x)] &+ p'_2(x) + \int_0^x p_2(\nu)d\nu = 4x(x-1) + 2x^x + \frac{x^4}{4} - \frac{x^3}{3}, \\
p_2(0) & = p'_2(0) = 0,
\end{aligned}
\quad (49) $$

the exact solutions of Eqs.(48),(49) are $p_1(x) = x^2(x+1)$, $p_2(x) = x^2(x-1)$. This example is solved by the proposed method and the results of the approximate and exact solutions are shown in the Figures 3 and 4. In the Tables 3 and 4, the absolute error between the approximate and analytical solutions are displayed.

## 4.3 Example 3.

Next, we describe the nonlinear coupled differential-integral equation

$$ 
\begin{aligned}
{}^C\mathbf{D}_x^{\rho_1(x)}[p_1(x)] &+ p'_1(x) + \int_0^x p_1(\nu)d\nu = 2x^2 + \frac{x^4}{12} + \frac{x^9}{27} - p_1^3(x), \\
p_1(0) & = p'_1(0) = 0,
\end{aligned}
\quad (50) $$

$$ 
\begin{aligned}
{}^C\mathbf{D}_x^{\rho_2(x)}[p_2(x)] &+ p'_2(x) + \int_0^x p_2(\nu)d\nu = 4x + \frac{7x^3}{3} - p'_2(x)p_2(x), \\
p_2(0) & = p'_2(0) = 0,
\end{aligned}
\quad (51) $$
<figure><img src="image_8.png" /><figcaption>Figure 5: Numerical and exact solutions of the Example 3 with $\rho_1(x) = e^x - x^2$ and for the values of the parameters $n, \mu = \frac{1}{2}, \nu = \frac{1}{4}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Numerical and exact solutions of the Example 3 with $\rho_2(x) = \sin^2(x)$ and for the values of the parameters $n = \mu = \frac{1}{2}, \nu = \frac{1}{4}$.</figcaption></figure>
<table><caption>Table 4: The absolute error for the values of the parameters $n, \mu = \frac{1}{2}, \nu = \frac{1}{4}$ and $\varrho_2(x) = \frac{x}{3}$ for Example 2.</caption><thead><tr><th></th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th></tr><tr><th>x</th><th>n=8</th><th>n=16</th><th>n=32</th><th>n=64</th></tr></thead><tbody><tr><td>0.1</td><td>2.18198e-15</td><td>2.26887e-16</td><td>2.35518e-17</td><td>2.44075e-18</td></tr><tr><td>0.2</td><td>1.83178e-15</td><td>1.91941e-16</td><td>2.00709e-17</td><td>2.09466e-18</td></tr><tr><td>0.3</td><td>1.48495e-15</td><td>1.57079e-16</td><td>1.65732e-17</td><td>1.74436e-18</td></tr><tr><td>0.4</td><td>1.15160e-15</td><td>1.23312e-16</td><td>1.31596e-17</td><td>1.39996e-18</td></tr><tr><td>0.5</td><td>0.84181e-15</td><td>0.91649e-16</td><td>0.99312e-17</td><td>1.07154e-18</td></tr><tr><td>0.6</td><td>0.56569e-15</td><td>0.63101e-16</td><td>0.69891e-17</td><td>0.76923e-18</td></tr><tr><td>0.7</td><td>0.33335e-15</td><td>0.38678e-16</td><td>0.44342e-17</td><td>0.50311e-18</td></tr><tr><td>0.8</td><td>0.15489e-15</td><td>0.19390e-16</td><td>0.23676e-17</td><td>0.28329e-18</td></tr><tr><td>0.9</td><td>0.04040e-15</td><td>0.06247e-16</td><td>0.08902e-17</td><td>0.11987e-18</td></tr></tbody></table>

<table><caption>Table 5: The absolute error for the values of the parameters $n = \mu = \frac{1}{2}, \nu = \frac{1}{4}$ and $\varrho_1(x) = e^x - x^2$ for Example 3.</caption><thead><tr><th></th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th></tr><tr><th>x</th><th>n=8</th><th>n=16</th><th>n=32</th><th>n=64</th></tr></thead><tbody><tr><td>0.1</td><td>1.35242e-15</td><td>1.42118e-16</td><td>1.49223e-17</td><td>1.56560e-18</td></tr><tr><td>0.2</td><td>1.09957e-15</td><td>1.15953e-16</td><td>1.22164e-17</td><td>1.28592e-18</td></tr><tr><td>0.3</td><td>0.88038e-15</td><td>0.93215e-16</td><td>0.98592e-17</td><td>1.04171e-18</td></tr><tr><td>0.4</td><td>0.69244e-15</td><td>0.73662e-16</td><td>0.78265e-17</td><td>0.83056e-18</td></tr><tr><td>0.5</td><td>0.53335e-15</td><td>0.57055e-16</td><td>0.60944e-17</td><td>0.65006e-18</td></tr><tr><td>0.6</td><td>0.40071e-15</td><td>0.43153e-16</td><td>0.46388e-17</td><td>0.49781e-18</td></tr><tr><td>0.7</td><td>0.29212e-15</td><td>0.31715e-16</td><td>0.34356e-17</td><td>0.37141e-18</td></tr><tr><td>0.8</td><td>0.20516e-15</td><td>0.22501e-16</td><td>0.24609e-17</td><td>0.26845e-18</td></tr><tr><td>0.9</td><td>0.13744e-15</td><td>0.15270e-16</td><td>0.16905e-17</td><td>0.18652e-18</td></tr></tbody></table>

with $\varrho_1(x) = e^x - x^2$ and $\varrho_2(x) = \sin^2(x)$. The exact solutions of Eqs.(50),(51) are $p_1(x) = \frac{1}{3}x^3, p_2(x) = x^2$. This example is solved by the proposed method and the results of the approximate and exact solutions are shown in the Figures 5 and 6. In the Tables 5 and 6, the absolute error between the approximate and analytical solutions are displayed.

## 5 Conclusion

This article is focused on the coupled differential-integral equation including the Caputo fractional operator of variable-orders and it is proposed a numerical method based on the shifted Jacobi-Gauss collocation scheme to obtain the solution of the coupled differential-integral equation. Also, in this paper about the convergence and an upper bound on the error are discussed. Some numerical Examples have been showed in order to display the high exactness of the suggested method.
<table><caption>Table 6: The absolute error for the values of the parameters $n = \mu = \frac{1}{2}$, $\nu = \frac{1}{4}$ and $\varrho_2(x) = \sin^2(x)$ for Example 3.</caption><thead><tr><th></th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th><th colspan="2">|p<sub>2</sub>(x) - p<sub>2,n</sub>(x)|</th></tr><tr><th>x</th><th>n=8</th><th>n=16</th><th>n=32</th><th>n=64</th></tr></thead><tbody><tr><td>0.1</td><td>1.96991e-15</td><td>2.02502e-16</td><td>2.08087e-17</td><td>2.13749e-18</td></tr><tr><td>0.2</td><td>1.75711e-15</td><td>1.80917e-16</td><td>1.86199e-17</td><td>1.91558e-18</td></tr><tr><td>0.3</td><td>1.55647e-15</td><td>1.60549e-16</td><td>1.65527e-17</td><td>1.70581e-18</td></tr><tr><td>0.4</td><td>1.36799e-15</td><td>1.41398e-16</td><td>1.46072e-17</td><td>1.50822e-18</td></tr><tr><td>0.5</td><td>1.19168e-15</td><td>1.23462e-16</td><td>1.27832e-17</td><td>1.32277e-18</td></tr><tr><td>0.6</td><td>1.02752e-15</td><td>1.06742e-16</td><td>1.10808e-17</td><td>1.14950e-18</td></tr><tr><td>0.7</td><td>0.87552e-15</td><td>0.91238e-16</td><td>0.95000e-17</td><td>0.98838e-18</td></tr><tr><td>0.8</td><td>0.73568e-15</td><td>0.76950e-16</td><td>0.80408e-17</td><td>0.83942e-18</td></tr><tr><td>0.9</td><td>0.60800e-15</td><td>0.63878e-16</td><td>0.67032e-17</td><td>0.70262e-18</td></tr></tbody></table>

## References

[1] Abdou M. A.: Fredholm–Volterra integral equation of the first kind and contact problem. *Applied Mathematics and Computation* 125 (2-3) (2002) 177–193.

[2] Amiri S, Hajipour M, Baleanu D.: A spectral collocation method with piecewise trigonometric basis functions for nonlinear Volterra–Fredholm integral equations. *Applied Mathematics and Computation* 370 (2020) 124915.

[3] Argyros I. K.: Quadratic equations and applications to Chandrasekhar’s and related equations. *Bulletin of the Australian Mathematical Society* 32 (2) (1985) 275–292.

[4] Chen Y. M., Wei Y. Q., Liu D. Y., Yu H.: Numerical solution for a class of nonlinear variable order fractional differential equations with Legendre wavelets. *Applied Mathematics Letters* (2015) 83–88.

[5] Coronel-Escamilla A., Gómez-Aguilar J. F.: A novel predictor-corrector scheme for solving variable-order fractional delay differential equations involving operators with Mittag-Leffler kernel. *Discrete & Continuous Dynamical Systems-S* 13 (3) (2020) 561.

[6] Dabiri A., Butcher E. A., Nazari M.: Coefficient of restitution in fractional viscoelastic compliant impacts using fractional Chebyshev collocation. *Journal of Sound and Vibration* 388 (2017) 230–244.

[7] Doha E. H., Abdelkawy M. A., Amin A. Z., Lopes A. M.: On spectral methods for solving variable-order fractional integro-differential equations. *Computational and Applied Mathematics* 37 (3) (2018) 3937–3950.

[8] Ganji R. M., Jafari H.: Numerical solution of variable order integro-differential equations. *Advanced Mathematical Models & Applications* 4 (1) (2019) 64–69.

[9] Gómez-Aguilar J. F.: Analytical and Numerical solutions of a nonlinear alcoholism model via variable-order fractional differential equations. *Physics A: Statistical Mechanics and its Applications* 494 (2018) 52–75.

[10] Hariharan G., Kannan K.: Review of wavelet methods for the solution of reaction-diffusion problems in science and engineering. *Applied Mathematical Modelling* 38 (3) (2014) 799–813.
[11] Heydari M. H.: Chebyshev cardinal wavelets for nonlinear variable-order fractional quadratic integral equations. *Applied Numerical Mathematics* 144 (2019) 190–203.

[12] Holmåker K.: Global asymptotic stability for a stationary solution of a system of integro-differential equations describing the formation of liver zones. *SIAM Journal on Mathematical Analysis* 24 (1) (1993) 116–128.

[13] Kaucher E. W., Miranker W. L.: Self-validating numerics for function space problems: Computation with guarantees for differential and integral equations. *Elsevier* (2014).

[14] Keshi F. K., Moghaddam B. P., Aghili A.: A numerical approach for solving a class of variable-order fractional functional integral equations. *Computational and Applied Mathematics* 37 (4) (2018) 4821–4834.

[15] Li P.: Singular integral equations of convolution type with Cauchy kernel in the class of exponentially increasing functions. *Applied Mathematics and Computation* 344 (2019) 116–127.

[16] Li Y., Jiang M., Liu F.: Time fractional super-diffusion model and its application in peak-preserving smoothing. *Chemometrics and Intelligent Laboratory Systems* 175 (2018) 13–19.

[17] Mashayekhi S., Ordokhani Y., Razzaghi M.: A hybrid functions approach for the Duffing equation. *Physica Scripta* 88 (2) (2013) 025002.

[18] Moghaddam B. P., Machado J. A.: A computational approach for the solution of a class of variable-order fractional integro-differential equations with weakly singular kernels. *Fractional Calculus and Applied Analysis* 20 (4) (2017) 1023–1042.

[19] Mohebbi A., Saffarian M.: Implicit RBF Meshless Method for the Solution of Two-dimensional Variable Order Fractional Cable Equation. *Journal of Applied and Computational Mechanics* 6 (2) (2020) 235–247.

[20] Mousavi B. K., Heydari M. H.: Wilson wavelets method for solving nonlinear fractional Fredholm-Hammerstein integro-differential equations. *International Journal of Computer Mathematics* (2019) 1–3.

[21] Mvogo A., Ben-Bolie G. H., Kofané T. C.: Fractional nonlinear dynamics of DNA breathing. *Communications in Nonlinear Science and Numerical Simulation* 48 (2017) 258–269.

[22] Rong S. I.: Theory of stochastic differential equations with jumps and applications: mathematical and analytical techniques with applications to engineering. Springer Science & Business Media (2006).

[23] Soradi-Zeid S., Jahanshahi H., Yousefpour A., Bekiros S.: King algorithm: A novel optimization approach based on variable-order fractional calculus with application in chaotic financial systems. *Chaos, Solitons & Fractals* 132 (2020) 109569.

[24] Xie J., Wang T., Ren Z., Zhang J., Quan L.: Haar wavelet method for approximating the solution of a coupled system of fractional-order integral-differential equations. *Mathematics and Computers in Simulation* 163 (2019) 80–89.

[25] Xu K., Wang J.: Weighted fractional permutation entropy and fractional sample entropy for nonlinear Potts financial dynamics. *Physics Letters A* 381 (8) (2017) 767–779.
[26] Zhuang P., Liu F., Anh V., Turner I.: Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term. *SIAM Journal on Numerical Analysis* 47 (3) (2009) 1760–1781.

Received: May 15, 2021

Accepted for publication: August 14, 2021

Communicated by: Yong-Xin Guo