A Spectral Method for Functional-Order Diffusion–Wave Models Based on Generalized Chelyshkov Wavelets

1. Introduction

Unlike the classical derivative, the use of fractional derivatives in many real-world phenomena can help to capture information of the current state from all its previous states. With this non-local property, fractional differential equations (FDEs) are used in many practical problems, such as modeling heat transfer [1], earthquakes [2], physics and mechanics [3,4,5], chemistry and dynamics in biological tissues [6], engineering [7], control theory [8,9,10], and option pricing [11].

In recent decades, FDEs have been extended to variable-order FDEs (VFDEs). In this equation, the orders of the derivatives can be non-constant functions [12,13,14]. This generalization enhances the fractional derivative’s ability to capture non-local effects, thus, popularizes the applications of VFDEs (see [15] and references therein). The advantages of VFDEs in comparison with FDEs have been discussed in [16], while some results on the consistency of this kind of equation and the uniqueness of the analytical solutions of VFDEs have been given in [17,18,19]. The VFDEs occurring in the applications usually do not have a solution which can be expressed via elementary functions. Therefore, it is crucial to design efficient methods for approximating solutions of VFDEs. The available methods can be seen in [20,21,22,23].

Recently, wavelets have become an important tool for solving fractional calculus problems. For example, Legendre wavelets have been used efficiently in various methods [24]. Chebyshev wavelets demonstrated high accuracy in many methods [25,26]. Some other useful wavelets are the Bernoulli and Taylor wavelets [27,28].

In [29] and [30], fractional-order Laguerre and Legendre functions are introduced to improve the accuracy of solutions for FDEs that contain power terms of fractional-order. In addition, generalized Legendre and Bernoulli wavelets were given in [31] and [32]. In these methods, finding the Riemann-Liouville integral operator with functional-order(RLI-FO) of the bases was an important step. The RLI-FO of hybrid Bernoulli polynomials was derived in [33,34] by using Laplace transforms. The same argument was used for Taylor wavelets in [28,35]. However, the Laplace transform cannot be used to find the closed form formula for the RLI-FO of any generalized form in both cases of polynomials and wavelets.

In this article, we provide a numerical technique for the solution of FO-FDWEs given in [36] by

$$\begin{array}{cc}\hfill & D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+\sum _{i=1}^m{a}_i{D}_{\mathbf{t}}^{{\omega }_i(\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+a{\displaystyle \frac{\partial }{\partial \mathbf{x}}}p(\mathbf{x},\mathbf{t})+b{p}^n(\mathbf{x},\mathbf{t})=c{\displaystyle \frac{{\partial }^2}{\partial {\mathbf{x}}^2}}p(\mathbf{x},\mathbf{t})+q(\mathbf{x},\mathbf{t}),\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & p(\mathbf{x},0)=r_0\left(\mathbf{x}\right),p_{\mathbf{t}}(\mathbf{x},0)=r_1\left(\mathbf{x}\right),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & p(0,\mathbf{t})=s_0\left(\mathbf{t}\right),p(1,\mathbf{t})=s_1\left(\mathbf{t}\right),\hfill \end{array}$$

(3)

where $(\mathbf{x},\mathbf{t})\in {[0,1]}^2$; $a,b,c,a_i$ are constants; $1<{\omega }_m(\mathbf{x},\mathbf{t})<\dots <{\omega }_1(\mathbf{x},\mathbf{t})<\omega (\mathbf{x},\mathbf{t})\le 2$; $r_0\left(\mathbf{x}\right),r_1\left(\mathbf{x}\right),s_0\left(\mathbf{t}\right)$ and $s_1\left(\mathbf{t}\right)$ are known functions. Here, $D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}$ and $D_{\mathbf{t}}^{{\omega }_i(\mathbf{x},\mathbf{t})}$ denote the Caputo functional-order derivative (CFD) defined by

$$\begin{array}{cc}\hfill & D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\Gamma (2-\omega (\mathbf{x},\mathbf{t}\left)\right)}}{\int }_0^{\mathbf{t}}{(\mathbf{t}-s)}^{1-\omega (\mathbf{x},\mathbf{t})}{\displaystyle \frac{{\partial }^{2}p(\mathbf{x},s)}{\partial s^2}}ds,\hfill & 1<\omega (\mathbf{x},\mathbf{t})<2,\hfill \\ p_{\mathbf{tt}}(\mathbf{x},\mathbf{t}),\hfill & \omega (\mathbf{x},\mathbf{t})=2.\hfill \end{array}\right.\hfill \\ \hfill & D_{\mathbf{t}}^{{\omega }_i(\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})={\displaystyle \frac{1}{\Gamma (2-{\omega }_i(\mathbf{x},\mathbf{t}))}}{\int }_0^{\mathbf{t}}{(\mathbf{t}-s)}^{1-{\omega }_i(\mathbf{x},\mathbf{t})}{\displaystyle \frac{{\partial }^{2}p(\mathbf{x},s)}{\partial s^2}}ds,i=1,\dots ,m.\hfill \end{array}$$

In the current paper, we numerically solve the FO-FDWE given in Eqs. (1)-(3) by using FOCW. The FOCW technique provides several benefits, including the orthogonality and the spectral accuracy properties [37,38]. The regularized incomplete Beta functions are used to derive an exact value for the RLI-FO of order $\omega \left(x\right)$ for the FOCW. This exact value can be used to simplify Eqs. (1)-(3) into a system of algebraic equations. In general, this system can be computationally complex and require a large storage space for large systems. However, by using FOCW, the computational complexity can be reduced [39]. The structure of this article is arranged as follows. The RLI-FO, the CFD, and the regularized incomplete Beta functions are recalled in Section 2. In Section 3, we construct the FOCW and find the exact value of their RLI-FO by using the regularized incomplete Beta functions. In Section 4 and Section 5, we present the idea of the method and its error estimation. In Section 6, several examples are included. In addition, in Example 1, we will get the exact solution which was not obtained in [36]. In example2,example3,ex7, we will point out that we will obtain better solutions than other methods in the literature.

2. Preliminaries

This section presents the fundamental concepts of the RLI-FO and CFD, which serve as the theoretical basis for our numerical approach. Unless otherwise stated, we assume $p\left(x\right)=0$ for all $x<0$.

Definition 1. 

The RLI-FO of order $\omega \left(x\right)$ of a function $p\left(\mathbf{x}\right)$ is defined as [40]

$$\begin{array}{ccc}\hfill I^{\omega \left(\mathbf{x}\right)}p\left(\mathbf{x}\right)& =& {\displaystyle \frac{1}{\Gamma \left(\omega \right(\mathbf{x}\left)\right)}}{\int }_0^{\mathbf{x}}p\left(\mathbf{t}\right){(\mathbf{x}-\mathbf{t})}^{\omega \left(\mathbf{x}\right)-1}d\mathbf{t},\mathbf{x}\ge 0.\hfill \end{array}$$

(4)

The CFD of order $\omega \left(x\right)$ of a function $p\left(\mathbf{x}\right)$ is defined as [40]

$$\begin{array}{c}\hfill D^{\omega \left(\mathbf{x}\right)}p\left(\mathbf{x}\right)={\displaystyle \frac{1}{\Gamma (\nu -\omega (\mathbf{x}\left)\right)}}{\int }_0^{\mathbf{x}}{p}^{\left(\nu \right)}\left(\mathbf{t}\right){(\mathbf{x}-\mathbf{t})}^{\nu -\omega \left(\mathbf{x}\right)-1}d\mathbf{t},\nu -1<\omega \left(\mathbf{x}\right)\le \nu ,\nu \in \mathbb{N}.\end{array}$$

(5)

The CFD and the RLI-FO are linear operators. Moreover, we have

(i)

For $\gamma >-1$, we have

$$\begin{array}{c}\hfill I^{\omega \left(\mathbf{x}\right)}{\mathbf{x}}^{\gamma }={\displaystyle \frac{\Gamma (1+\gamma )}{\Gamma (1+\gamma +\omega (\mathbf{x}\left)\right)}}{x}^{\gamma +\omega \left(\mathbf{x}\right)},\end{array}$$

(ii)

For $\gamma >\nu -1$,

$$\begin{array}{c}\hfill D^{\omega \left(\mathbf{x}\right)}{\mathbf{x}}^{\gamma }={\displaystyle \frac{\Gamma (1+\gamma )}{\Gamma (1+\gamma -\omega (\mathbf{x}\left)\right)}}{\mathbf{x}}^{\gamma -\omega \left(\mathbf{x}\right)}.\end{array}$$

The Regularized Beta functions, denoted by $\mathbb{I}(\mathbf{x};c_1,c_2)$, are defined as [41]

$$\mathbb{I}(\mathbf{x};c_1,c_2)={\displaystyle \frac{\Gamma (c_1+c_2)}{\Gamma \left(c_1\right)\Gamma \left(c_2\right)}}{\int }_0^{\mathbf{x}}{\mathbf{t}}^{c_1-1}{\left(1-\mathbf{t}\right)}^{c_2-1}d\mathbf{t},c_1>0,c_2>0,0\le \mathbf{x}\le 1,$$

where

$$\Gamma \left(u\right)={\int }_0^{\infty }{e}^{-\mathbf{t}}{\mathbf{t}}^{u-1}d\mathbf{t},u>0$$

is the Gamma function.

By using the regularized incomplete Beta functions, we get the following lemma:

Lemma 1. 

Let $c\ge 0,\gamma \ge -1$, we have

$$\begin{array}{ccc}\hfill I^{\omega \left(\mathbf{x}\right)}\left({\mathbf{x}}^{\gamma }{\mu }_c\left(\mathbf{x}\right)\right)& =& {\displaystyle \frac{\Gamma (1+\gamma ){\mathbf{x}}^{\gamma +\omega \left(\mathbf{x}\right)}}{\Gamma (1+\gamma +\omega (\mathbf{x}\left)\right)}}\left[1-\mathbb{I}\left({\displaystyle \frac{c}{\mathbf{x}}};1+\gamma ,\omega \left(\mathbf{x}\right)\right)\right]{\mu }_c\left(\mathbf{x}\right),\hfill \end{array}$$

(6)

where

$$\begin{array}{c}\hfill {\mu }_c\left(\mathbf{x}\right)=\left\{\begin{array}{cc}0,& \mathrm{if}\mathbf{x}<c,\\ 1,& \mathrm{if}\mathbf{x}\ge c.\end{array}\right.\end{array}$$

Proof. 

For $\mathbf{x}\ge c$, by applying Eq. (4), the left-hand-side of Eq. (6) becomes

$$\begin{array}{ccc}\hfill I^{\omega \left(\mathbf{x}\right)}\left({\mathbf{x}}^{\gamma }{\mu }_c\left(\mathbf{x}\right)\right)& =& {\displaystyle \frac{1}{\Gamma \left(\omega \right(\mathbf{x}\left)\right)}}{\int }_0^{\mathbf{x}}{s}^{\gamma }{\mu }_c\left(s\right){(\mathbf{x}-s)}^{\omega \left(\mathbf{x}\right)-1}ds\hfill \\ & =& {\displaystyle \frac{1}{\Gamma \left(\omega \right(\mathbf{x}\left)\right)}}{\int }_0^{\mathbf{x}}{s}^{\gamma }{(\mathbf{x}-s)}^{\omega \left(\mathbf{x}\right)-1}ds-{\displaystyle \frac{1}{\Gamma \left(\omega \right(\mathbf{x}\left)\right)}}{\int }_0^c{s}^{\gamma }{(\mathbf{x}-s)}^{\omega \left(\mathbf{x}\right)-1}ds\hfill \end{array}$$

Thus

$$\begin{array}{ccc}\hfill I^{\omega \left(\mathbf{x}\right)}\left({\mathbf{x}}^{\gamma }{\mu }_c\left(\mathbf{x}\right)\right)& =& I^{\omega \left(\mathbf{x}\right)}\left({\mathbf{x}}^{\gamma }\right)-{\displaystyle \frac{1}{\Gamma \left(\omega \right(\mathbf{x}\left)\right)}}{\int }_0^c{\left({\displaystyle \frac{s}{\mathbf{x}}}\right)}^{\gamma }{\mathbf{x}}^{\gamma }{\left(1-{\displaystyle \frac{s}{\mathbf{x}}}\right)}^{\omega \left(\mathbf{x}\right)-1}{\mathbf{x}}^{\omega \left(\mathbf{x}\right)-1}\mathbf{x}d\left({\displaystyle \frac{s}{\mathbf{x}}}\right)\hfill \\ & =& I^{\omega \left(\mathbf{x}\right)}\left({\mathbf{x}}^{\gamma }\right)-{\displaystyle \frac{{\mathbf{x}}^{\gamma +\omega \left(\mathbf{x}\right)}}{\Gamma \left(\omega \right(\mathbf{x}\left)\right)}}{\int }_0^{{\displaystyle \frac{c}{\mathbf{x}}}}{\mathbf{t}}^{\gamma }{\left(1-\mathbf{t}\right)}^{\omega \left(\mathbf{x}\right)-1}d\mathbf{t}\hfill \\ & =& I^{\omega \left(\mathbf{x}\right)}\left({\mathbf{x}}^{\gamma }\right)-{\displaystyle \frac{{\mathbf{x}}^{\gamma +\omega \left(\mathbf{x}\right)}}{\Gamma \left(\omega \right(\mathbf{x}\left)\right)}}\mathbb{I}\left({\displaystyle \frac{c}{\mathbf{x}}};1+\gamma ,\omega \left(\mathbf{x}\right)\right){\displaystyle \frac{\Gamma (1+\gamma )\Gamma \left(\omega \right(\mathbf{x}\left)\right)}{\Gamma (1+\gamma +\omega (\mathbf{x}\left)\right)}}\hfill \\ & =& {\displaystyle \frac{\Gamma (1+\gamma )}{\Gamma (1+\gamma +\omega (\mathbf{x}\left)\right)}}{\mathbf{x}}^{\gamma +\omega \left(\mathbf{x}\right)}\left(1-\mathbb{I}\left({\displaystyle \frac{c}{\mathbf{x}}};1+\gamma ,\omega \left(\mathbf{x}\right)\right)\right).\hfill \end{array}$$

The lemma is then followed. □

3. Fractional-Order Chelyshkov Wavelets and Function Approximations in Two Dimensional Space

The Chelyshkov polynomials are defined as

$$\begin{array}{c}\hfill {\rho }_{m,}\left(\mathbf{x}\right)=\sum _{l=0}^{\mathbb{M}-m}{a}_{l,m}.{\mathbf{x}}^{m+l},m=0,1,\dots ,\mathbb{M},\end{array}$$

where

$$\begin{array}{c}\hfill a_{l,m}={(-1)}^l\left(\begin{array}{c}\mathbb{M}-m\\ l\end{array}\right)\left(\begin{array}{c}\mathbb{M}+m+l+1\\ \mathbb{M}-m\end{array}\right).\end{array}$$

Let $\mathbb{k},\mathbb{M}$ be nonnegative integers. The Chelyshkov wavelets ${\psi }_{n,m}\left(\mathbf{x}\right)$ for $n=0,1,\dots ,2^{\mathbb{k}}-1$ and $m=0,1,\dots ,\mathbb{M}$ are defined on $[0,1)$ as

$$\begin{array}{c}\hfill {\psi }_{n,m}\left(\mathbf{x}\right)=\left\{\begin{array}{cc}\sqrt{2^{\mathbb{k}}(2m+1)}{\rho }_{m,\mathbb{M}}(2^{\mathbb{k}}\mathbf{x}-n),\hfill & \mathrm{if}\mathbf{x}\in \left[{\displaystyle \frac{n}{2^{\mathbb{k}}}},{\displaystyle \frac{n+1}{2^{\mathbb{k}}}}\right),\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

The FOCW are obtained by using the transformation $\mathbf{x}={\mathbf{t}}^{\alpha }$ for some $\alpha >0$, as

$$\begin{array}{c}\hfill {\psi }_{n,m}^{\alpha }\left(\mathbf{t}\right)=\left\{\begin{array}{cc}\sqrt{2^{\mathbb{k}}(2m+1)}{\rho }_{m,\mathbb{M}}(2^{\mathbb{k}}{\mathbf{t}}^{\alpha }-n),\hfill & \mathrm{if}t\in \left[{\left({\displaystyle \frac{n}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}},{\left({\displaystyle \frac{n+1}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}\right),\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

The FOCW form an orthonormal basis for $L_{w^{\alpha }}^2[0,1]$, i.e.

$$\begin{array}{c}〈{\psi }_{n_1,m_1}^{\alpha }\left(t\right),{\psi }_{n_2,m_2}^{\alpha }\left(\mathbf{t}\right)〉={\int }_0^1{\psi }_{n_1,m_1}^{\alpha }\left(t\right){\psi }_{n_2,m_2}^{\alpha }\left(\mathbf{t}\right)w^{\alpha }\left(\mathbf{t}\right)dt={\delta }_{n_1,n_2}·{\delta }_{m_1,m_2},\end{array}$$

where $w^{\alpha }\left(\mathbf{t}\right)=\alpha {\mathbf{t}}^{\alpha -1}$ is the weight function. Therefore a function $p\left(\mathbf{t}\right)\in L_{w^{\alpha }}^2\left([0,1]\right)$ can be expressed as

$$p\left(\mathbf{t}\right)\approx \sum _{n=0}^{\mathbb{M}}\sum _{n=0}^{2^{\mathbb{k}}-1}{u}_{n,m}{\psi }_{n,m}^{\alpha }\left(\mathbf{t}\right)=C^T{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(t\right),$$

where

$$\begin{array}{c}{C}^T=\left[u_{0,0},\dots ,u_{0,\mathbb{M}},u_{1,0},\dots ,u_{1,\mathbb{M}},\dots ,u_{2^{\mathbb{k}}-1,0},\dots ,u_{2^{\mathbb{k}}-1,\mathbb{M}}\right],\end{array}$$

and

$$\begin{array}{c}{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(\mathbf{t}\right)={\left[{\psi }_{0,0}^{\alpha }\left(\mathbf{t}\right),\dots ,{\psi }_{0,\mathbb{M}}^{\alpha }\left(\mathbf{t}\right),{\psi }_{1,0}^{\alpha }\left(\mathbf{t}\right),\dots ,{\psi }_{1,\mathbb{M}}^{\alpha }\left(\mathbf{t}\right)\dots ,{\psi }_{2^{\mathbb{k}}-1,0}^{\alpha }\left(\mathbf{t}\right),\dots ,{\psi }_{2^{\mathbb{k}}-1,\mathbb{M}}^{\alpha }\left(\mathbf{t}\right)\right]}^T.\end{array}$$

The theorem below establishes a closed-form expression for the variable-order Riemann–Liouville integral of fractional-order Chelyshkov wavelets. This result serves as the foundation for constructing a novel numerical scheme for solving variable-order fractional diffusion–wave equations in the subsequent section.

Theorem 1. 

Let $\alpha >0$, we have

$$I^{\omega \left(\mathbf{t}\right)}{\psi }_{n,m}^{\alpha }\left(\mathbf{t}\right)=\left\{\begin{array}{cccc}\hfill & 0,\hfill & \hfill & \mathrm{if}\mathbf{t}\le c_1,\hfill \\ \hfill & A\left(\mathbf{t}\right),\hfill & \hfill & \mathrm{if}{c}_1\le \mathbf{t}<c_2,\hfill \\ \hfill & B\left(\mathbf{t}\right),\hfill & \hfill & \mathrm{if}\mathbf{t}\ge c_2,\hfill \end{array}\right.$$

where

$$\begin{array}{ccc}\hfill A\left(\mathbf{t}\right)& =& \sqrt{2^{\mathbb{k}}(2m+1)}\sum _{l=0}^{\mathbb{M}-m}\sum _{k=0}^{m+l}{a}_{l,m}\left(\begin{array}{c}m+l\\ k\end{array}\right)2^{\mathbb{k}l}{(-n)}^{m+l-k}\hfill \\ & & \times {\displaystyle \frac{\Gamma (1+\alpha k)}{\Gamma (1+\alpha k+\omega (\mathbf{t}\left)\right)}}{\mathbf{t}}^{\alpha k+\omega \left(\mathbf{t}\right)}\times \left[1-\mathbb{I}\left({\displaystyle \frac{c_1}{\mathbf{t}}};1+\alpha k,\omega \left(\mathbf{t}\right)\right)\right],\hfill \end{array}$$

$$\begin{array}{ccc}\hfill B\left(\mathbf{t}\right)& =& \sqrt{2^{\mathbb{k}}(2m+1)}\sum _{l=0}^{\mathbb{M}-m}\sum _{k=0}^{k+l}{a}_{l,m}\left(\begin{array}{c}m+l\\ k\end{array}\right)2^{\mathbb{k}k}{(-n)}^{n+l-k}\hfill \\ & & \times {\displaystyle \frac{\Gamma (1+\alpha k)}{\Gamma (1+\alpha k+\omega (\mathbf{t}\left)\right)}}{\mathbf{t}}^{\alpha k+\omega \left(\mathbf{t}\right)}\times \mathbb{I}\left({\displaystyle \frac{c_2}{\mathbf{t}}},{\displaystyle \frac{c_1}{\mathbf{t}}};1+\alpha k,\omega \left(\mathbf{t}\right)\right),\hfill \end{array}$$

and $c_1,c_2$ are two parameters defined by

$$\begin{array}{c}\hfill c_1={\left({\displaystyle \frac{n}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}},c_2={\left({\displaystyle \frac{n+1}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}.\end{array}$$

Proof. 

First, we reformulate the FOCW in the following form:

$$\begin{array}{ccc}\hfill {\psi }_{m,n}^{\alpha }\left(\mathbf{t}\right)& =& \sqrt{2^{\mathbb{k}}(2m+1)}{\rho }_{m,\mathbb{M}}\left(2^{\mathbb{k}}{\mathbf{t}}^{\alpha }-n\right)\left({\mu }_{c_1}\left(\mathbf{t}\right)-{\mu }_{c_2}\left(\mathbf{t}\right)\right)\hfill \\ & =& \sqrt{2^{\mathbb{k}}(2m+1)}\left(\sum _{l=0}^{\mathbb{M}-m}{a}_{l,m}{\left(2^{\mathbb{k}}{\mathbf{t}}^{\alpha }-n\right)}^{m+l}\left[{\mu }_{c_1}\left(\mathbf{t}\right)-{\mu }_{c_2}\left(\mathbf{t}\right)\right]\right)\hfill \\ & =& \sqrt{2^{\mathbb{k}}(2m+1)}\sum _{l=0}^{\mathbb{M}-m}\sum _{k=0}^{m+l}{a}_{l,m}\left(\begin{array}{c}m+l\\ k\end{array}\right)2^{\mathbb{k}k}{(-n)}^{m+l-k}\hfill \\ & & \times \left[{\mathbf{t}}^{\alpha k}{\mu }_{c_1}\left(\mathbf{t}\right)-{\mathbf{t}}^{\alpha k}{\mu }_{c_2}\left(\mathbf{t}\right)\right],\hfill \end{array}$$

Then

$$\begin{array}{ccc}\hfill I^{\omega \left(\mathbf{t}\right)}\left({\psi }_{n,m}^{\alpha }\left(\mathbf{t}\right)\right)& =& \sqrt{2^{\mathbb{k}}(2m+1)}\sum _{l=0}^{\mathbb{M}-m}\sum _{k=0}^{m+l}{a}_{l,m}\left(\begin{array}{c}m+l\\ k\end{array}\right)2^{\mathbb{k}k}{(-n)}^{m+l-k}\hfill \\ & \times & \left[I^{\omega \left(\mathbf{t}\right)}\left({}_{\alpha k}\mu _c^1\left(\mathbf{t}\right)\right)-I^{\omega \left(\mathbf{t}\right)}\left({\mathbf{t}}^{\alpha k}{\mu }_{c_2}\left(\mathbf{t}\right)\right)\right].\hfill \end{array}$$

After that, we apply Lemma 1 to obtain

$$\begin{array}{ccc}\hfill I^{\omega \left(\mathbf{t}\right)}\left({\psi }_{n,m}^{\alpha }\left(\mathbf{t}\right)\right)& =& \sqrt{2^{\mathbb{k}}(2m+1)}\sum _{l=0}^{\mathbb{M}-m}\sum _{k=0}^{m+l}{a}_{l,m}\left(\begin{array}{c}m+l\\ k\end{array}\right)2^{\mathbb{k}k}{(-n)}^{m+l-k}\hfill \\ & & \times \left[{\displaystyle \frac{\Gamma (1+\alpha k)}{\Gamma (1+\alpha k+\omega (\mathbf{t}\left)\right)}}{\mathbf{t}}^{\alpha k+\omega \left(\mathbf{t}\right)}\left(1-\mathbb{I}\left({\displaystyle \frac{c_1}{\mathbf{t}}};1+\alpha k,\omega \left(\mathbf{t}\right)\right)\right){\mu }_{c_1}\left(\mathbf{t}\right)\right.\hfill \\ & & \left.-{\displaystyle \frac{\Gamma (1+\alpha k)}{\Gamma (1+\alpha k+\omega (\mathbf{t}\left)\right)}}{\mathbf{t}}^{\alpha k+\omega \left(\mathbf{t}\right)}\left(1-\mathbb{I}\left({\displaystyle \frac{c_2}{\mathbf{t}}};1+\alpha k,\omega \left(\mathbf{t}\right)\right)\right){\mu }_{c_2}\left(\mathbf{t}\right)\right].\hfill \end{array}$$

If $\mathbf{t}<c_1$ then $I^{\omega \left(\mathbf{t}\right)}\left[{\psi }_{n,m}^{\alpha }\left(\mathbf{t}\right)\right]=0$. In case $c_1\le \mathbf{t}<c_2$, then

$$\begin{array}{ccc}\hfill I^{\omega \left(\mathbf{t}\right)}\left({\psi }_{n,m}^{\alpha }\left(\mathbf{t}\right)\right)& =& \sqrt{2^{\mathbb{k}}(2m+1)}\sum _{l=0}^{\mathbb{M}-m}\sum _{k=0}^{m+l}{a}_{l,m}\left(\begin{array}{c}m+l\\ k\end{array}\right)2^{\mathbb{k}k}{(-n)}^{m+l-k}\hfill \\ & & \times \left[{\displaystyle \frac{\Gamma (1+\alpha k)}{\Gamma (1+\alpha k+\omega (\mathbf{t}\left)\right)}}{\mathbf{t}}^{\alpha k+\omega \left(\mathbf{t}\right)}\left(1-\mathbb{I}\left({\displaystyle \frac{c_1}{\mathbf{t}}};1+\alpha k,\omega \left(\mathbf{t}\right)\right)\right)\right].\hfill \end{array}$$

For $\mathbf{t}\ge c_2$, we have

$$\begin{array}{c}\hfill I^{\omega \left(\mathbf{t}\right)}\left({\psi }_{n,m}^{\alpha }\left(\mathbf{t}\right)\right)=\sqrt{2^{\mathbb{k}}(2m+1)}\sum _{l=0}^{\mathbb{M}-m}\sum _{k=0}^{m+l}{a}_{l,m}\left(\begin{array}{c}m+l\\ k\end{array}\right)2^{\mathbb{k}k}{(-n)}^{m+l-k}\\ \hfill \times {\displaystyle \frac{\Gamma (1+\alpha k)}{\Gamma (1+\alpha k+\omega (\mathbf{t}\left)\right)}}{\mathbf{t}}^{\alpha k+\omega \left(\mathbf{t}\right)}\left[\mathbb{I}\left({\displaystyle \frac{c_2}{\mathbf{t}}};1+\alpha k,\omega \left(\mathbf{t}\right)\right)-\mathbb{I}\left({\displaystyle \frac{c_1}{\mathbf{t}}};1+\alpha k,\omega \left(\mathbf{t}\right)\right)\right]\\ \hfill =\sqrt{2^{\mathbb{k}}(2m+1)}\sum _{l=0}^{\mathbb{M}-m}\sum _{k=0}^{m+l}{a}_{l,m}\left(\begin{array}{c}m+l\\ k\end{array}\right)2^{\mathbb{k}k}{(-n)}^{m+l-k}\\ \hfill \times {\displaystyle \frac{\Gamma (1+\alpha k)}{\Gamma (1+\alpha k+\omega (\mathbf{t}\left)\right)}}{\mathbf{t}}^{\alpha k+\omega \left(\mathbf{t}\right)}\left[\mathbb{I}\left({\displaystyle \frac{c_2}{\mathbf{t}}},{\displaystyle \frac{c_1}{\mathbf{t}}};1+\alpha k,\omega \left(\mathbf{t}\right)\right)\right].\end{array}$$

The theorem is then proved. □

In the two-dimensional case, for $\alpha ,{\alpha }^{\prime }>0$, and non-negative integers $\mathbb{k},{\mathbb{k}}^{\prime },\mathbb{M},{\mathbb{M}}^{\prime }$ the set

$$\begin{array}{c}\mathcal{O}:=\left\{({\psi }_{m,n}^{\alpha }·{\psi }_{m^{\prime },n^{\prime }}^{{\alpha }^{\prime }})\right|0\le m\le 2^{\mathbb{k}}-1,0\le n\le \mathbb{M};0\le m^{\prime }\le 2^{{\mathbb{k}}^{\prime }}-1,0\le n^{\prime }\le {\mathbb{M}}^{\prime }\}\end{array}$$

forms an orthonormal set in $L_{w^{*}}^2([0,1]\times [0,1])$. Here, the norm ${\parallel .\parallel }_{w^{*}}$ is defined by

$${\parallel p(\mathbf{x},\mathbf{t})\parallel }_{w^{*}}^2={\int }_0^1{\int }_0^1{\left[p(\mathbf{x},\mathbf{t})\right]}^2{w}^{\alpha }\left(\mathbf{x}\right)w^{{\alpha }^{\prime }}\left(\mathbf{t}\right)d\mathbf{x}d\mathbf{t},$$

where $w^{\alpha }\left(\mathbf{x}\right)=\alpha {\mathbf{x}}^{\alpha -1}$ and $w^{{\alpha }^{\prime }}\left(\mathbf{x}\right)={\alpha }^{\prime }{\mathbf{x}}^{{\alpha }^{\prime -1}}$. To create a new arrangement, let ${\psi }_i^{\alpha }\left(\mathbf{x}\right)={\psi }_{m,n}^{\alpha }\left(\mathbf{x}\right),$ and ${\psi }_j^{{\alpha }^{\prime }}\left(\mathbf{t}\right)={\psi }_{m,n}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)$, where $i=n(\mathbb{M}+1)+m+1,\mathrm{and}j=n^{\prime }({\mathbb{M}}^{\prime }+1)+m^{\prime }+1$, then the orthonormal set $\mathcal{O}$ becomes

$$\{{\psi }_i^{\alpha }\left(\mathbf{x}\right)·{\psi }_j^{{\alpha }^{\prime }}\left(\mathbf{t}\right)|i=1,\dots \widehat{\mathbb{M}},j=1,\dots ,{\widehat{\mathbb{M}}}^{\prime }\},$$

where $\widehat{\mathbb{M}}=2^{\mathbb{k}}(\mathbb{M}+1)\mathrm{and}{\widehat{\mathbb{M}}}^{\prime }=2^{{\mathbb{k}}^{\prime }}({\mathbb{M}}^{\prime }+1)$. Any function $p(\mathbf{x},\mathbf{t})\in L_{w^{*}}^2\left({[0,1]}^2\right)$ can be represented by

$$\begin{array}{c}p(\mathbf{x},\mathbf{t})\approx \overline{p}(\mathbf{x},\mathbf{t})=\sum _{i=1}^{\widehat{\mathbb{M}}}\sum _{j=1}^{{\widehat{\mathbb{M}}}^{\prime }}{c}_{i,j}{\psi }_i^{\alpha }\left(\mathbf{x}\right){\psi }_j^{{\alpha }^{\prime }}\left(\mathbf{t}\right)={\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right),\end{array}$$

where ${\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(\mathbf{x}\right)$ and ${\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)$ are given as

$$\begin{array}{c}{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(\mathbf{x}\right)={\left[{\psi }_1^{\alpha }\left(\mathbf{x}\right),{\psi }_2^{\alpha }\left(\mathbf{x}\right),\dots ,{\psi }_{\widehat{\mathbb{M}}}^{\alpha }\left(\mathbf{x}\right)\right]}^T,\end{array}$$

(7)

$$\begin{array}{c}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)={\left[{\psi }_1^{{\alpha }^{\prime }}\left(\mathbf{t}\right),{\psi }_2^{{\alpha }^{\prime }}\left(\mathbf{t}\right),\dots ,{\psi }_{{\widehat{\mathbb{M}}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\right]}^T,\end{array}$$

(8)

and $\Delta =\left[u_{i,j}\right]$ is an $\widehat{\mathbb{M}}\times {\widehat{\mathbb{M}}}^{\prime }$ matrix with

$$\begin{array}{ccc}\hfill u_{i,j}& =& {〈{\psi }_i^{\alpha }\left(\mathbf{x}\right),{〈p(\mathbf{x},\mathbf{t}),{\psi }_j^{{\alpha }^{\prime }}\left(\mathbf{t}\right)〉}_{w^{{\alpha }^{\prime }}\left(\mathbf{t}\right)}〉}_{w^{\alpha }\left(\mathbf{x}\right)}\hfill \\ & =& {\int }_0^1{\int }_0^1{\psi }_i^{\alpha }\left(\mathbf{x}\right){\psi }_j^{{\alpha }^{\prime }}\left(\mathbf{t}\right)p(\mathbf{x},\mathbf{t}){\omega }^{\alpha }\left(\mathbf{x}\right){\omega }^{{\alpha }^{\prime }}\left(\mathbf{t}\right)d\mathbf{x}d\mathbf{t}.\hfill \end{array}$$

4. Error Bound

Theorem 2. 

Let $\Omega =[0,1]\times [0,1]$ and let $p:\Omega \to \mathbb{R}$ be a smooth function. Suppose that $\overline{p}$ is the best approximation of p out of

$$Span\{{\psi }_i^{\alpha }\left(\mathbf{x}\right){\psi }_j^{{\alpha }^{\prime }}\left(\mathbf{t}\right),1\le i\le \widehat{\mathbb{M}},1\le j\le {\widehat{\mathbb{M}}}^{\prime }\}.$$

with respect to the norm in $L_{w^{*}}^2(\Omega )$. Then

$$\begin{array}{ccc}\hfill \parallel p(\mathbf{x},\mathbf{t})-\overline{p}{(\mathbf{x},\mathbf{t})\parallel }_{w^{*}}^2& \le & {\displaystyle \frac{2^{1-(\mathbb{k}+2)(\mathbb{M}+1)}}{{\alpha }^{\mathbb{M}+1}·(\mathbb{M}+1)!}}·L_1+{\displaystyle \frac{2^{1-({\mathbb{k}}^{\prime }+2)({\mathbb{M}}^{\prime }+1)}}{{{\alpha }^{\prime }}^{{\mathbb{M}}^{\prime }+1}·({\mathbb{M}}^{\prime }+1)!}}·L_2\hfill \\ & +& {\displaystyle \frac{2^{2-(\mathbb{k}+2)(\mathbb{M}+1)-({\mathbb{k}}^{\prime }+2)({\mathbb{M}}^{\prime }+1)}}{{\alpha }^{\mathbb{M}+1}·{{\alpha }^{\prime }}^{{\mathbb{M}}^{\prime }+1}·(\mathbb{M}+1)!({\mathbb{M}}^{\prime }+1)!}}·L_3,\hfill \end{array}$$

where $L_1,L_2,\mathrm{and}{L}_2$ are constants satisfying

$$\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{\mathbb{M}+1}p(\mathbf{x},\mathbf{t})}{\partial {\mathbf{x}}^{\mathbb{M}+1}}}\right|\le L_1,\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{{\mathbb{M}}^{\prime }+1}p(\mathbf{x},\mathbf{t})}{\partial {\mathbf{t}}^{{\mathbb{M}}^{\prime }+1}}}\right|\le L_2,\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{\mathbb{M}+{\mathbb{M}}^{\prime }+2}p(\mathbf{x},\mathbf{t})}{\partial {\mathbf{x}}^{\mathbb{M}+1}\partial {\mathbf{t}}^{{\mathbb{M}}^{\prime }+1}}}\right|\le L_3.$$

Proof. 

We denote $I_{n,\mathbb{k}}^{\alpha }=\left[{\left({\displaystyle \frac{n}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}},{\left({\displaystyle \frac{n+1}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}\right)$ and $I_{n^{\prime },{\mathbb{k}}^{\prime }}^{{\alpha }^{\prime }}=\left[{\left({\displaystyle \frac{n^{\prime }}{2^{{\mathbb{k}}^{\prime }}}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}},{\left({\displaystyle \frac{n^{\prime }+1}{2^{{\mathbb{k}}^{\prime }}}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}}\right)$. Suppose that the polynomial $v(\mathbf{x},\mathbf{t})$ interpolates the function $p(\mathbf{x},\mathbf{t})$ at the points $(x_i,t_j)$, where $x_i$s and $t_{j}s$ are the roots of the Chebyshev polynomial of $(\mathbb{M}+1)$-degree in $I_{n,\mathbb{k}}^{\alpha }$ and $({\mathbb{M}}^{\prime }+1)$-degree in $I_{n^{\prime },{\mathbb{k}}^{\prime }}^{{\alpha }^{\prime }}$ correspondingly.

Thus

$$\begin{array}{c}\hfill \parallel p(\mathbf{x},\mathbf{t})-\overline{p}{(\mathbf{x},\mathbf{t})\parallel }_{\infty }\le {\parallel p(\mathbf{x},\mathbf{t})-v(\mathbf{x},\mathbf{t})\parallel }_{\infty }.\end{array}$$

(9)

Similar to the methods in References [36,42] and [43,44], we have

$$\begin{array}{ccc}\hfill p(\mathbf{x},\mathbf{t})-v(\mathbf{x},\mathbf{t})& =& {\displaystyle \frac{1}{(\mathbb{M}+1)!}}·{\displaystyle \frac{{\partial }^{\mathbb{M}+1}u(\eta ,\mathbf{t})}{\partial {\mathbf{x}}^{\mathbb{M}+1}}}·\prod _{i=0}^{\mathbb{M}}(\mathbf{x}-x_i)\hfill \\ \hfill & +& {\displaystyle \frac{1}{({\mathbb{M}}^{\prime }+1)!}}·{\displaystyle \frac{{\partial }^{{\mathbb{M}}^{\prime }+1}u(\mathbf{x},\xi )}{\partial {\mathbf{t}}^{{\mathbb{M}}^{\prime }+1}}}·\prod _{j=0}^{{\mathbb{M}}^{\prime }}(\mathbf{t}-t_j)\hfill \\ \hfill & -& {\displaystyle \frac{1}{(\mathbb{M}+1)!({\mathbb{M}}^{\prime }+1)!}}·{\displaystyle \frac{{\partial }^{\mathbb{M}+{\mathbb{M}}^{\prime }+2}u(\overline{\eta },\overline{\xi })}{\partial {\mathbf{x}}^{\mathbb{M}+1}\partial {\mathbf{t}}^{{\mathbb{M}}^{\prime }+1}}}·\prod _{i=0}^{\mathbb{M}}(\mathbf{x}-x_i)·\prod _{j=0}^{{\mathbb{M}}^{\prime }}(\mathbf{t}-t_j),\hfill \\ \hfill \end{array}$$

where $(\eta ,\xi ),(\overline{\eta },\overline{\xi })\in [0,1]\times [0,1]$. So we obtain

$$\begin{array}{ccc}\hfill {\parallel p(\mathbf{x},\mathbf{t})-v(\mathbf{x},\mathbf{t})\parallel }_{\infty }& \le & {\displaystyle \frac{1}{(\mathbb{M}+1)!}}·\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{\mathbb{M}+1}u(\eta ,\mathbf{t})}{\partial {\mathbf{x}}^{\mathbb{M}+1}}}\right|·\parallel \prod _{i=0}^{\mathbb{M}}(\mathbf{x}-x_i){\parallel }_{\infty }\hfill \\ \hfill & +& {\displaystyle \frac{1}{({\mathbb{M}}^{\prime }+1)!}}·\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{{\mathbb{M}}^{\prime }+1}u(\mathbf{x},\xi )}{\partial {\mathbf{t}}^{{\mathbb{M}}^{\prime }+1}}}\right|·\parallel \prod _{j=0}^{{\mathbb{M}}^{\prime }}(\mathbf{t}-t_j){\parallel }_{\infty }\hfill \\ \hfill & +& {\displaystyle \frac{1}{(\mathbb{M}+1)!({\mathbb{M}}^{\prime }+1)!}}·\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{\mathbb{M}+{\mathbb{M}}^{\prime }+2}u(\overline{\eta },\overline{\xi })}{\partial {\mathbf{x}}^{\mathbb{M}+1}\partial {\mathbf{t}}^{{\mathbb{M}}^{\prime }+1}}}\right|\hfill \\ \hfill & ·& \parallel \prod _{i=0}^{\mathbb{M}}(\mathbf{x}-x_i){\parallel }_{\infty }·\parallel \prod _{j=0}^{{\mathbb{M}}^{\prime }}(\mathbf{t}-t_j){\parallel }_{\infty }.\hfill \\ \hfill \end{array}$$

(10)

Since $p(\mathbf{x},\mathbf{t})$ is a smooth function on the compact domain $\Omega$, then there exist $L_1,L_2,\mathrm{and}{L}_3$ satisfying

$$\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{\mathbb{M}+1}p(\mathbf{x},\mathbf{t})}{\partial {\mathbf{x}}^{\mathbb{M}+1}}}\right|\le L_1,\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{{\mathbb{M}}^{\prime }+1}p(\mathbf{x},\mathbf{t})}{\partial {\mathbf{t}}^{{\mathbb{M}}^{\prime }+1}}}\right|\le L_2,\underset{(\mathbf{x},\mathbf{t})\in \Omega }{max}\left|{\displaystyle \frac{{\partial }^{\mathbb{M}+{\mathbb{M}}^{\prime }+2}p(\mathbf{x},\mathbf{t})}{\partial {\mathbf{x}}^{\mathbb{M}+1}\partial {\mathbf{t}}^{{\mathbb{M}}^{\prime }+1}}}\right|\le L_3.$$

(11)

To evaluate $\parallel {\prod }_{i=0}^{\mathbb{M}}(\mathbf{x}-x_i){\parallel }_{\infty }$, we consider the mapping

$$\mathbf{x}={\displaystyle \frac{a}{2^{\frac{\mathbb{k}}{\alpha }+1}}}·\mathbf{z}+{\displaystyle \frac{b}{2^{\frac{\mathbb{k}}{\alpha }+1}}},$$

where $a={(n+1)}^{{\displaystyle \frac{1}{\alpha }}}-n^{{\displaystyle \frac{1}{\alpha }}}$, and $b={(n+1)}^{{\displaystyle \frac{1}{\alpha }}}+n^{{\displaystyle \frac{1}{\alpha }}}$ between the intervals $[-1,1]$ and $\left[{\left({\displaystyle \frac{n}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}},{\left({\displaystyle \frac{n+1}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}\right)$, we get

$$\begin{array}{ccc}\hfill \underset{x_i\in [0,1]}{min}\underset{\mathbf{x}\in [0,1]}{max}\left|\prod _{i=0}^{\mathbb{M}}(\mathbf{x}-x_i)\right|& =& {\displaystyle \frac{a^{\mathbb{M}+1}}{2^{(\mathbb{M}+1)(\frac{\mathbb{k}}{\alpha }+1)}}}·\underset{z_i\in [-1,1]}{min}\underset{z\in [-1,1]}{max}\left|\prod _{i=0}^{\mathbb{M}}(\mathbf{z}-z_i)\right|\hfill \\ \hfill & \le & {\displaystyle \frac{a^{\mathbb{M}+1}}{2^{(\mathbb{M}+1)(\frac{\mathbb{k}}{\alpha }+2)-1}}},\hfill \end{array}$$

where $z_0,z_1,\dots ,z_{\mathbb{M}}$ are the solutions of the Chebyshev polynomial of degree $\mathbb{M}+1$.

Furthermore, by applying the mean value theorem,

$$\begin{array}{c}\hfill a={(n+1)}^{{\displaystyle \frac{1}{\alpha }}}-n^{{\displaystyle \frac{1}{\alpha }}}={\displaystyle \frac{1}{\alpha }}·{\zeta }^{{\displaystyle \frac{1}{\alpha }}-1},\zeta \in (n,n+1).\end{array}$$

We get

$$\begin{array}{c}\hfill a^{\mathbb{M}+1}={\displaystyle \frac{1}{{\alpha }^{\mathbb{M}+1}}}·{\zeta }^{(\mathbb{M}+1)({\displaystyle \frac{1}{\alpha }}-1)}\le {\displaystyle \frac{1}{{\alpha }^{\mathbb{M}+1}}}·2^{\mathbb{k}(\mathbb{M}+1)({\displaystyle \frac{1}{\alpha }}-1)}.\end{array}$$

So we obtain

$$\begin{array}{c}\hfill \underset{x_i\in [0,1]}{min}\underset{\mathbf{x}\in [0,1]}{max}\left|\prod _{i=0}^{\mathbb{M}}(\mathbf{x}-x_i)\right|\le {\displaystyle \frac{2^{1-(\mathbb{k}+2)(\mathbb{M}+1)}}{{\alpha }^{\mathbb{M}+1}}}.\end{array}$$

(12)

Similarly, we receive

$$\underset{t_j\in [0,1]}{min}\underset{\mathbf{t}\in [0,1]}{max}\left|\prod _{j=0}^{{\mathbb{M}}^{\prime }}(\mathbf{t}-t_j)\right|\le {\displaystyle \frac{2^{1-({\mathbb{k}}^{\prime }+2)({\mathbb{M}}^{\prime }+1)}}{{{\alpha }^{\prime }}^{{\mathbb{M}}^{\prime }+1}}}.$$

(13)

From Eqs. (10)-(13) we obtain

$$\begin{array}{c}\hfill {\parallel p(\mathbf{x},\mathbf{t})-v(\mathbf{x},\mathbf{t})\parallel }_{\infty }\le \mathbf{L},\end{array}$$

(14)

where

$$\begin{array}{ccc}\hfill \mathbf{L}& =& {\displaystyle \frac{2^{1-(\mathbb{k}+2)(\mathbb{M}+1)}}{{\alpha }^{\mathbb{M}+1}·(\mathbb{M}+1)!}}·L_1+{\displaystyle \frac{2^{1-({\mathbb{k}}^{\prime }+2)({\mathbb{M}}^{\prime }+1)}}{{{\alpha }^{\prime }}^{{\mathbb{M}}^{\prime }+1}·({\mathbb{M}}^{\prime }+1)!}}·L_2\hfill \\ & +& {\displaystyle \frac{2^{2-(\mathbb{k}+2)(\mathbb{M}+1)-({\mathbb{k}}^{\prime }+2)({\mathbb{M}}^{\prime }+1)}}{{\alpha }^{\mathbb{M}+1}·{{\alpha }^{\prime }}^{{\mathbb{M}}^{\prime }+1}·(\mathbb{M}+1)!({\mathbb{M}}^{\prime }+1)!}}·L_3.\hfill \end{array}$$

Note that $\Omega$ is divided to $2^{\mathbb{k}+{\mathbb{k}}^{\prime }}$ subdomains ${\Omega }^{\prime }=I_{n,\mathbb{k}}^{\alpha }\times I_{n^{\prime },{\mathbb{k}}^{\prime }}^{{\alpha }^{\prime }}$, we have

$$\begin{array}{ccc}\hfill \parallel p(\mathbf{x},\mathbf{t})-\overline{p}{(\mathbf{x},\mathbf{t})\parallel }_{w^{*}}^2& =& \sum _{n=0}^{2^{\mathbb{k}}-1}\sum _{n^{\prime }=0}^{2^{{\mathbb{k}}^{\prime }}-1}{\int }_{{\left({\displaystyle \frac{n}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}}^{{\left({\displaystyle \frac{n+1}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}}{\int }_{{\left({\displaystyle \frac{n^{\prime }}{2^{{\mathbb{k}}^{\prime }}}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}}}^{{\left({\displaystyle \frac{n^{\prime }+1}{2^{{\mathbb{k}}^{\prime }}}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}}}{w}^{\alpha }\left(\mathbf{x}\right)w^{{\alpha }^{\prime }}\left(\mathbf{t}\right){(p(\mathbf{x},\mathbf{t})-\overline{p}(\mathbf{x},\mathbf{t}))}^{2}d\mathbf{x}d\mathbf{t}.\hfill \end{array}$$

From Eqs. (9) and (14), we get

$$\begin{array}{ccc}\hfill \parallel p(\mathbf{x},\mathbf{t})-\overline{p}{(\mathbf{x},\mathbf{t})\parallel }_{w^{*}}^2& \le & \mathbf{L}·\sum _{n=0}^{2^{\mathbb{k}}-1}\sum _{n^{\prime }=0}^{2^{{\mathbb{k}}^{\prime }}-1}{\int }_{{\left({\displaystyle \frac{n}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}}^{{\left({\displaystyle \frac{n+1}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}}{\int }_{{\left({\displaystyle \frac{n^{\prime }}{2^{{\mathbb{k}}^{\prime }}}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}}}^{{\left({\displaystyle \frac{n^{\prime }+1}{2^{{\mathbb{k}}^{\prime }}}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}}}{w}^{\alpha }\left(\mathbf{x}\right)w^{{\alpha }^{\prime }}\left(\mathbf{t}\right)d\mathbf{x}d\mathbf{t}\hfill \\ \hfill & \le & \mathbf{L}·\sum _{n=0}^{2^{\mathbb{k}}-1}\sum _{n^{\prime }=0}^{2^{{\mathbb{k}}^{\prime }}-1}{\int }_{{\left({\displaystyle \frac{n}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}}^{{\left({\displaystyle \frac{b+1}{2^{\mathbb{k}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}}}{\int }_{{\left({\displaystyle \frac{n^{\prime }}{2^{{\mathbb{k}}^{\prime }}}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}}}^{{\left({\displaystyle \frac{b^{\prime }+1}{2^{{\mathbb{k}}^{\prime }}}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}}}\alpha {\mathbf{x}}^{\alpha -1}·{\alpha }^{\prime }{\mathbf{t}}^{{\alpha }^{\prime }-1}d\mathbf{x}d\mathbf{t}\hfill \\ \hfill & =& \mathbf{L}.\hfill \end{array}$$

The theorem is proved. □

5. Numerical Method

We use FOCW to establish a numerical method for solving the Eq. (1) associated with the conditions (2) and (3). First, we approximate

$$\begin{array}{c}{\displaystyle \frac{{\partial }^4}{\partial {\mathbf{x}}^2\partial {\mathbf{t}}^2}}p(\mathbf{x},\mathbf{t})\approx {\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T·\Delta ·{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right),\end{array}$$

where ${\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(\mathbf{x}\right),{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)$ are given in Eqs. (7) and (8). By integrating both sides of the above equation with respect to t twice and applying the Eq. (2)

$$\begin{array}{c}{\displaystyle \frac{{\partial }^2}{\partial {\mathbf{x}}^2}}p(\mathbf{x},\mathbf{t})\approx {\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T·\Delta ·I_{\mathbf{t}}^2{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)+{\displaystyle \frac{d^2}{d{\mathbf{x}}^2}}{r}_0\left(\mathbf{x}\right)+{\displaystyle \frac{d^2}{d{\mathbf{x}}^2}}{r}_1\left(\mathbf{x}\right)\mathbf{t}.\end{array}$$

(15)

Integrating both sides of Eq. (15) with respect to $\mathbf{x}$ twice and using the conditions (3) we get

$$\begin{array}{ccc}\hfill p(\mathbf{x},\mathbf{t})& \approx & I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .I_{\mathbf{t}}^2{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)-\mathbf{x}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(1\right)}^T.\Delta .I_{\mathbf{t}}^2{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\hfill \\ & +& (1-\mathbf{x})\left[s_0\left(\mathbf{t}\right)+(r_0\left(\mathbf{x}\right)-r_0\left(0\right))+(r_1\left(\mathbf{x}\right)-r_1\left(0\right))\mathbf{t}\right]+\mathbf{x}{s}_1\left(\mathbf{t}\right),\hfill \end{array}$$

(16)

where

$$\begin{array}{c}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(\mathbf{x}\right)={\left[I_{\mathbf{x}}^2{\psi }_1^{\alpha }\left(\mathbf{x}\right),I_{\mathbf{x}}^2{\psi }_2^{\alpha }\left(\mathbf{x}\right),\dots ,I_{\mathbf{x}}^2{\psi }_{\mathbf{M}}^{\alpha }\left(\mathbf{x}\right)\right]}^T,\end{array}$$

and

$$\begin{array}{c}{I}_{\mathbf{t}}^2{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)={\left[I_{\mathbf{t}}^2{\psi }_1^{{\alpha }^{\prime }}\left(\mathbf{t}\right),I_{\mathbf{t}}^2{\psi }_2^{{\alpha }^{\prime }}\left(\mathbf{t}\right),\dots ,I_{\mathbf{t}}^2{\psi }_{{\mathbf{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\right]}^T.\end{array}$$

Therefore

$$\begin{array}{ccc}\hfill D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})& \approx & I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)-\mathbf{x}·I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(1\right)}^T.\Delta .I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\hfill \\ & +& (1-\mathbf{x})\left[D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}{s}_0\left(\mathbf{t}\right)+{\displaystyle \frac{r_1\left(\mathbf{t}\right)-r_1\left(0\right)}{\Gamma (2-\omega (\mathbf{x},\mathbf{t})}}{\mathbf{t}}^{1-\omega (\mathbf{x},\mathbf{t})}\right]+\mathbf{x}{D}_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}{s}_1\left(\mathbf{t}\right),\hfill \end{array}$$

(17)

and

$$\begin{array}{ccc}\hfill D_{\mathbf{t}}^{{\omega }_i(\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})& \approx & I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .I_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\hfill \\ & -& \mathbf{x}·I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(1\right)}^T.\Delta .I_t^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\hfill \\ & +& (1-\mathbf{x})\left[D_{\mathbf{t}}^{{\omega }_i(\mathbf{x},\mathbf{t})}{s}_0\left(\mathbf{t}\right)+{\displaystyle \frac{r_1\left(\mathbf{x}\right)-r_1\left(0\right)}{\Gamma (2-{\omega }_i(\mathbf{x},\mathbf{t})}}{\mathbf{t}}^{1-{\omega }_i(\mathbf{x},\mathbf{t})}\right]+\mathbf{x}{D}_{\mathbf{t}}^{{\omega }_i(\mathbf{x},\mathbf{t})}{s}_1\left(\mathbf{t}\right),\hfill \end{array}$$

(18)

where

$$\begin{array}{c}{I}_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)={\left[I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\psi }_1^{{\alpha }^{\prime }}\left(\mathbf{t}\right),I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\psi }_2^{{\alpha }^{\prime }}\left(\mathbf{t}\right),\dots ,I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\psi }_{{\mathbf{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\right]}^T,\end{array}$$

$$\begin{array}{c}{I}_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)={\left[I_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\psi }_1^{{\alpha }^{\prime }}\left(\mathbf{t}\right),I_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\psi }_2^{{\alpha }^{\prime }}\left(\mathbf{t}\right),\dots ,I_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\psi }_{{\mathbf{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\right]}^T,\end{array}$$

and $D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}{s}_i\left(\mathbf{t}\right),D_{\mathbf{t}}^{{\omega }_i(\mathbf{x},\mathbf{t})}{s}_i\left(\mathbf{t}\right);i=0,1,$ are given in Eq. (5). From Eq. (16) we obtain

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial p(\mathbf{x},\mathbf{t})}{\partial \mathbf{t}}}& =& I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T·\Delta ·I_{\mathbf{t}}^1{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)-\mathbf{x}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(1\right)}^T·\Delta ·I_{\mathbf{t}}^1{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\hfill \\ & +& (1-\mathbf{x})[s_0^{\prime }\left(\mathbf{t}\right)+r_1\left(\mathbf{x}\right)-r_1\left(0\right)]+\mathbf{x}{s}_1^{\prime }\left(\mathbf{t}\right).\hfill \end{array}$$

(19)

By substituting Eqs. (15)-(19) in Eq. (1), we get an algebraic equation. By collocating at $\widehat{\mathbb{M}}·{\widehat{\mathbb{M}}}^{\prime }$ nodes given by

$$(x_i,t_j)=\left({\left({\displaystyle \frac{2i-1}{2\widehat{\mathbb{M}}}}\right)}^{{\displaystyle \frac{1}{\alpha }}},{\left({\displaystyle \frac{2j-1}{2{\widehat{\mathbb{M}}}^{\prime }}}\right)}^{{\displaystyle \frac{1}{{\alpha }^{\prime }}}}\right),i=1,\dots ,\widehat{\mathbb{M}};j=1,\dots ,{\widehat{\mathbb{M}}}^{\prime },$$

we can find the $\widehat{\mathbb{M}}·{\widehat{\mathbb{M}}}^{\prime }$ unknowns $u_{i,j}$.

6. Illustrative Examples

We compare the experimental results obtained by using our method and with known methods, such as shifted Chebyshev cardinal functions [42], Chebyshev wavelets [36], and generalized polynomials [45] in the following examples:

Example 1 

(See [36, Example 1]). Consider the following FO-FDWE:

$$D_t^{\omega (\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+D_t^{{\omega }_1(\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+p_{\mathbf{t}}(\mathbf{x},\mathbf{t})=p_{\mathbf{xx}}(\mathbf{x},\mathbf{t})+q(\mathbf{x},\mathbf{t}),$$

(20)

where $\omega (\mathbf{x},\mathbf{t})=2-0.3e^{-\mathbf{xt}},{\omega }_1(\mathbf{x},\mathbf{t})=2-0.6e^{-\mathbf{xt}}$,

$$\begin{array}{ccc}\hfill q(\mathbf{x},\mathbf{t})& =& 100{\mathbf{x}}^{3-a}{\mathbf{t}}^3(1-\mathbf{x})\left({\displaystyle \frac{6}{\Gamma (4-\omega (\mathbf{x},\mathbf{t}\left)\right)}}{\mathbf{t}}^{-\omega (\mathbf{x},\mathbf{t})}-{\displaystyle \frac{24}{\Gamma (5-\omega (\mathbf{x},\mathbf{t}\left)\right)}}{\mathbf{t}}^{1-\omega (\mathbf{x},\mathbf{t})}\right)\hfill \\ & +& 100{\mathbf{x}}^{3-a}{\mathbf{t}}^3(1-\mathbf{x})\left({\displaystyle \frac{6}{\Gamma (4-{\omega }_1(\mathbf{x},\mathbf{t}))}}{\mathbf{t}}^{-{\omega }_1(\mathbf{x},\mathbf{t})}-{\displaystyle \frac{24}{\Gamma (5-{\omega }_1(\mathbf{x},\mathbf{t}))}}{\mathbf{t}}^{1-{\omega }_1(\mathbf{x},\mathbf{t})}\right)\hfill \\ & +& 100{\mathbf{x}}^{3-a}(1-\mathbf{x})(3{\mathbf{t}}^2-4{\mathbf{t}}^3)-100{\mathbf{t}}^3(1-\mathbf{t})((3-a)(2-1){\mathbf{x}}^{1-a}-(4-a)(3-a){\mathbf{x}}^{2-a}),\hfill \end{array}$$

and the conditions in Eqs. (2) and (3) are given so that

$$p(\mathbf{x},\mathbf{t})=100{\mathbf{x}}^{3-a}{\mathbf{t}}^3(1-\mathbf{x})(1-\mathbf{t}).$$

This example was considered in [36] with $a=0$. Here, we applied our method with $\mathbb{k}={\mathbb{k}}^{\prime }=1,\mathbb{M}={\mathbb{M}}^{\prime }=2,\mathrm{and}\alpha ={\alpha }^{\prime }=1$. First, we evaluate

$$\begin{array}{c}{\displaystyle \frac{{\partial }^4}{\partial {\mathbf{x}}^2\partial {\mathbf{t}}^2}}p(\mathbf{x},\mathbf{t})\approx {\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)\end{array}$$

where

$$\begin{array}{c}{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(\mathbf{x}\right)={\left[{\psi }_1^1\left(\mathbf{x}\right),{\psi }_2^1\left(\mathbf{x}\right),{\psi }_3^1\left(\mathbf{x}\right),{\psi }_4^1\left(\mathbf{x}\right),{\psi }_5^1\left(\mathbf{x}\right),{\psi }_6^1\left(\mathbf{x}\right)\right]}^T,\end{array}$$

$$\begin{array}{c}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)={\left[{\psi }_1^1\left(\mathbf{t}\right),{\psi }_2^1\left(\mathbf{t}\right),{\psi }_3^1\left(\mathbf{t}\right),{\psi }_4^1\left(\mathbf{t}\right),{\psi }_5^1\left(\mathbf{t}\right),{\psi }_6^1\left(\mathbf{t}\right)\right]}^T,\end{array}$$

and

$$\Delta ={\left(u_{i,j}\right)}_{i,j=1,...,6}$$

From Eqs. (15) and (16) we obtain the following.

$${\displaystyle \frac{{\partial }^2}{\partial {\mathbf{x}}^2}}p(\mathbf{x},\mathbf{t})\approx {\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .I_{\mathbf{t}}^2{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right),$$

(21)

$$p(\mathbf{x},\mathbf{t})=I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .I_{\mathbf{t}}^2{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)-\mathbf{x}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(1\right)}^T.\Delta .I_{\mathbf{t}}^2{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right),$$

(22)

where

$$\begin{array}{c}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(\mathbf{x}\right)={\left[I_{\mathbf{x}}^2{\psi }_1^1\left(\mathbf{x}\right),I_{\mathbf{x}}^2{\psi }_2^1\left(\mathbf{x}\right),I_{\mathbf{x}}^2{\psi }_3^1\left(\mathbf{x}\right),I_{\mathbf{x}}^2{\psi }_4^1\left(\mathbf{x}\right),I_{\mathbf{x}}^2{\psi }_5^1\left(\mathbf{x}\right),I_{\mathbf{x}}^2{\psi }_6^1\left(\mathbf{x}\right)\right]}^T,\end{array}$$

$$\begin{array}{c}{I}_{\mathbf{t}}^2{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)={\left[I_{\mathbf{t}}^2{\psi }_1^1\left(\mathbf{t}\right),I_{\mathbf{t}}^2{\psi }_2^1\left(\mathbf{t}\right),I_{\mathbf{t}}^2{\psi }_3^1\left(\mathbf{t}\right),I_{\mathbf{t}}^2{\psi }_4^1\left(\mathbf{t}\right),I_{\mathbf{t}}^2{\psi }_5^1\left(\mathbf{t}\right),I_{\mathbf{t}}^2{\psi }_6^1\left(\mathbf{t}\right)\right]}^T,\end{array}$$

$$\begin{array}{c}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }\left(1\right)={\left[I_{\mathbf{x}}^2{\psi }_1^1\left(1\right),I_{\mathbf{x}}^2{\psi }_2^1\left(1\right),I_{\mathbf{x}}^2{\psi }_3^1\left(1\right),I_{\mathbf{x}}^2{\psi }_4^1\left(1\right),I_{\mathbf{x}}^2{\psi }_5^1\left(1\right),I_{\mathbf{x}}^2{\psi }_6^1\left(1\right)\right]}^T.\end{array}$$

By using Eqs. (17) and (18) we get

$$\begin{array}{c}{D}_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})\approx I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)-\mathbf{x}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(1\right)}^T.\Delta .I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right),\end{array}$$

(23)

and

$$\begin{array}{c}{D}_{\mathbf{t}}^{{\omega }_1(\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})\approx I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T.\Delta .I_{\mathbf{t}}^{2-{\omega }_1(\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)-\mathbf{x}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(1\right)}^T.\Delta .I_{\mathbf{t}}^{2-{\omega }_1(\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right),\end{array}$$

(24)

where

$$\begin{array}{ccc}\hfill I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)& =& {\left[I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\psi }_0^1\left(\mathbf{t}\right),I_{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\psi }_1^1\left(\mathbf{t}\right),\dots ,I{\mathbf{t}}^{2-\omega (\mathbf{x},\mathbf{t})}{\psi }_6^1\left(\mathbf{t}\right)\right]}^T,\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill I_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)& =& {\left[I_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\psi }_0^1\left(\mathbf{t}\right),I_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\psi }_1^1\left(\mathbf{t}\right)\dots ,I_{\mathbf{t}}^{2-{\omega }_i(\mathbf{x},\mathbf{t})}{\psi }_6^1\left(\mathbf{t}\right)\right]}^T.\hfill \end{array}$$

From Eq. (19) we get

$$p_{\mathbf{t}}(\mathbf{x},\mathbf{t})\approx I_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(\mathbf{x}\right)}^T·\Delta ·I_{\mathbf{t}}^1{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)-\mathbf{x}{I}_{\mathbf{x}}^2{\Psi }_{\mathbb{k},\mathbb{M}}^{\alpha }{\left(1\right)}^T·\Delta ·I_{\mathbf{t}}^1{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right),$$

(25)

where

$$\begin{array}{ccc}\hfill I_{\mathbf{t}}^1{\Psi }_{{\mathbb{k}}^{\prime },{\mathbb{M}}^{\prime }}^{{\alpha }^{\prime }}\left(\mathbf{t}\right)& =& {\left[I_{\mathbf{t}}^1{\psi }_1^1\left(\mathbf{t}\right),I_{\mathbf{t}}^1{\psi }_2^1\left(\mathbf{t}\right),I_{\mathbf{t}}^1{\psi }_3^1\left(\mathbf{t}\right),I_{\mathbf{t}}^1{\psi }_4^1\left(\mathbf{t}\right),I_{\mathbf{t}}^1{\psi }_5^1\left(\mathbf{t}\right),I_{\mathbf{t}}^1{\psi }_6^1\left(\mathbf{t}\right)\right]}^T.\hfill \end{array}$$

By substituting Eqs. (21) and (23)-(25) in Eq. (20) and collocating the resulted equation at the nodes

$$\begin{array}{c}\left(x_i,y_j\right)=\left({\displaystyle \frac{2i-1}{12}},{\displaystyle \frac{2j-1}{12}}\right),i,j=\overline{1,6}.\end{array}$$

we obtain a system of 36 algebraic equations in 36 unknowns $u_{i,j}$. The solution of the system is

$$\begin{array}{c}\Delta ={\displaystyle \frac{1}{8}}\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 75& 15\sqrt{15}& 0& -75& -135\sqrt{15}\\ 0& 15\sqrt{15}& 45& 0& 15\sqrt{15}& -405\\ 0& 0& 0& 0& 0& 0\\ 0& -75& 15\sqrt{15}& 0& 75& -135\sqrt{15}\\ 0& -135\sqrt{15}& -405& 0& -135\sqrt{15}& 3645\end{array}\right).\end{array}$$

By substituting this solution into Eq. (22), we get the exact solution $p(\mathbf{x},\mathbf{t})=100{\mathbf{x}}^3{\mathbf{t}}^3(1-\mathbf{x})(1-\mathbf{t})$. This solution was not obtained in [36].

In case $a=0.5$, by using the same argument as above with $\mathbb{k}={\mathbb{k}}^{\prime }=0,\mathbb{M}=3,{\mathbb{M}}^{\prime }=2,\alpha =0.5\mathrm{and}{\alpha }^{\prime }=1$, we obtain

$$\begin{array}{c}\hfill \Delta =\left(\begin{array}{ccc}0& 0& 0\\ 0& {\displaystyle \frac{75}{4}}& -{\displaystyle \frac{45\sqrt{15}}{4}}\\ 0& {\displaystyle \frac{75\sqrt{15}}{4}}& -{\displaystyle \frac{675}{4}}\\ 0& -25\sqrt{21}& 45\sqrt{35}\end{array}\right)\end{array}$$

which yields the exact solution $p(\mathbf{x},\mathbf{t})=100{\mathbf{x}}^{2.5}{\mathbf{t}}^3(1-\mathbf{x})(1-\mathbf{t})$. This exact solution cannot be obtained by using integer-order polynomial wavelets.

Example 2 

(See [36, Example 2]). Consider

$$D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+a_1{D}_{\mathbf{t}}^{{\omega }_1(\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+a_2{D}_{\mathbf{t}}^{{\omega }_2(\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+p_{\mathbf{t}}(\mathbf{x},\mathbf{t})=b{p}_{\mathbf{xx}}(\mathbf{x},\mathbf{t})+q(\mathbf{x},\mathbf{t}),(\mathbf{x},\mathbf{t})\in {[0,1]}^2,$$

where

$$q(\mathbf{x},\mathbf{t})=\left[{\displaystyle \frac{6{\mathbf{t}}^{-\omega (\mathbf{x},\mathbf{t})}}{\Gamma \left[4-\omega (\mathbf{x},\mathbf{t})\right]}}+{\displaystyle \frac{6a_1{\mathbf{t}}^{-{\omega }_1(\mathbf{x},\mathbf{t})}}{\Gamma \left[4-{\omega }_1(\mathbf{x},\mathbf{t})\right]}}+{\displaystyle \frac{6a_2{\mathbf{t}}^{-{\omega }_2(\mathbf{x},\mathbf{t})}}{\Gamma \left[4-{\omega }_2(\mathbf{x},\mathbf{t})\right]}}\right]{\mathbf{t}}^3{e}^{\mathbf{x}}+\left(3-b\mathbf{t}\right){\mathbf{t}}^2{e}^{\mathbf{x}},$$

$\omega (\mathbf{x},\mathbf{t})=0.3sin\left(\mathbf{xt}\right)+1.7$, ${\omega }_1(\mathbf{x},\mathbf{t})=0.3sin\left(\mathbf{xt}\right)+1.5$, ${\omega }_2(\mathbf{x},\mathbf{t})=0.3sin\left(\mathbf{xt}\right)+1.3$, $a_1={\displaystyle \frac{1}{2}},a_2={\displaystyle \frac{1}{3}}$ and $b={\displaystyle \frac{1}{2}}$. The conditions in Eqs. (2) and (3) are given such that $p(\mathbf{x},\mathbf{t})={\mathbf{t}}^3{e}^{\mathbf{x}}$.

First, we solve this problem with $\mathbb{k}={\mathbb{k}}^{\prime }=1,\mathbb{M}={\mathbb{M}}^{\prime }=3$. Table 1 presents the maximum of the relative errors (REs) in various pairs of fractional orders $(\alpha ,{\alpha }^{\prime })$, and the minimal of these errors is reached at $(\alpha ,{\alpha }^{\prime })=(1,{\displaystyle \frac{1}{2}})$. In Table 2, we compare the REs for this method with those in [36]. The comparison is based on the same number of elements in the bases of the function approximation. The table suggests that the present method provides approximate solutions with lower error. In Table 2, M is the degree of Chebyshev polynomials. Figure 1 illustrates the graphs of the approximated solution and the absolute error function (AEF) for $\mathbb{k}={\mathbb{k}}^{\prime }=1,\mathbb{M}={\mathbb{M}}^{\prime }=9,\mathrm{and}\alpha =1,{\alpha }^{\prime }={\displaystyle \frac{1}{2}}$.

Example 3 

(See [36, Example 3]). Consider the nonlinear FO-FDWE

$$D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+\sum _{i=1}^4{a}_i{D}_{\mathbf{t}}^{{\omega }_i(\mathbf{x},\mathbf{t})}p(\mathbf{x},\mathbf{t})+u_{\mathbf{t}}(\mathbf{x},\mathbf{t})=b{p}_{\mathbf{xx}}(\mathbf{x},\mathbf{t})+q(\mathbf{x},\mathbf{t}),(\mathbf{x},\mathbf{t})\in {[0,1]}^2,$$

where

$$\begin{array}{c}q(\mathbf{x},\mathbf{t})=\left({\displaystyle \frac{6{\mathbf{t}}^{-\omega (\mathbf{x},\mathbf{t})}}{\Gamma (4-\omega (\mathbf{x},\mathbf{t}\left)\right)}}+\sum _{i=1}^4{a}_i{\displaystyle \frac{6{\mathbf{t}}^{-{\omega }_i(\mathbf{x},\mathbf{t})}}{\Gamma (4-{\omega }_i(\mathbf{x},\mathbf{t}))}}\right){\mathbf{t}}^{3}sin\left(\pi \mathbf{x}\right)+(3+a{\pi }^2\mathbf{t}){\mathbf{t}}^{2}sin\left(\pi \mathbf{x}\right),\end{array}$$

$\omega (\mathbf{x},\mathbf{t})=2-0.3e^{-\mathbf{xt}};{\omega }_i(\mathbf{x},\mathbf{t})=2-{\displaystyle \frac{i+3}{10}}{e}^{-\mathbf{xt}},a_i={\displaystyle \frac{1}{i+1}},i=1,\dots ,4;$ and $b={\displaystyle \frac{1}{3}}$. The conditions in Eqs. (2) and (3) are given such that $p(\mathbf{x},\mathbf{t})={\mathbf{t}}^{3}sin\left(\pi \mathbf{x}\right)$.

Table 3 provides the maximal AEs of the approximated solutions using our method in the case $\mathbb{k}={\mathbb{k}}^{\prime }=1,\mathbb{M}={\mathbb{M}}^{\prime }=3$ and various pairs of fractional orders $(\alpha ,{\alpha }^{\prime })$. The table shows that the minimum occurred at $(\alpha ,{\alpha }^{\prime })=(1.25,0.5)$. Figure 2 shows the graph of the solution and AEF in cases $\mathbb{k}={\mathbb{k}}^{\prime }=1,\mathbb{M}={\mathbb{M}}^{\prime }=8$ and $\alpha =1.25,{\alpha }^{\prime }=0.5$. Table 4 and Table 5 compare the REs of numerical solutions using our method with the method in [36] with different values of $\mathbb{M}={\mathbb{M}}^{\prime }$. These tables indicate that the present method yields more accurate solutions with fewer errors.

Example 4. 

In particle physics, anomalous diffusion processes were formulated in terms of diffusion equations by using fractional derivatives and distributed-order derivatives in [46]. As a generalization, we consider the following diffusion equation with variable-order derivative:

$$\begin{array}{c}{D}_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}Q(\mathbf{x},\mathbf{t})=\nu {\displaystyle \frac{{\partial }^{2}Q(\mathbf{x},\mathbf{t})}{\partial {\mathbf{x}}^2}}\end{array}$$

(26)

where $(\mathbf{x},\mathbf{t})\in (-\infty ,+\infty )\times (0,+\infty )$, $Q(\mathbf{x},\mathbf{t})$ is the pdf (probability density function) of the diffusion-process, and ν is the diffusion-coefficient of physical dimension (cm²/sec). Following [46], we consider this diffusion equation:

$$\begin{array}{cc}\hfill & Q(\mathbf{x},0)=\underset{\mathbf{t}\to 0}{lim}Q(\mathbf{x},\mathbf{t})=1,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & Q(\pm \infty ,\mathbf{t})=0.\hfill \end{array}$$

(28)

We scale the problem to finite intervals by setting:

$$\begin{array}{c}\mathbf{x}=ln\left({\displaystyle \frac{\mathbf{z}}{1-\mathbf{z}}}\right),\mathrm{or}\mathrm{equivalently}\mathbf{z}={\displaystyle \frac{1}{1+e^{-\mathbf{x}}}}.\end{array}$$

(29)

Then $\mathbf{z}\in (0,1)$ and Eqs. (26), (27), (28) become

$$\begin{array}{cc}\hfill & D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}Q(\mathbf{z},\mathbf{t})=\nu {\displaystyle \frac{\partial Q(\mathbf{z},\mathbf{t})}{\partial \mathbf{x}}}\mathbf{z}(1-\mathbf{z})(1-2\mathbf{z})+\nu {\displaystyle \frac{{\partial }^{2}Q(\mathbf{z},\mathbf{t})}{\partial {\mathbf{x}}^2}}{\mathbf{z}}^2{(1-\mathbf{z})}^2\hfill \\ \hfill & Q(1,\mathbf{t})=0,Q(0,\mathbf{t})=0,Q(\mathbf{z},0)=\underset{\mathbf{t}\to 0}{lim}Q(\mathbf{z},\mathbf{t})=1,\hfill \end{array}$$

(30)

where

$$\begin{array}{c}p(\mathbf{z},\mathbf{t})=Q(\mathbf{z},\mathbf{t})=Q\left(ln(\mathbf{z}/(1-\mathbf{z})),\mathbf{t}\right).\end{array}$$

(31)

Following [46], we consider this equation with $\nu ={\displaystyle \frac{1}{2}}$. In addition, we use $k_1=k_2=1$, $({\alpha }_1,{\alpha }_2)=(1,0.5)$ and different values of $\mathbb{M}={\mathbb{M}}^{\prime }$. The numerical solution is illustrated in Figure 3. Table 6 reports the maximum residual errors along with the convergence order. Here, for each numerical solution $\tilde{u}$, the maximal residual error is

$$\mathrm{MaxRE}=\underset{(\mathbf{z},\mathbf{t})\in {[0,1]}^2}{max}\left|D_{\mathbf{t}}^{\omega (\mathbf{x},\mathbf{t})}\tilde{p}(\mathbf{z},\mathbf{t})-\nu {\displaystyle \frac{\partial \tilde{p}(\mathbf{z},\mathbf{t})}{\partial \mathbf{x}}}\mathbf{z}(1-\mathbf{z})(1-2\mathbf{z})-\nu {\displaystyle \frac{{\partial }^2\tilde{p}(\mathbf{z},\mathbf{t})}{\partial {\mathbf{z}}^2}}{\mathbf{z}}^2{(1-\mathbf{z})}^2\right|.$$

7. Conclusions

This study introduced a spectral wavelet-based approach utilizing generalized fractional-order Chelyshkov wavelets to solve functional-order fractional diffusion–wave equations. The method demonstrated high computational efficiency and superior accuracy, achieving significantly lower numerical errors than conventional schemes in benchmark comparisons. The integration of regularized incomplete Beta functions allowed for an exact representation of the fractional integral operator, contributing to the stability and effectiveness of the proposed framework. Notably, the method successfully captured the physical behavior in Example 4, illustrating its capability to handle realworld models with high fidelity. Given its flexibility and precision, the proposed strategy offers a promising foundation for solving a broader class of fractional partial differential equations with variable orders. Future work may extend this framework to diverse applications in physics and engineering where memory effects and anomalous diffusion are prominent.