The Novel Analytical–Numerical Method for Multi-Dimensional Multi-Term Time-Fractional Equations with General Boundary Conditions

Abstract

This article presents a simple but effective two-step analytical–numerical algorithm for solving multi-dimensional multi-term time-fractional equations. The first step is to derive an analytic representation that satisfies boundary requirements for 1D, 2D, and 3D problems, respectively. The second step is the meshless approximation where the Müntz polynomials are used to form the approximate solution and the unknown parameters are obtained by imposing the approximation for the governing equations. We illustrate first the detailed derivation of the analytic approximation and then the numerical implementation of the solution procedure. Several numerical examples are provided to verify the accuracy, efficiency, and adaptability to problems with general boundary conditions. The numerical results are compared with exact solutions and numerical methods reported in the literature, showing that the algorithm has great potential for multi-dimensional multi-term time-fractional equations with various boundary conditions.

1. Introduction

In recent years, many mathematical models with time-fractional multi-term derivatives have been studied in physics, hydrology, chemistry, etc. Hilfer [1] provided a collection of fractional calculus applications in physics. Molz [2] reviewed some fundamental properties of fractional motion and applications in hydrology. Singh [3] analyzed a chemical kinetics system pertaining to a fractional derivative. Furthermore, the applications of fractional equations can also be found in [4]. The well-known cable equations with fractional-order temporal operators belong to this class. They were introduced to model electronic properties of spiny neuronal dendrites [5,6]. The time-fractional partial differential equations (TFPDEs) can also include the Sobolev equations, which have been used to model many phenomena, such as the migration of moisture in soil, thermodynamics, and the motion of fluid in various media [7,8]. The equations of different models of heat transfer, whether the classical or dual-phase-lagging ones, undoubtedly also fall into this group [9]. Fractional differential equations have broad applications for the fact that fractional operators can describe physical phenomena more precisely than classical integral operators for some practical problems. The collection of real world applications of fractional differential equations can be seen in [10] and references therein.

Solutions to fractional partial differential equations are crucial for representing physical phenomena. Some analytical methods have been proposed that can be useful for parametric study. Liu carried out an analytical study for magnetohydrodynamic flow using fractional derivatives [11]. Ming derived the analytical solution for the problems containing multi-term time-fractional diffusion [12]. Ding proposed an analytical solution to the multi-term TFPDEs considering the non-local damping term [13] and also the fractional delay PDEs with mixed boundary conditions [14]. Jong used the analytic expression of the multi-term fractional integral operators to obtain the analytical expressions for the fractional equations [15]. Jiang investigated the multi-term fractional diffusion equations and obtained the analytical solutions by the method of separation variables [16]. In [17], the Laplace transform method was used to derive the solution to the time-fractional distributed-order heat conduction law. The Adomian decomposition method was also being used in [18]. Despite the fact that so many analytical techniques have been developed, the explicit forms of analytical or semi-analytical solutions are rare only for some problems under certain idealized conditions. For the further prospects of engineering applications, the numerical methods are still necessary and useful tools in this field. Numerical methods have already been used to observe some important mathematical models. Liu investigated the fluid mechanics for semiconductor circuit breakers based on finite element analysis [19]. Yang, Liu, and Xu applied functional differential equations to analyze the problems in financial accounting [20,21,22]. A computational heuristic was designed to solve the nonlinear Lienard differential model [23]. The nonuniform difference scheme was applied to study the distributed-order fractional parabolic equations with fractional Laplacian in [24]. A fast Fourier spectral exponential time-differencing method was used to solve time-fractional mobile–immobile equations [25]. A fast difference scheme was proposed to solve the fractional equations considering non-smooth data [26]. Some other works can be found in [27,28,29] and references therein. Among them, the most popular methods are mesh-based methods. Some of the references and studies are listed below. Dehghan [30] proposed a high-order numerical algorithm based on the finite difference scheme for multi-term time-fractional diffusion wave problems. The Galerkin finite element method [31] was proposed for the approximation of the multi-term time-fractional diffusion equations. The finite difference and the finite element method were used to solve the multi-term time-fractional equations that are mixed by the sub-diffusion and diffusion-wave equation [32]. The fast algorithm combined with the finite difference method was used to solve multi-term time-fractional reaction–diffusion wave equations with stability analysis and error analysis [33]. The second-order numerical method was proposed for the problems with non-smooth solutions [34]. As implied by the name, mesh-based methods require the mesh of the whole domain and also information about the nodal topology that may introduce some unreasonable constraints on the problems. The automatic and efficient approach to constructing mesh for 3D complicated domains has long been the challenge for computational mechanics. Spectral methods [35,36] and spectral-based methods, especially those based on the collocation method, have been published recently [37,38]. The spectral-based method has also been used to solve the distributed order time-fractional diffusion equations in [39]. A pseudo-spectral method based on the reproducing kernel has been proposed to study the time-fractional diffusion-wave equation [40].

In this paper, considering the advantages of analytical and numerical methods, a novel analytical–numerical method is proposed for solving multi-term TFPDEs with boundary conditions of general kinds. First, we apply the Fourier method, which can also be regarded as an expansion method over the eigenfunctions, in order to remove the partial derivatives with respect to the space variables and transform the original TFPDE into fractional ordinary differential equations (FODEs) without truncation error. In the general case of the time-dependent non-homogeneous boundary conditions, the solutions with the features of separation of spatial variables are not available naturally. The time-dependent non-homogeneous term in the equation also poses a problem for the application of the Fourier method. To deal with the time dependence boundary conditions and the source term, the Green function method and the operational methods such as the Laplace transform method are used [41]. A similar technique has been proposed for the time-fractional PDEs [42,43]. Then, we apply the recently proposed meshless collocation method, the backward substitution method (BSM) [44,45,46], to solve each FODE individually. In [44], the BSM was also applied to solve systems of FODEs. The BSM technique can also deal with the non-homogeneous time-dependent source term. The BSM is a newly developed meshless method. By introducing special analytical functions or numerical approximations that satisfy the boundary conditions, the original problem is degenerated into a homogeneous one. Then, the BSM attempts to form an orthogonal basis system that satisfies the homogeneous boundary conditions in a general way. The approximated solution is formed using the proposed basis system where the weighted parameters are determined by backward substituting the approximation into the governing equations. This improvement has significantly increased the accuracy and stability of the usual collocation methods. Recently, some variants of the BSM have been proposed, such as the space-time BSM [47], the localized BSM [48], and the Fourier-based BSM [49]. In order to apply the BSM for the fractional differential equations, some special bases should be used. Firstly, the solution of the fractional equations can contain the fractional-power terms where the common bases cannot be used for such purpose. Secondly, in terms of the BSM, the collocation method is applied. In order to apply the collocation method, it is critical that derivatives of the trial functions should be approximated by the same trial bases where the common polynomial bases cannot be used. In order to match the two requirements, the Müntz polynomials can be considered as alternative bases. The reasons behind this is to apply the critical feature that a fractional derivative of a Müntz polynomial is again a Müntz polynomial. Therefore, we can hope to obtain a good approximation for the fractional derivatives by the Müntz polynomial approximation. Due to these outstanding features, Müntz polynomials have been widely used for the solution of fractional equations in literature. Esmaeili [50] provided the solution for fractional differential equations with the Müntz polynomial collocation method. Mokhtary [51] solved the fractional problems with the Müntz polynomial Tau method. Bahmanpour discussed the Müntz polynomial wavelets collocation method for fractional equations [52]. Recently, Maleknejad discussed the Müntz–Legendre wavelet approach [53]. The Müntz polynomial has also been absorbed in the BSM to solve fractional equations [46].

The remainder of this paper is organized as follows. Section 2 contains a brief definition of the problems to be solved and also a brief description of the solution process. Section 3 contains the derivation of the analytical approximations satisfying the general boundary conditions. This technique for the orthogonal basis is described in detail in Section 4. Following the main algorithm in Section 5, numerical examples that illustrate the presented procedure are placed in Section 6. Finally, a brief conclusion is drawn in Section 7.

2. Preliminaries

In the present work, our goal is to find an effective solution to the following multi-dimensional multi-term time-fractional partial differential equations (TFPDEs):

$${\mathcal{L}}_t\left[u\right]={\mathcal{M}}_t\left[{\nabla }^{2}u\right]+f\left(\mathbf{x},t\right),t\in \left[0,T\right],\mathbf{x}\in {\mathsf{\Omega }}^d={\left[0,1\right]}^d,d=1,2,3,$$

(1)

where

$${\mathcal{L}}_t=D_t^{\left(\mu \right)}+\sum _{k=1}^I{a}_k\left(t\right)D_t^{\left({\mu }_k\right)},{\mathcal{M}}_t=\sum _{k=I+1}^K{a}_k\left(t\right)D_t^{\left({\mu }_k\right)},$$

(2)

in which $\mu \in \left(l-1,l\right]$, $0\le {\mu }_k<\mu$, $D_t^{\left(0\right)}\left[\phi \right]\equiv \phi$ is the identical operator, and $a_k\left(t\right)$, $k=1,\dots ,K$. Some initial conditions (ICs) should be prescribed in advance for the time-dependent problems

$$u\left(\mathbf{x},0\right)=u_0\left(\mathbf{x}\right),\frac{\partial u\left(\mathbf{x},0\right)}{\partial t}=u_1\left(\mathbf{x}\right),\dots ,\frac{{\partial }^{l-1}u\left(\mathbf{x},0\right)}{\partial t^{l-1}}=u_{l-1}\left(\mathbf{x}\right),$$

(3)

where l is the highest integer derivative of the considered problem.

The operator $D_t^{\left(\nu \right)}$, which has the following form,

$$D_t^{\left(\nu \right)}\left[\xi (x,t)\right]=\left\{\begin{array}{ccc}\frac{1}{\mathsf{\Gamma }\left(n-\nu \right)}{\int }_0^x\frac{{\partial }_t^{\left(n\right)}\xi \left(x,\tau \right)d\tau }{{\left(t-\tau \right)}^{\nu -n+1}},\hfill & \hfill & n-1<\nu <n,\hfill \\ {\partial }_t^{\left(n\right)}\xi \left(x,t\right),\hfill & \hfill & \nu =n,\hfill \end{array}\right.$$

(4)

is the Caputo fractional derivative of the order $\nu$. In particular, if $\xi (x,t)$ is the power function $t^z$, we have

$$D_t^{\left(\nu \right)}\left[t^z\right]=0,$$

(5)

$\mathrm{if}z\in {\mathcal{N}}_0\mathrm{and}z<n$, and

$$D_t^{\left(\nu \right)}\left[t^z\right]=\frac{\mathsf{\Gamma }\left(z+1\right)}{\mathsf{\Gamma }\left(z+1-\nu \right)}{t}^{z-\nu },$$

(6)

if $z\in {\mathcal{N}}_0\mathrm{and}z\ge n\mathrm{or}z\notin {\mathcal{N}}_0\mathrm{and}z>n-1$. Further, ${\mathcal{N}}_0$ denotes the set of all non-negative integers.

In order to solve Equation (1), suitable boundary conditions have to be prescribed to ensure the solvability of the problems. In the present work, we address general forms of boundary conditions:

$$B_i\left[u\right]={\alpha }_i\frac{\partial u\left(\mathbf{x},t\right)}{\partial n_i}+{\beta }_{i}u\left(\mathbf{x},t\right)=g_i\left(\mathbf{x},t\right),\mathbf{x}\in \partial {\mathsf{\Omega }}^d,{\mathsf{\Omega }}^d={\left[0,1\right]}^d,{\alpha }_i^2+{\beta }_i^2\ne 0,$$

(7)

where $\partial u\left(\mathbf{x},t\right)/\partial n_i$ denotes differentiation along the outward unit normal direction of the surface boundary, and $i=1,2$, for $d=1$, $i=1,\dots ,4$ for $d=2$, and $i=1,\dots ,6$ for $d=3.$

It is known to all that in mathematics and engineering applications, the Fourier expansion in eigenfunctions of the differential operator is an efficient numerical method in the case of homogeneous BCs when the solution can be represented as a linear combination of eigenfunctions ${\psi }_{n_1}^{\left(1\right)}\left(x_1\right),{\psi }_{n_2}^{\left(2\right)}\left(x_2\right),{\psi }_{n_3}^{\left(3\right)}\left(x_3\right)$,

$$u\left(\mathbf{x},t\right)=\sum _{n_1,n_2,n_3=1}^{\infty }{U}_{n_1,n_2,n_3}\left(t\right){\psi }_{n_1}^{\left(1\right)}\left(x_1\right){\psi }_{n_2}^{\left(2\right)}\left(x_2\right){\psi }_{n_3}^{\left(3\right)}\left(x_3\right),$$

(8)

and the unknowns can be determined by substituting the expression into the initial conditions using the orthogonality property of eigenfunctions with different eigenvalues. It should be emphasized that, for the application of the Fourier method, we have to transform the original problem into a homogeneous problem. In this case, our main goal is to calculate the analytic function $v_g\left(\mathbf{x},t\right)$ that exactly satisfies the boundary conditions of Equation (7) for any given ${\alpha }_i$, ${\beta }_i$, $g_i\left(\mathbf{x},t\right)$ at each t in Equation (7). This function can be used to solve the problem of the non-homogeneous boundary conditions cardinally. Suppose that the solution can be approximated by the following approximation:

$$u\left(\mathbf{x},t\right)=v_g\left(\mathbf{x},t\right)+w\left(\mathbf{x},t\right).$$

(9)

Substituting the above equation into Equations (1), (3), and (7), we have:

$${\mathcal{L}}_t\left[w\right]={\mathcal{M}}_t\left[{\nabla }^{2}w\right]+f_1\left(\mathbf{x},t\right),t\in \left[0,T\right],\mathbf{x}\in {\mathsf{\Omega }}^d,$$

(10)

$$\frac{{\partial }^{i}w\left(\mathbf{x},0\right)}{\partial t^i}=w_i\left(\mathbf{x}\right),i=0,\dots ,l-1,$$

(11)

$$B_i\left[w\right]={\alpha }_i\frac{\partial w\left(\mathbf{x},t\right)}{\partial n_i}+{\beta }_{i}w\left(\mathbf{x},t\right)=0,\mathbf{x}\in \partial {\mathsf{\Omega }}^d.$$

(12)

It is evident that the boundary conditions have been transformed into the homogeneous one, which makes it possible to use the Fourier-series expansion as follows:

$$w\left(\mathbf{x},t\right)=\sum _{n=1}^{\infty }{w}_{\mathbf{n}}\left(t\right){\psi }_{\mathbf{n}}\left(\mathbf{x}\right),$$

(13)

where the orthonormal basis ${\psi }_{\mathbf{n}}\left(\mathbf{x}\right)$ is corresponding to the BC Equation (12), which satisfies

$${\nabla }^2{\psi }_{\mathbf{n}}\left(\mathbf{x}\right)=-{\lambda }_{\mathbf{n}}^2{\psi }_{\mathbf{n}}\left(\mathbf{x}\right),\mathbf{x}\in {\mathsf{\Omega }}^d,B_i\left[{\psi }_{\mathbf{n}}\left(\mathbf{x}\right)\right]=0,\mathbf{x}\in \partial {\mathsf{\Omega }}^d.$$

(14)

This orthonormal basis will be described in Section 4. By substituting Equation (13) into Equation (10) and projecting $〈\dots ,{\psi }_{\mathbf{n}}〉$, we have

$${\mathcal{L}}_t\left[w_{\mathbf{n}}\left(t\right)\right]=-{\lambda }_{\mathbf{n}}^2{\mathcal{M}}_t\left[w_{\mathbf{n}}\left(t\right)\right]+{\theta }_{\mathbf{n}}\left(t\right),t\in \left[0,T\right]$$

(15)

for the approximation of each $w_{\mathbf{n}}$, which will be described in detail in the following several sections.

3. The Algorithm for Computing ${\mathit{v}}_{\mathit{g}}\left(\mathbf{x},\mathit{t}\right)$

Construction of the analytic function $v_g\left(\mathbf{x},t\right)$, which exactly satisfies the BCs of the original problem, is the main subject of the proposed method. In this section, we will propose an approach to derive $v_g\left(\mathbf{x},t\right)$ for the $(1+1)$-dimensional problems, $(2+1)$-dimensional problems, and $(3+1)$-dimensional problems, respectively.

3.1. $(1+1)$-Dimensional Problems

In this case, the problem of finding the function $v_g\left(x,t\right)$, which conforms the BC

$$\begin{array}{c}{L}_W\left[v_g\left(x,t\right)\right]\left(x=0\right)={\alpha }_W\frac{\partial v_g}{\partial x}\left(0,t\right)+{\beta }_W{v}_g\left(0,t\right)=g_W\left(t\right),\end{array}$$

(16)

$$\begin{array}{c}{L}_E\left[v_g\left(x,t\right)\right]\left(x=1\right)={\alpha }_E\frac{\partial v_g}{\partial x}\left(1,t\right)+{\beta }_E{v}_g\left(1,t\right)=g_E\left(t\right),\end{array}$$

(17)

at the endpoints of the interval ${\mathsf{\Omega }}^1=\left[0,1\right]$ is a trivial one. Indeed, one can prove easily that the following functions

$$\begin{array}{c}{\theta }_E\left(x\right)=\frac{{\alpha }_W-{\beta }_{W}x}{{\alpha }_W{\beta }_E-{\beta }_W\left({\alpha }_E+{\beta }_E\right)},\end{array}$$

(18)

$$\begin{array}{c}{\theta }_W\left(x\right)=\frac{{\beta }_{E}x-\left({\alpha }_E+{\beta }_E\right)}{{\alpha }_W{\beta }_E-{\beta }_W\left({\alpha }_E+{\beta }_E\right)},\end{array}$$

(19)

satisfy the conditions

$$\begin{array}{c}{L}_W\left(x\right)\left[{\theta }_W\left(x\right)\right]\left(0\right)=1,L_E\left(x\right)\left[{\theta }_W\left(x\right)\right]\left(1\right)=0,\end{array}$$

(20)

$$\begin{array}{c}{L}_W\left(x\right)\left[{\theta }_E\left(x\right)\right]\left(0\right)=0,L_E\left(x\right)\left[{\theta }_E\left(x\right)\right]\left(1\right)=1.\end{array}$$

(21)

Then, the function

$$v_g\left(x,t\right)={\theta }_E\left(x\right)g_E\left(t\right)+{\theta }_W\left(x\right)g_W\left(t\right)$$

(22)

satisfies the BCs of Equations (16) and (17).

3.2. $(2+1)$-Dimensional Problems

For the $(2+1)$-dimensional problem, we intend to obtain the $v_g$ satisfying the BCs in a square domain; for example:

$$\begin{array}{c}{L}_W\left(x\right)\left[v_g\right]\equiv {\alpha }_W\frac{\partial v_g}{\partial x}+{\beta }_W{v}_g=g_W\left(y,t\right),x=0,0\le y\le 1,\end{array}$$

(23)

$$\begin{array}{c}{L}_E\left(x\right)\left[v_g\right]\equiv {\alpha }_E\frac{\partial v_g}{\partial x}+{\beta }_E{v}_g=g_E\left(y,t\right),x=1,0\le y\le 1,\end{array}$$

(24)

$$\begin{array}{c}{L}_S\left(y\right)\left[v_g\right]\equiv {\alpha }_S\frac{\partial v_g}{\partial y}+{\beta }_S{v}_g=g_S\left(x,t\right),0\le x\le 1,y=0,\end{array}$$

(25)

$$\begin{array}{c}{L}_N\left(y\right)\left[v_g\right]\equiv {\alpha }_N\frac{\partial v_g}{\partial y}+{\beta }_N{v}_g=g_N\left(x,t\right),0\le x\le 1,y=1.\end{array}$$

(26)

We assume that the desired function $v_g\left(x,y,t\right)$ is smooth enough and so the functions $g_W\left(y,t\right)$, $g_E\left(y,t\right)$, $g_S\left(x,t\right)$, and $g_N\left(x,t\right)$ guarantee continuity condition at the apexes $\left(0,0\right)$, $\left(0,1\right)$, $\left(1,0\right)$, and $\left(1,1\right)$ of the unit square ${\left[0,1\right]}^2$:

$$\begin{array}{c}{L}_S\left(y\right)\left[g_W\left(y,t\right)\right]=L_W\left(x\right)\left[g_S\left(x,t\right)\right]\mathrm{at}\left(0,0\right),\end{array}$$

(27)

$$\begin{array}{c}{L}_N\left(y\right)\left[g_W\left(y,t\right)\right]=L_W\left(x\right)\left[g_N\left(x,t\right)\right]\mathrm{at}\left(0,1\right),\end{array}$$

(28)

$$\begin{array}{c}{L}_S\left(y\right)\left[g_E\left(y,t\right)\right]=L_E\left(x\right)\left[g_S\left(x,t\right)\right]\mathrm{at}\left(1,0\right),\end{array}$$

(29)

$$\begin{array}{c}{L}_N\left(y\right)\left[g_E\left(y,t\right)\right]=L_E\left(x\right)\left[g_N\left(x,t\right)\right]\mathrm{at}\left(1,1\right).\end{array}$$

(30)

Let us define the functions

$${\theta }_N\left(y\right)=\frac{{\alpha }_S-{\beta }_{S}y}{{\alpha }_S{\beta }_N-{\beta }_S({\alpha }_N+{\beta }_N)},$$

(31)

$${\theta }_S\left(y\right)=\frac{{\beta }_{N}y-({\alpha }_N+{\beta }_N)}{{\alpha }_S{\beta }_N-{\beta }_S({\alpha }_N+{\beta }_N)},$$

(32)

which are similar to the functions ${\theta }_E\left(x\right)$, ${\theta }_W\left(x\right)$. They satisfy the boundary conditions

$$\begin{array}{c}{L}_S\left(y\right)\left[{\theta }_S\left(y\right)\right]\left(0\right)=1,L_N\left(y\right)\left[{\theta }_S\left(y\right)\right]\left(1\right)=0,\end{array}$$

(33)

$$\begin{array}{c}{L}_S\left(y\right)\left[{\theta }_N\left(y\right)\right]\left(0\right)=0,L_N\left(y\right)\left[{\theta }_N\left(y\right)\right]\left(1\right)=1.\end{array}$$

(34)

Let us define the function

$$v_1={\theta }_E\left(x\right)g_E+{\theta }_W\left(x\right)g_W.$$

(35)

One can easily prove that $v_1$ satisfies Equations (23) and (24):

$$\begin{array}{c}{L}_W\left(x\right){\left[v_1\left(x,y,t\right)\right]}_{x=0}=g_W\left(y,t\right),0\le y\le 1,\end{array}$$

(36)

$$\begin{array}{c}{L}_E\left(x\right){\left[v_1\left(x,y,t\right)\right]}_{x=1}=g_E\left(y,t\right),0\le y\le 1.\end{array}$$

(37)

This follows directly from Definitions (20) and (21). Additionally, we define $g_{N1}\left(x,t\right)$ and $g_{S1}\left(x,t\right)$ as follows:

$$g_{N1}\left(x,t\right)=g_N\left(x,t\right)-L_N\left(y\right){\left[v_1\left(x,y,t\right)\right]}_{y=1},$$

(38)

$$g_{S1}\left(x,t\right)=g_S\left(x,t\right)-L_S\left(y\right){\left[v_1\left(x,y,t\right)\right]}_{y=0}.$$

(39)

Finally, we can prove that the following combination,

$$v_g\left(x,y,t\right)=v_1\left(x,y,t\right)+{\theta }_N\left(y\right)g_{N1}\left(x,t\right)+{\theta }_S\left(y\right)g_{S1}\left(x,t\right),$$

(40)

satisfies Equations (23)–(26).

3.3. $(3+1)$-Dimensional Problems

For $(3+1)$-dimensional problems, we try to seek a smooth analytical function $v_g$ that satisfies

$$\begin{array}{c}{L}_W\left(x\right)\left[v_g\right]\equiv {\alpha }_W\frac{\partial v_g}{\partial x}+{\beta }_W{v}_g=g_W\left(y,z,t\right),x=0,\end{array}$$

(41)

$$\begin{array}{c}{L}_E\left(x\right)\left[v_g\right]\equiv {\alpha }_E\frac{\partial v_g}{\partial x}+{\beta }_E{v}_g=g_E\left(y,z,t\right),x=1,\end{array}$$

(42)

$$\begin{array}{c}{L}_S\left(y\right)\left[v_g\right]\equiv {\alpha }_S\frac{\partial v_g}{\partial y}+{\beta }_S{v}_g=g_S\left(x,z,t\right),y=0,\end{array}$$

(43)

$$\begin{array}{c}{L}_N\left(y\right)\left[v_g\right]\equiv {\alpha }_N\frac{\partial v_g}{\partial y}+{\beta }_N{v}_g=g_N\left(x,z,t\right),y=1,\end{array}$$

(44)

$$\begin{array}{c}{L}_T\left(z\right)\left[v_g\right]\equiv {\alpha }_T\frac{\partial v_g}{\partial z}+{\beta }_T{v}_g=g_T\left(x,y,t\right),z=1,\end{array}$$

(45)

$$\begin{array}{c}{L}_B\left(z\right)\left[v_g\right]\equiv {\alpha }_B\frac{\partial v_g}{\partial z}+{\beta }_B{v}_g=g_B\left(x,y,t\right),z=0,\end{array}$$

(46)

where $0\le x,y,z\le 1$ if the variables $x,y,z$ are not defined in the above equations. Without loss of generality, the operators $L_E\left(x\right)$ and $L_N\left(y\right)$ hold the general property $L_E\left(x\right)\left[L_N\left(y\right)\left[u\right]\right]=L_N\left(y\right)\left[L_E\left(x\right)\left[u\right]\right]$ in the presence of any given smooth function u for $x=1,y=1,0\le z\le 1$. It follows that

$$\begin{array}{c}{L}_E\left(x\right){\left[L_N\left(y\right)\left[u\right]\right]}_{x=1,y=1}=L_N\left(y\right){\left[L_E\left(x\right)\left[u\right]\right]}_{x=1,y=1}\Rightarrow L_E\left(x\right){\left[L_N\left(y\right){\left[u\right]}_{y=1}\right]}_{x=1}\hfill \\ \hfill =L_N\left(y\right){\left[L_E\left(x\right){\left[u\right]}_{x=1}\right]}_{y=1}\Rightarrow L_E\left(x\right){\left[g_N\left(x,z,t\right)\right]}_{x=1}=L_N\left(y\right){\left[g_E\left(y,z,t\right)\right]}_{y=1}.\end{array}$$

(47)

The above condition is obviously fulfilled for the remaining 11 edges. Let us define the functions ${\theta }_T\left(z\right),$${\theta }_B\left(z\right)$ as

$${\theta }_T\left(z\right)=\frac{{\alpha }_B-{\beta }_{B}z}{{\beta }_T{\alpha }_B-{\beta }_B({\alpha }_T+{\beta }_T)},$$

(48)

$${\theta }_B\left(z\right)=\frac{{\beta }_{T}z-({\alpha }_T+{\beta }_T)}{{\beta }_T{\alpha }_B-{\beta }_B({\alpha }_T+{\beta }_T)},$$

(49)

similar to Equations (18), (19), (31), and (32), which satisfy the boundary conditions

$$\begin{array}{c}{L}_T\left(z\right)\left[{\theta }_T\left(z\right)\right]\left(1\right)=1,L_B\left(z\right)\left[{\theta }_T\left(z\right)\right]\left(0\right)=0,\end{array}$$

(50)

$$\begin{array}{c}{L}_T\left(z\right)\left[{\theta }_B\left(z\right)\right]\left(1\right)=0,L_B\left(z\right)\left[{\theta }_B\left(z\right)\right]\left(0\right)=1.\end{array}$$

(51)

Let us now try to construct the set of auxiliary functions along with their linear combinations

$$\begin{array}{c}{v}_1={\theta }_T\left(z\right)g_T\left(x,y,t\right)+{\theta }_B\left(z\right)g_B\left(x,y,t\right),\end{array}$$

(52)

$$\begin{array}{c}{g}_{N1}=g_N\left(x,z,t\right)-L_N\left(y\right){\left[v_1\left(x,y,z,t\right)\right]}_{y=1},\end{array}$$

(53)

$$\begin{array}{c}{g}_{S1}=g_S\left(x,z,t\right)-L_S\left(y\right){\left[v_1\left(x,y,z,t\right)\right]}_{y=0},\end{array}$$

(54)

$$v_2=v_1+{\theta }_N\left(y\right)g_{N1}+{\theta }_S\left(y\right)g_{S1},$$

(55)

where ${\theta }_N\left(y\right)$ and ${\theta }_S\left(y\right)$ can be defined according to Equations (31) and (32). Then, let us define

$$\begin{array}{c}{g}_{E1}\left(y,z,t\right)=g_E\left(y,z,t\right)-L_E\left(x\right){\left[v_2\left(x,y,z,t\right)\right]}_{x=1},\end{array}$$

(56)

$$\begin{array}{c}{g}_{W1}\left(y,z,t\right)=g_W\left(y,z,t\right)-L_W\left(x\right){\left[v_2\left(x,y,z,t\right)\right]}_{x=0},\end{array}$$

(57)

and, finally,

$$v_g(x,y,z,t)=v_2(x,y,z,t)+{\theta }_E\left(x\right)g_{E1}\left(y,z,t\right)+{\theta }_W\left(x\right)g_{W1}\left(y,z,t\right),$$

(58)

which determines the objective function that satisfies the BC Equations (41)–(46). Note also that there are not any problems in the derivation of $v_g$ for various boundary conditions, including the boundary conditions of the third kind. Furthermore, it is noteworthy that the $v_g$ obtained in this section is not unique. In the present work, we show a simple form of the formulas for easy application. Other methods such as numerical methods can also be used for such purposes that are not reported in the present work.

4. Orthogonal Basis for the Fourier Expansion

This section deals with the application of the Fourier method to the solution of corresponding homogeneous problems. As mentioned earlier, once we obtain the function $v_g(\mathbf{x},t)$, the original equation can be transformed into the following homogeneous system by substituting the $v_g(\mathbf{x},t)$ into the governing equations and boundary conditions

$${\mathcal{L}}_t\left[w\right]={\mathcal{M}}_t\left[{\nabla }^{2}w\right]+f_1\left(\mathbf{x},t\right),f_1\left(\mathbf{x},t\right)=f\left(\mathbf{x},t\right)+{\mathcal{M}}_t\left[{\nabla }^2{v}_g\left(\mathbf{x},t\right)\right]-{\mathcal{L}}_t\left[v_g\left(\mathbf{x},t\right)\right],$$

(59)

which satisfies homogeneous BC

$$B_i\left[w\right]={\alpha }_i\frac{\partial w\left(\mathbf{x},t\right)}{\partial n_i}+{\beta }_{i}w\left(\mathbf{x},t\right)=0,\mathbf{x}\in \partial {\mathsf{\Omega }}^d,$$

(60)

where the solution can now be approximated with the following functional form,

$$w\left(\mathbf{x},t\right)=\sum _{n=1}^{\infty }{w}_{\mathbf{n}}\left(t\right){\psi }_{\mathbf{n}}\left(\mathbf{x}\right),$$

(61)

over the orthonormal basis ${\psi }_{\mathbf{n}}\left(\mathbf{x}\right)$, corresponding to BC Equation (60),

$${\nabla }^2{\psi }_{\mathbf{n}}\left(\mathbf{x}\right)=-{\lambda }_{\mathbf{n}}^2{\psi }_{\mathbf{n}}\left(\mathbf{x}\right),B_i\left[{\psi }_{\mathbf{n}}\left(\mathbf{x}\right)\right]=0,\mathbf{x}\in \partial {\mathsf{\Omega }}^d.$$

(62)

In this section, we describe the orthonormal basis ${\psi }_{\mathbf{n}}\left(\mathbf{x}\right)$ and begin with the $(1+1)$-dimensional problem as an example.

4.1. (1 + 1)-Dimensional Problems

Let us consider the following Sturm–Liouville problem:

$$\begin{array}{c}\frac{d^2\psi }{d{x}^2}=-\mu \psi ,\end{array}$$

(63)

$$\begin{array}{c}{L}_W\left(x\right){\left[\psi \right]}_{x=0}={\left({\alpha }_W\frac{d\psi }{dx}+{\beta }_W\psi \right)}_{x=0}=0,L_E\left(x\right){\left[\psi \right]}_{x=1}={\left({\alpha }_E\frac{d\psi }{dx}+{\beta }_E\psi \right)}_{x=1}=0,\end{array}$$

(64)

where we describe some possible forms of the eigenfunctions corresponding to Equations (63) and (64).

• First, let us consider the general case ${\alpha }_W\ne 0$, ${\alpha }_E\ne 0$. We write the boundary conditions in the traditional form for transport problems:

  $$\frac{d\psi }{dx}-b_1\psi {|}_{x=0}=0,\frac{d\psi }{dx}+b_2\psi {|}_{x=1}=0,b_1,b_2\ge 0.$$

  (65)

The change of the sign is connected to the different direction of the outward normal at the endpoints $x=0$ and $x=1$ [41].

From Equation (63), it follows:

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \mu {\int }_0^1{\psi }^{2}dx& =\hfill & \hfill -{\int }_0^1\frac{d^2\psi }{d{x}^2}\psi dx={\int }_0^1\left[{\left(\frac{d\psi }{dx}\right)}^{2}dx-d\left(\psi \frac{d\psi }{dx}\right)\right]\\ & =\hfill & \hfill -\left(\psi \frac{d\psi }{dx}\right){|}_0^1+{\int }_0^1{\left(\frac{d\psi }{dx}\right)}^{2}dx\\ & =\hfill & \hfill \psi \left(0\right)\frac{d\psi \left(0\right)}{dx}-\psi \left(1\right)\frac{d\psi \left(1\right)}{dx}+{\int }_0^1{\left(\frac{d\psi }{dx}\right)}^{2}dx.\end{array}\end{array}$$

(66)

Using Equation (65), one gets:

$$\mu {\int }_0^1{\psi }^{2}dx=b_1{\psi }^2\left(0\right)+b_2{\psi }^2\left(1\right)+{\int }_0^1{\left(\frac{d\psi }{dx}\right)}^{2}dx.$$

(67)

Thus, under the condition $b_1,b_2\ge 0$, the values of $\mu$ are positive, and we denote $\mu ={\lambda }^2$. The function

$${\psi }_{1,n}\left(x\right)=\frac{1}{{\mathcal{R}}_n}\left[cos\left({\lambda }_{n}x\right)+\frac{b_1}{{\lambda }_n}sin\left({\lambda }_{n}x\right)\right],$$

(68)

where

$${\mathcal{R}}_n=\sqrt{\frac{1}{2}\left(1+\frac{b_1^2}{{\lambda }_n^2}\right)+\frac{b_1}{2{\lambda }_n^2}+\frac{b_2}{2{\lambda }_n^2}\frac{{\lambda }_n^2+b_1^2}{{\lambda }_n^2+b_2^2}}$$

(69)

and where ${\lambda }_n$ is the $n^{th}$ solution that can be obtained by solving the following system of transcendental equation

$$\left({\lambda }^2-b_1{b}_2\right)tan\left(\lambda \right)=\left(b_1+b_2\right)\lambda ,$$

(70)

where ${\psi }_{1,n}\left(x\right)$ constructs an orthonormal basis in the Hilbert space $L\left(\left[0,1\right]\right)$, and the following identity holds for different bases

$$〈{\psi }_{1,n},{\psi }_{1,m}〉={\int }_0^1{\psi }_{1,n}\left(x\right){\psi }_{1,m}\left(x\right)dx={\delta }_{n,m}.$$

(71)

2.

Let us consider the case ${\alpha }_W=0$, ${\alpha }_E\ne 0$. We get the Dirichlet condition at the left endpoint $x=0$:

$${\psi |}_{x=0}=0,\frac{d\psi }{dx}+b_2{\psi |}_{x=1}=0,b_2\ge 0.$$

(72)

The function

$${\psi }_{2,n}\left(x\right)=\frac{1}{{\mathcal{R}}_n}sin\left({\lambda }_{n}x\right),$$

(73)

is an eigenfunction of Sturm–Liouville Problems (63) and (72). Here

$${\mathcal{R}}_n=\sqrt{\frac{1}{2}\left(1-\frac{sin\left(2{\lambda }_n\right)}{2{\lambda }_n}\right)},$$

(74)

and ${\lambda }_n$ is the $n^{th}$ solution of the transcendental equation

$$\lambda cos\left(\lambda \right)+b_{2}sin\left(\lambda \right)=0.$$

(75)

Function ${\psi }_{2,n}\left(x\right)$ satisfies $〈{\psi }_{2,n},{\psi }_{2,m}〉={\delta }_{n,m}$.

3.

Let us consider the case ${\alpha }_W\ne 0$, ${\alpha }_E=0$. We get the Dirichlet condition at the right endpoint $x=1$:

$$\frac{d\psi }{dx}-b_1\psi {{|}_{x=0}=0,\psi |}_{x=1}=0,b_1\ge 0.$$

(76)

The function

$${\psi }_{3,n}\left(x\right)=\frac{1}{{\mathcal{R}}_n}\left[sin\left({\lambda }_{n}x\right)+\frac{{\lambda }_n}{b_1}cos\left({\lambda }_{n}x\right)\right]$$

(77)

is an eigenfunction of Sturm–Liouville Problems (63) and (76). Here

$${\mathcal{R}}_n=\sqrt{\frac{1}{2}\left(1-\frac{sin\left(2{\lambda }_n\right)}{2{\lambda }_n}\right)+\frac{1}{2b_1}\left(1-cos\left(2{\lambda }_n\right)\right)+\frac{{\lambda }_n^2}{2b_1^2}\left(1+\frac{sin\left(2{\lambda }_n\right)}{2{\lambda }_n}\right)},$$

(78)

and ${\lambda }_n$ is the $n^{th}$ solution of the transcendental equation

$$\lambda cos\left(\lambda \right)+b_{1}sin\left(\lambda \right)=0.$$

(79)

Function ${\psi }_{3,n}\left(x\right)$ satisfies $〈{\psi }_{3,n},{\psi }_{3,m}〉={\delta }_{n,m}$.

4.

The case ${\alpha }_W=0$, ${\alpha }_E=0$ corresponds to the Dirichlet conditions at both endpoints of the interval $\left[0,1\right]$. The function

$${\psi }_{4,n}\left(x\right)=\sqrt{2}sinn\pi x,n=1,2,\dots ,$$

(80)

also satisfies the property $〈{\psi }_{4,n},{\psi }_{4,m}〉={\delta }_{n,m}$.

5.

Finally, let us consider the following case ${\alpha }_W$, ${\alpha }_E$, ${\beta }_E\ne 0$, ${\beta }_W=0$. Using Equation (65), the boundary conditions can be expressed as

$$\frac{d\psi }{dx}=0,\frac{d\psi }{dx}+b_2{\psi |}_{x=1}=0,b_1,b_2\ge 0.$$

(81)

Thus, we get the Neumann condition at the left-hand side endpoint. The function

$${\psi }_{5,n}\left(x\right)=\frac{1}{{\mathcal{R}}_n}cos\left({\lambda }_{n}x\right),$$

(82)

satisfies $〈{\psi }_{5,n},{\psi }_{5,m}〉={\delta }_{n,m}$. Here

$${\mathcal{R}}_n=\sqrt{\frac{1}{2}\left(1+\frac{b_2}{{\lambda }_n^2+b_2^2}\right)},$$

(83)

and ${\lambda }_n$ is the $n^{th}$ solution of the transcendental equation

$$\lambda sin\left(\lambda \right)=b_{2}cos\left(\lambda \right).$$

(84)

4.2. $(2+1)$ and $(3+1)$-Dimensional Problems

We use the products ${\psi }_{i_1,n_1}\left(x\right){\psi }_{i_2,n_2}\left(y\right)$ and ${\psi }_{i_1,n_1}\left(x\right){\psi }_{i_2,n_2}\left(y\right){\psi }_{i_3,n_3}\left(z\right)$ as the basis function for solving $(2+1)$ and $(3+1)$-dimensional problems. Here, the first index $i_1,i_2,i_3=1,2,3,4$, 5 indicates the type of the basis function as described in the last subsection; the second index $n_1,n_2,n_3$ indicates the harmonic number. These functions satisfy the equations:

$${\nabla }_{x,y}^2{\psi }_{i_1,n_1}\left(x\right){\psi }_{i_2,n_2}\left(y\right)=-\left({\lambda }_{{}_{i_1,n_1}}^2+{\lambda }_{{}_{i_2,n_2}}^2\right){\psi }_{i_1,n_1}\left(x\right){\psi }_{i_2,n_2}\left(y\right),$$

(85)

$${\nabla }_{x,y,z}^2{\psi }_{i_1,n_1}\left(x\right){\psi }_{i_2,n_2}\left(y\right){\psi }_{i_3,n_3}\left(z\right)=-\left({\lambda }_{{}_{i_1,n_1}}^2+{\lambda }_{{}_{i_2,n_2}}^2+{\lambda }_{{}_{i_3,n_3}}^2\right){\psi }_{i_1,n_1}\left(x\right){\psi }_{i_2,n_2}\left(y\right){\psi }_{i_3,n_3}\left(z\right).$$

(86)

Below they are used in the Fourier transformation of governing Equation (59).

5. The Solution Procedure

In this section, we demonstrate the solution procedure for the proposed method for the multi-dimensional fractional equations by utilizing the analytic function $v_g$ and the orthogonal basis derived in the last section. At first, we discuss the solution for the one-dimensional problem with a single harmonic.

5.1. $(1+1)$ Problem with a Single Spatial Harmonic

Now, for the $(1+1)$-dimensional TFPDE with a specific source term,

$$L_t\left[w\right]=M_t\left[\frac{{\partial }^{2}w}{\partial x^2}\right]+{\theta }_n\left(t\right){\psi }_n\left(x\right),$$

(87)

${\psi }_n\left(x\right)$ denotes the eigenfunctions given by Equations (68), (73), (77), (80), or (82). Suppose that the equation is subjected to the BCs and ICs, which conform the chosen eigenfunction ${\psi }_n\left(x\right)$:

$$L_W\left(x\right){\left[w\right]}_{x=0}=0,L_E\left(x\right){\left[w\right]}_{x=1}=0,$$

(88)

$$w\left(x,0\right)=w_0{\psi }_n\left(x\right),\frac{\partial w\left(x,0\right)}{\partial t}=w_1{\psi }_n\left(x\right),\dots ,\frac{{\partial }^{l-1}w\left(x,0\right)}{\partial t^{l-1}}=w_{l-1}{\psi }_n\left(x\right).$$

(89)

We see the solution as

$$w\left(x,t\right)=w_n\left(t\right){\psi }_n\left(x\right).$$

(90)

Then, we get the multi-term FODE

$$L_t\left[w_n\right]=-{\lambda }_n^2{M}_t\left[w_n\right]+{\theta }_n\left(t\right).$$

(91)

Using Equation (2), we rewrite the equation in the expanded form

$$D_t^{\left(\mu \right)}\left[w_n\right]=-\sum _{k=1}^K{b}_k\left(t\right)D_t^{\left({\mu }_k\right)}\left[w_n\right]+{\theta }_n\left(t\right),$$

(92)

where $b_k\left(t\right)=a_k\left(t\right)$, $k=1,\dots ,I$, $b_k\left(t\right)={\lambda }_n^2{a}_k\left(t\right)$, $k=I+1,\dots ,K$, and $\mu \in \left(l-1,l\right],0\le {\mu }_k<\mu$. The considered initial conditions from Equation (89) are as follows:

$$w_n\left(0\right)=w_0,\frac{d{w}_n\left(0\right)}{dt}=w_1,\dots ,\frac{d^{l-1}{w}_n\left(0\right)}{d{t}^{l-1}}=w_{l-1}.$$

(93)

Supposing that the right-hand side of Equation (92) can be approximated by the linear combination of $\phi$, as shown below,

$$-\sum _{k=1}^K{b}_k\left(t\right)D_t^{\left({\mu }_k\right)}\left[w_n\right]+{\theta }_n\left(t\right)=\sum _{m=1}^{\infty }{q}_m{\phi }_m\left(t\right),$$

(94)

where the $\phi$ is the chosen basis and $q_m$ are unknown coefficients to be determined. In this way, Equation (92) is transformed into

$$D_t^{\left(\mu \right)}\left[w_n\left(t\right)\right]=\sum _{m=1}^{\infty }{q}_m{\phi }_m\left(t\right),$$

(95)

where ${\phi }_m\left(t\right)$ and ${\varphi }_m\left(t\right)$ holds

$$D_t^{\left(\mu \right)}\left[{\varphi }_m\left(t\right)\right]={\phi }_m\left(t\right).$$

(96)

It is easy to verify that the Müntz polynomial basis (MPB) proposed in [46,50,51,52,53] meets the requirement given in Equation (96). The explicit expression of such functions is as follows:

$${\phi }_m\left(t\right)=t^{{\delta }_m},{\delta }_m=\sigma \left(m-1\right),$$

(97)

where $\sigma$ is an auxiliary variable within $(0,1]$, and m is a positive integer. From the preceding numerical examples in [44], we set $\sigma =0.25$ for this paper.

It follows that the function

$${\varphi }_m\left(t\right)=\frac{\mathsf{\Gamma }\left({\delta }_m+1\right)}{\mathsf{\Gamma }\left({\delta }_m+\mu +1\right)}{t}^{{\delta }_m+\mu }$$

(98)

satisfies the following equality:

$$D^{\left(\mu \right)}\left[{\varphi }_m\left(t\right)\right]={\phi }_m\left(t\right).$$

(99)

As far as $l-1<\mu \le l$, the function ${\varphi }_m\left(t\right)$ satisfies zero ICs:

$${\varphi }_m^{\left(i\right)}\left(0\right)=0,i=0,1,\dots ,l-1.$$

(100)

From Equations (97), (99), and (100) it can be seen that the series

$$w_{h,n}\left(t,\mathbf{q}\right)=\sum _{m=1}^{\infty }{q}_m{\varphi }_m\left(t\right)$$

(101)

satisfies Equation (95) with any ${\left\{q_m\right\}}_{m=1}^{\infty }\equiv \mathbf{q}$.

Let us define the following approximation:

$$w_{p,n}\left(t\right)=w_0+w_{1}t+\cdots +w_{l-1}\frac{t^{l-1}}{\left(l-1\right)!}=\sum _{i=0}^{l-1}\frac{w_i}{i!}{t}^i.$$

(102)

The function $w_{p,n}\left(t\right)$ satisfies IC Equation (93). In this way, the following approximation

$$w_n\left(t,\mathbf{q}\right)=w_{p,n}\left(t\right)+w_{h,n}\left(t,\mathbf{q}\right)=w_{p,n}\left(t\right)+\sum _{m=1}^{\infty }{q}_m{\varphi }_m\left(t\right)$$

(103)

satisfies Equations (95) and (93) with any ${\left\{q_m\right\}}_{m=1}^{\infty }$. Additionally, the unknown weighted parameter $q_m$ is determined by backward substituting the $w_n\left(t,\mathbf{q}\right)$ into Equation (94):

$$\sum _{m=1}^{\infty }{q}_m\left[{\phi }_m\left(t\right)+\sum _{k=1}^K{b}_k\left(t\right)D_t^{\left({\mu }_k\right)}\left[{\varphi }_m\left(t\right)\right]\right]={\theta }_n\left(t\right)-\sum _{k=1}^K{b}_k\left(t\right)D_t^{\left({\mu }_k\right)}\left[w_p\left(t\right)\right]\equiv F_n\left(t\right),$$

(104)

where we denote

$$\begin{array}{c}{D}_t^{\left({\mu }_k\right)}\left[{\varphi }_m\left(t\right)\right]=\frac{\mathsf{\Gamma }\left({\delta }_m+1\right)D_t^{\left({\mu }_k\right)}\left[t^{{\delta }_m+\mu }\right]}{\mathsf{\Gamma }\left({\delta }_m+\mu +1\right)}\hfill \\ =\frac{\mathsf{\Gamma }\left({\delta }_m+1\right)\mathsf{\Gamma }\left({\delta }_m+\mu +1\right)t^{{\delta }_m+\mu -{\mu }_k}}{\mathsf{\Gamma }\left({\delta }_m+\mu +1\right)\mathsf{\Gamma }\left({\delta }_m+\mu +1-{\mu }_k\right)}=\frac{\mathsf{\Gamma }\left({\delta }_m+1\right)t^{{\delta }_m+\mu -{\mu }_k}}{\mathsf{\Gamma }\left({\delta }_m+\mu +1-{\mu }_k\right)},\hfill \end{array}$$

(105)

$$D_t^{\left({\mu }_k\right)}\left[w_p\left(t\right)\right]=\sum _{i=1}^{l-1}\frac{w_i}{i!}{D}_t^{\left({\mu }_k\right)}\left[t^i\right]=\sum _{i=1}^{l-1}\frac{w_i\mathsf{\Gamma }\left(i+1\right)t^{i-{\mu }_k}}{i!\mathsf{\Gamma }\left(i+1-{\mu }_k\right)}=\sum _{i\ge {\mu }_k}^{l-1}\frac{w_i{t}^{i-{\mu }_k}}{\mathsf{\Gamma }\left(i+1-{\mu }_k\right)}.$$

(106)

We consider the truncated series of Equation (103),

$$w_n\left(t,M,\mathbf{q}\right)=w_{p,n}\left(t\right)+\sum _{m=1}^M{q}_m{\varphi }_m\left(t\right),$$

(107)

which also satisfies the modified Equation (95):

$$D^{\left(\mu \right)}\left[w_n\left(t,M,\mathbf{q}\right)\right]=\sum _{m=1}^M{q}_m{\phi }_m\left(t\right).$$

(108)

The unknown parameters $q_1,\dots ,q_M$ should satisfy the truncated version of Equation (104),

$$\sum _{m=1}^M{q}_m\left[{\phi }_m\left(t\right)+\sum _{k=1}^K{b}_k\left(t\right)D_t^{\left({\mu }_k\right)}\left[{\varphi }_m\left(t\right)\right]\right]=F_n\left(t\right),$$

(109)

by the collocation method as follows:

$$\sum _{m=1}^M{q}_m\left[{\phi }_m\left(t_j\right)+\sum _{k=1}^K{b}_k\left(t_j\right)D_t^{\left({\mu }_k\right)}\left[{\varphi }_m\left(t_j\right)\right]\right]=F_n\left(t_j\right),$$

(110)

where

$$t_j=0.5T\left[1+cos\left(\pi \left(2j-1\right)/2N_c\right)\right]\in \left[0,T\right],j=1,2,\dots N_c\ge M,$$

(111)

are the Gauss–Chebyshev (GC) collocation points on the time interval $\left[0,T\right]$. It is important to note that in the framework of the presented method, the $F_n\left(t\right)$ are required at several time steps $t_j$ given in Equation (111) only. As it follows from Equation (104), the same is true for time function ${\theta }_n\left(t\right)$. As a result, the Fourier expansion of the $f_1\left(\mathbf{x},t\right)$ of Equation (10) also should be performed at the same time moments $t_j$ only.

5.2. $(1+1)$-Dimensional Problems of the General Case

Consider the following equation:

$$L_t\left[u\right]=M_t\left[\frac{{\partial }^{2}u}{\partial x^2}\right]+f\left(x,t\right),x\in \left[0,1\right],t\in \left[0,T\right].$$

(112)

The substitution

$$u=v_g+w,$$

(113)

with $v_g\left(x,t\right)$, gives us

$$L_t\left[w\right]=M_t\left[\frac{{\partial }^{2}w}{\partial x^2}\right]+f_1\left(x,t\right),$$

(114)

$$L_W\left(x\right){\left[w\right]}_{x=0}=0,L_E\left(x\right){\left[w\right]}_{x=1}=0,$$

(115)

and IC

$$\frac{{\partial }^{i}w\left(x,0\right)}{\partial t^i}=u_i\left(x\right)-\frac{{\partial }^i{v}_g\left(x,0\right)}{{\partial }^{i}t}=w_i\left(x\right),i=0,\dots ,l-1.$$

(116)

Here,

$$f_1\left(x,t\right)=f\left(x,t\right)-L_t\left[v_g\left(x,t\right)\right]+M_t\left[\frac{{\partial }^2{v}_g\left(x,t\right)}{\partial x^2}\right].$$

(117)

We seek the solution of the problem in Equations (114)–(116) using the following linear combination,

$$w_N\left(x,t\right)=\sum _{n=1}^N{w}_n\left(t\right){\psi }_n\left(x\right),$$

(118)

where ${\psi }_n\left(x\right)$ denotes one of the eigenfunctions given by Equations (68), (73), (77), (80), and (82). Substituting Equation (118) into Equation (114), we have

$$L_t\left[w_n\right]=-{\lambda }_n^2{M}_t\left[w_n\right]+{\theta }_n\left(t\right),n=1,\dots ,N,$$

(119)

where

$${\theta }_n\left(t\right)={\int }_0^1{f}_1\left(x,t\right){\psi }_n\left(x\right)dx.$$

(120)

Further, $w_n$ holds

$$\frac{d^i{w}_n\left(0\right)}{d{t}^i}={\int }_0^1{w}_i\left(x\right){\psi }_n\left(x\right)dx=w_{i,n},i=0,\dots ,l-1,$$

(121)

which follows from Equation (116). The solution of $w_n$ can be given by

$$w_n\left(t,M,{\mathbf{q}}_n\right)=w_{p,n}\left(t\right)+\sum _{m=1}^M{q}_{n,m}{\varphi }_m\left(t\right),{\mathbf{q}}_n={\left\{q_{n,m}\right\}}_{m=1}^M.$$

(122)

Therefore, the solution $u_{appro}\left(x,t\right)$ to the origin problem is given by

$$u_{appro}=v_g\left(x,t\right)+\sum _{n=1}^N{w}_{p,n}\left(t\right){\psi }_n\left(x\right)+\sum _{m=1}^M{Q}_m\left(x\right){\varphi }_m\left(t\right),$$

(123)

where the $Q_m\left(x\right)$ is given by

$$Q_m\left(x\right)=\sum _{n=1}^N{q}_{n,m}{\psi }_n\left(x\right).$$

(124)

5.3. $(2+1)$-Dimensional Problems

Consider the following $(2+1)$-dimensional problem:

$$L_t\left[u\right]=M_t\left[{\nabla }^{2}u\right]+f\left(x,y,t\right),\left(x,y\right)\in {\left[0,1\right]}^2,t\in \left[0,T\right].$$

(125)

The substitution

$$u=v_g+w,$$

(126)

with $v_g\left(x,y,t\right)$ given by Equation (40) gives us

$$L_t\left[w\right]=M_t\left[{\nabla }^{2}w\right]+f_1\left(x,y,t\right),$$

(127)

and the homogeneous BCs

$$L_W\left(x\right){\left[w\right]}_{x=0}=0,L_E\left(x\right){\left[w\right]}_{x=1}=0,L_S\left(y\right){\left[u\right]}_{y=0}=0,L_N\left(y\right){\left[u\right]}_{y=1}=0,$$

(128)

and ICs

$$\frac{{\partial }^{i}w\left(x,y,0\right)}{\partial t^i}=u_i\left(x,y\right)-\frac{{\partial }^i{v}_g\left(x,y,0\right)}{\partial t^i}=w_i\left(x,y\right),i=0,\dots ,l-1.$$

(129)

The approximate solution to Problems (127)–(129) is approximated as follows:

$$w_N\left(\mathbf{x},t\right)=\sum _{n_1,n_2=1}^N{w}_{n_1,n_2}\left(t\right){\psi }_{n_1,n_2}\left(x_1,x_2\right)=\sum _{\mathbf{n}=1}^N{w}_{\mathbf{n}}\left(t\right){\psi }_{\mathbf{n}}\left(\mathbf{x}\right).$$

(130)

Here, ${\psi }_{n_1,n_2}\left(x_1,x_2\right)={\psi }_{\mathbf{n}}\left(\mathbf{x}\right)={\psi }_{n_1}\left(x_1\right){\psi }_{n_2}\left(x_2\right)$ is the product of the eigenfunctions given by Equations (68), (73), (77), (80), or (82), and we use the following short notations: $\mathbf{n}=\left(n_1,n_2\right)$, $\mathbf{x}=\left(x_1,x_2\right)=\left(x,y\right)$.

Substituting $w_N\left(\mathbf{x},t\right)$ in (127), we have

$${\mathcal{L}}_t\left[w_{\mathbf{n}}\left(t\right)\right]=-\left({\lambda }_{n_1}^2+{\lambda }_{n_2}^2\right){\mathcal{M}}_t\left[w_{\mathbf{n}}\left(t\right)\right]+{\theta }_{\mathbf{n}}\left(t\right),t\in \left[0,T\right],$$

(131)

where

$${\theta }_{\mathbf{n}}\left(t\right)={\int }_0^1{\int }_0^1{f}_1\left(x_1,x_2,t\right){\psi }_{n_1,n_2}\left(x_1,x_2\right)d{x}_{1}d{x}_2.$$

(132)

The harmonic $w_{\mathbf{n}}\left(t\right)=w_{n_1,n_2}\left(t\right)$ satisfies

$$\frac{{\partial }^i{w}_{0,\mathbf{n}}\left(0\right)}{\partial t^i}=w_{i,\mathbf{n}}={\int }_0^1{\int }_0^1{w}_i\left(x_1,x_2\right){\psi }_{n_1,n_2}\left(x_1,x_2\right)d{x}_{1}d{x}_2,i=0,1,\dots ,l-1,$$

(133)

which follows from Equation (129). The approximate solution

$$w_{\mathbf{n}}\left(t,M,{\mathbf{q}}_{\mathbf{n}}\right)=w_{p,\mathbf{n}}\left(t\right)+\sum _{m=1}^M{q}_{\mathbf{n},m}{\varphi }_m\left(t\right),{\mathbf{q}}_{\mathbf{n}}={\left\{{\mathbf{q}}_{\mathbf{n},m}\right\}}_{m=1}^M={\left\{q_{n_1,n_2,m}\right\}}_{m=1}^M,$$

(134)

can be obtained obviously for each harmonic. Then, the $u_{N_{,}M}\left(\mathbf{x},t\right)$ can be approximated as follows:

$$\begin{array}{c}{u}_{N_{,}M}\left(\mathbf{x},t\right)=v_g\left(\mathbf{x},t\right)+{\displaystyle \sum _{n_1,n_2=1}^N}{w}_{M,n_1,n_2}\left(t\right){\psi }_{n_1,n_2}\left(x_1,x_2\right)\hfill \\ =v_g\left(\mathbf{x},t\right)+{\displaystyle \sum _{n_1,n_2=1}^N}{w}_{p,n_1,n_2}\left(t\right){\psi }_{n_1,n_2}\left(x_1,x_2\right)+{\displaystyle \sum _{m=1}^M}{Q}_m\left(\mathbf{x}\right){\varphi }_m\left(t\right),\hfill \end{array}$$

(135)

where

$$Q_m\left(\mathbf{x}\right)=\sum _{n_1,n_2=1}^N{q}_{n_1,n_2,m}{\psi }_{n_1,n_2}\left(x_1,x_2\right).$$

(136)

5.4. $(3+1)$-Dimensional Problems

Now let us move to the $(3+1)$-dimensional problems

$$L_t\left[u\right]=M_t\left[{\nabla }^{2}u\right]+f,\left(x,y,z\right)\in {\left[0,1\right]}^3,t\in \left[0,T\right].$$

(137)

Using the substitution

$$u=v_g+w,$$

(138)

with $v_g$ given by Equation (58), we have

$$L_t\left[w\right]=M_t\left[{\nabla }^{2}w\right]+f_1\left(x,y,z,t\right),$$

(139)

subjected to

$$\begin{array}{c}{L}_W\left(x\right){\left[w\right]}_{x=0}=0,L_E\left(x\right){\left[w\right]}_{x=1}=0,L_S\left(y\right){\left[u\right]}_{y=0}=0,L_N\left(y\right){\left[u\right]}_{y=1}=0,\hfill \\ L_B\left(z\right){\left[u\right]}_{z=0}=0,L_T\left(z\right){\left[u\right]}_{z=1}=0,\hfill \end{array}$$

(140)

and IC

$$\frac{{\partial }^{i}w\left(x,y,z,0\right)}{\partial t^i}=u_i\left(x,y,z\right)-\frac{{\partial }^i{v}_g\left(x,y,z,0\right)}{\partial t^i}=w_i\left(x,y,z\right),i=0,\dots ,l-1,$$

(141)

in which

$$f_1=f-L_t\left[v_g\right]+M_t\left[{\nabla }^2{v}_g\right].$$

(142)

With the same technology, the approximated solution of $w_N\left(\mathbf{x},t\right)$ can be expressed in the following functional form,

$$w_N\left(\mathbf{x},t\right)=\sum _{\mathbf{n}=1}^N{w}_{\mathbf{n}}\left(t\right){\psi }_{\mathbf{n}}\left(\mathbf{x}\right),$$

(143)

where $\mathbf{n}=\left(n_1,n_2,n_3\right),$ $\mathbf{x}=\left(x,y,z\right)=\left(x_1,x_2,x_3\right)$, ${\psi }_{\mathbf{n}}\left(\mathbf{x}\right)={\psi }_{n_1,n_2,n_3}\left(x_1,x_2,x_3\right)$ is the product of the eigenfunctions ${\psi }_{n_1}\left(x\right){\psi }_{n_2}\left(y\right){\psi }_{n3}\left(z\right)$ given by Equations (68), (73), (77), (80), or (82).

Substituting $w_N\left(\mathbf{x},t\right)$ in Equation (127), we have

$${\mathcal{L}}_t\left[w_{\mathbf{n}}\left(t\right)\right]=-\left({\lambda }_{n_1}^2+{\lambda }_{n_2}^2+{\lambda }_{n_3}^2\right){\mathcal{M}}_t\left[w_{\mathbf{n}}\left(t\right)\right]+{\theta }_{\mathbf{n}}\left(t\right),t\in \left[0,T\right],$$

(144)

where

$${\theta }_{\mathbf{n}}\left(t\right)={\int }_0^1{\int }_0^1{\int }_0^1{f}_1\left(\mathbf{x},t\right){\psi }_{n_1,n_2,n_3}\left(\mathbf{x}\right)d\mathbf{x}.$$

(145)

It is important to note that in the context of the present method, the Fourier expansion of the source term $f_1\left(\mathbf{x},t\right)$ over the eigenfunction basis should be performed at several fixed time moments $t_j$ only. This follows from the algorithm of the backward substitution technique (see Equations (104) and (121)).

The harmonic $w_{\mathbf{n}}\left(t\right)=w_{n_1,n_2,n_3}\left(t\right)$ satisfies

$$\frac{{\partial }^i{w}_{0,\mathbf{n}}\left(0\right)}{\partial t^i}=w_{i,\mathbf{n}}={\int }_0^1{\int }_0^1{\int }_0^1{w}_i\left(\mathbf{x}\right){\psi }_{\mathbf{n}}\left(\mathbf{x}\right)d\mathbf{x},i=0,\dots ,l-1.$$

(146)

Then, the solution is approximated as

$$\begin{array}{c}{u}_{N_{,}M}\left(\mathbf{x},t\right)=v_g+{\displaystyle \sum _{\mathbf{n}=\left(1,1,1\right)}^N}{w}_{M,\mathbf{n}}\left(t\right){\psi }_{\mathbf{n}}\left(\mathbf{x}\right)\hfill \\ =v_g+{\displaystyle \sum _{\mathbf{n}=\left(1,1,1\right)}^N}{w}_{p,\mathbf{n}}\left(t\right){\psi }_{\mathbf{n}}\left(\mathbf{x}\right)+{\displaystyle \sum _{m=1}^M}{Q}_m\left(\mathbf{x}\right){\varphi }_m\left(t\right),\hfill \end{array}$$

(147)

where

$$Q_m\left(\mathbf{x}\right)=\sum _{n\mathbf{n}=\left(1,1,1\right)}^N{q}_{\mathbf{n},m}{\psi }_{\mathbf{n}}\left(\mathbf{x}\right).$$

(148)

6. Numerical Examples

In this section, the feasibility of the proposed method is experimentally verified. The maximum absolute error (MAE) and the $E_{H^1}\left(t\right)$ error containing the derivatives were used as numerical criteria, as shown below:

$$E_{max}\left(t\right)=\underset{1\le i\le N_t}{max}\left|u_{N,M}\left({\mathbf{x}}_i,t\right)-u_{ex}({\mathbf{x}}_i,t)\right|,$$

(149)

$$E_{H^1}\left(t\right)=\sqrt{\frac{1}{N_t}{\displaystyle \sum _{i=1}^N}\left[{\left(u_{N,M}\left({\mathbf{x}}_i,t\right)-u_{ex}\left({\mathbf{x}}_i,t\right)\right)}^2+{\left(\frac{\partial u_{N,M}}{\partial t}\left({\mathbf{x}}_i,t\right)-\frac{\partial u_{ex}}{\partial t}\left({\mathbf{x}}_i,t\right)\right)}^2\right]},$$

(150)

where $u_{N,M}\left({\mathbf{x}}_i,t\right)$ indicates the approximate solutions obtained by the presented analytical–numerical method for the compared solution $u_{ex}\left({\mathbf{x}}_i,t\right)$, and $N_t$ is the total number of test nodes.

For the 1D problems, we used the number of test nodes $N_t=4N$, where N is the number of spatial harmonics; i.e., we use 4 testing nodes per harmonic. For the 2D examples, the errors were carried using the $N_t=4000$ test nodes distributed in the solution domain. For the 3D problem, we have transformed it into the FODE analytically so that we only have to check the solution accuracy in the time domain.

As for the total number of collocation nodes, we have to illustrate that, in this paper, the derivation of $v_g$ can be done analytically. Therefore, we do not have to place nodes on the boundary. Using the approximate solution in the form of the Fourier series over the eigenfunction, we transform the TFPDE into the set of the FODEs for each of the Fourier harmonics. Therefore, we do not have to place collocation nodes inside the domain. The collocation nodes are placed in the time domain only.

The collocation points $t_j$, $j=1,2,\dots ,N_c,$ in the time interval $\left[0,T\right]$ are used to form the collocation system for solving each FODE. In all the examples, we use the number of collocation nodes $N_c=2M$ in the time domain, where M is the number of the Müntz polynomials that are used in the approximate solution of the FODEs. The parameter M defines the accuracy of the approximation in time.

6.1. $(1+1)$-Dimensional Problem

6.1.1. Example 1

For the first example, the following time fractional cable equation is studied under the Dirichlet boundary condition

$$\frac{\partial u\left(x,t\right)}{\partial t}=D_t^{(1-{\gamma }_1)}\left[\frac{{\partial }^{2}u}{\partial x^2}\left(x,t\right)\right]-D_t^{(1-{\gamma }_2)}\left[u\left(x,t\right)\right]+f\left(x,t\right),x\in \left(0,1\right),t\in \left(0,1\right),$$

(151)

$$u\left(x,0\right)=0,u\left(0,t\right)=0,u\left(1,t\right)=0,$$

(152)

where the source term $f(x,t)$ can be computed by substituting the analytic solution $u\left(x,t\right)$

$$u\left(x,t\right)=t^{2}sin\pi x,$$

(153)

into the governing equation, which yields

$$f\left(x,t\right)=2\left(t+\frac{{\pi }^2{t}^{1+{\gamma }_1}}{\mathsf{\Gamma }\left(2+{\gamma }_1\right)}+\frac{t^{1+{\gamma }_2}}{\mathsf{\Gamma }\left(2+{\gamma }_2\right)}\right)sin\pi x=F\left(t\right)sin\pi x.$$

(154)

Thus, the problem considered is a special 1D case with a single spatial harmonic, which was considered in Section 4.1. As it is shown there, the problem can be reduced to a single FODE,

$$\frac{dw\left(t\right)}{dt}=-{\pi }^2{D}_t^{(1-{\gamma }_1)}\left[w\left(t\right)\right]-D_t^{(1-{\gamma }_2)}\left[w\left(t\right)\right]+F\left(t\right),$$

(155)

for the sole harmonic.

To illustrate the effects of the error in M,

Table 1. The MAE versus the M at $t=1$ for different $\sigma$.

M	1	2	3	4	5	6	10
$\sigma$ = 0.25	1.19 × 10${}^{-2}$	1.93 × 10${}^{-3}$	1.97 × 10${}^{-4}$	3.73 × 10${}^{-5}$	1.22 × 10${}^{-16}$	1.22 × 10${}^{-16}$	-
$\sigma$ = 0.50	1.19 × 10${}^{-2}$	1.47 × 10${}^{-3}$	1.22 × 10${}^{-16}$	1.22 × 10${}^{-16}$	-	-	-
$\sigma$ = 0.33	1.19 × 10${}^{-2}$	1.83 × 10${}^{-3}$	1.21 × 10${}^{-4}$	9.35 × 10${}^{-7}$	5.44 × 10${}^{-8}$	6.38 × 10${}^{-9}$	3.12 × 10${}^{-1}$

shows the behavior of the maximum absolute error with respect to M for the approximation of the source terms in Equations (97), (98), and (107). The approximate solution is sought as a truncated series Equation (107) over the function ${\varphi }_m\left(t\right)$: $D^{\left(\mu \right)}\left[{\varphi }_m\left(t\right)\right]={\phi }_m\left(t\right)$ and so belongs to the linear span $S\left(\mu ,\sigma ,M\right)=Span{\left(1,t^{\mu +\sigma \left(m-1\right)}\right)}_{m=1}^M$. In the case considered, $\mu =1$ and $S\left(1,\sigma ,M\right)=Span{\left(1,t^{1+\sigma \left(m-1\right)}\right)}_{m=1}^M$. For $\sigma =0.25$, $S\left(1,0.25,M\right)=$$Span\left(1,t,t^{1.25},\dots ,t^{1+0.25\left(M-1\right)}\right).$ Therefore, as it comes from Equation (153), the exact solution $w\left(t\right)=t^2$ belongs to $S\left(1,0.25,M\right)$ for $M\ge 5$. The data in Table 1 demonstrate that, for this particular case, the results of the proposed analytical–numerical method converge to the exact solution for the parameter $M\ge 5$ and reach the computer rounding errors. Let us consider the case $\sigma =0.5$. Here, $w\left(t\right)=t^2$ belongs to $S\left(1,0.5,M\right)$ for $M\ge 3$. The data in Table 1 illustrate this situation.

It is easy to verify that for $\sigma =0.33$ (the bottom row of the table), there is no such M that the exact solution $w\left(t\right)$ belongs to $S\left(1,0.33,M\right)$. As a result, for $\sigma =0.33$, the error decreases gradually with the growth of M, while for $\sigma =0.25$ and $\sigma =0.50$, it decreases sharply when the exact solution belongs to the corresponding range $S\left(\alpha ,\sigma ,M\right)$. The calculations have been performed for $\left({\gamma }_1=0.5,{\gamma }_2=0.5\right)$, $\left({\gamma }_1=0.3,{\gamma }_2=0.9\right),$ and $\left({\gamma }_1=0.7,{\gamma }_2=0.6\right)$ which have produced the same results. Thus, if the parameters $\sigma$ of the Müntz polynomial basis are chosen in such a way that exact solution belongs to $S\left(\mu ,\sigma ,M\right),$ then the present method provides the exact solution up to the rounding errors of the computer. This problem has been considered by Yang, Jiang, and Zhang in [54] using the spectral Legendre–Tau method. The most accurate result that has been achieved there has the error E${}_{max}$ = 9.3019 × 10${}^{-6}$ when 13 Legendre’s polynomials were used for the spatial and temporal approximations.

6.1.2. Example 2

Let us consider the following problem that has been studied using the time-space spectral tau method in [54]:

$$\frac{\partial u\left(x,t\right)}{\partial t}=D_t^{(1-{\gamma }_1)}\left[\frac{{\partial }^{2}u}{\partial x^2}\left(x,t\right)\right]-D_t^{(1-{\gamma }_2)}\left[u\left(x,t\right)\right]+f\left(x,t\right),$$

(156)

subjected to the following conditions:

$$u\left(0,t\right)-\frac{\partial u\left(0,t\right)}{\partial x}=-\pi t^{2+{\gamma }_1},$$

(157)

$$u\left(1,t\right)+\frac{\partial u\left(1,t\right)}{\partial x}=-\pi t^{2+{\gamma }_1}+2e{t}^{1+{\gamma }_2}.$$

(158)

The BCs, the source term, and the IC can be computed from the corresponding exact solution:

$$u\left(x,t\right)=t^{2+{\gamma }_1}sin\pi x+t^{1+{\gamma }_2}{e}^x.$$

(159)

In order to show the effects of N and M on the accuracy,

Table 2. The MAE, $C{O}_N$, and the elapsed time versus the N using M = 7 (left) and $M=15$ (right).

	$\mathit{M}=7$			$\mathit{M}=15$
$\mathit{N}$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}$	${\mathit{C}\mathit{O}}_{\mathit{N}}$	$\mathit{C}\mathit{P}\mathit{U}$, $\mathit{s}$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}$	${\mathit{C}\mathit{O}}_{\mathit{N}}$	$\mathit{C}\mathit{P}\mathit{U}$, $\mathit{s}$
2	9.37 × 10${}^{-2}$	-	0.02	9.37 × 10${}^{-2}$	-	0.11
4	7.31 × 10${}^{-3}$	3.68	0.07	7.31 × 10${}^{-3}$	3.68	0.17
8	6.81 × 10${}^{-4}$	3.42	0.16	6.81 × 10${}^{-4}$	3.42	0.33
16	7.30 × 10${}^{-5}$	3.22	0.36	7.22× 10${}^{-5}$	3.24	0.66
32	9.09 × 10${}^{-6}$	3.01	0.61	8.30 × 10${}^{-6}$	3.12	1.0
64	1.78 × 10${}^{-6}$	2.35	0.98	9.94 × 10${}^{-7}$	3.06	1.9
128	9.97 × 10${}^{-7}$	0.84	1.84	1.22 × 10${}^{-7}$	3.03	4.1
256	9.94 × 10${}^{-7}$	0.04	3.18	1.51 × 10${}^{-8}$	3.01	9.0
512	9.94 × 10${}^{-7}$	<<1	6.5	1.89 × 10${}^{-9}$	2.99	19
1024	9.94 × 10${}^{-7}$	<<1	13	2.57 × 10${}^{-10}$	2.88	44
2048	9.94 × 10${}^{-7}$	<<1	24	5.36 × 10${}^{-11}$	2.26	74
[54], Table 2, $E_{max}=$3.2279 × 10${}^{-5}$

displays the MAE, elapsed computational time, and the order of convergence with respect to parameter N:

$$C{O}_N=\frac{log\left(E_{max}\left(N\right)/E_{max}\left(2N\right)\right)}{log2},$$

(160)

for $M=7$ and $M=15$ with $\sigma =0.25$, ${\gamma }_1=0.6$, and ${\gamma }_2=0.9$. From this table, it is clearly seen that, for the case $M=7$, the error decreases shapely with the increase of parameter N up to the value $N=128$. For larger $N>128$, the results of the proposed method do not change, and the solution accuracy remains at ${10}^{-7}$. In the case of $M=15$, the error decreases monotonically over the whole range of N, and the final accuracy comes to ${10}^{-11}$. The order of convergence is three. Table 2 tabulates the solutions given in [54] by the usage of the spectral Legendre–Tau method for comparison. The data correspond to the case where 13 Legendre’s polynomials are used for the spatial and temporal approximations. From the comparison, it can be seen obviously that the proposed analytical–numerical method leads to a better solution even for small values of M and N from the point of view of standard accuracy.

Table 3. The MAE and $C{O}_M$ versus the M with the fixed number of harmonics N.

	$\mathit{N}=32$		$\mathit{N}=128$		$\mathit{N}=512$
$\mathit{M}$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}$	${\mathit{C}\mathit{O}}_{\mathit{M}}$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}$	${\mathit{C}\mathit{O}}_{\mathit{M}}$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}$	${\mathit{C}\mathit{O}}_{\mathit{M}}$
2	1.87 × 10${}^{-1}$	-	1.87 × 10${}^{-1}$	-	1.87 × 10${}^{-1}$	-
3	4.02 × 10${}^{-2}$	3.79	4.02 × 10${}^{-2}$	3.79	4.02 × 10${}^{-2}$	3.79
4	6.69 × 10${}^{-3}$	6.23	6.69 × 10${}^{-3}$	6.23	6.69 × 10${}^{-3}$	6.23
5	6.66 × 10${}^{-4}$	10.3	6.66 × 10${}^{-4}$	10.3	6.66 × 10${}^{-4}$	10.3
6	1.82 × 10${}^{-5}$	19.7	1.24 × 10${}^{-5}$	21.8	1.24 × 10${}^{-5}$	21.8
7	9.09 × 10${}^{-6}$	4.5	9.97 × 10${}^{-7}$	16.4	9.97 × 10${}^{-7}$	16.4
8	8.33 × 10${}^{-6}$	0.66	1.51 × 10${}^{-7}$	14.1	3.72 × 10${}^{-8}$	24.6
9	8.30 × 10${}^{-6}$	<<1	1.27 × 10${}^{-7}$	1.47	7.42 ×10${}^{-9}$	13.7
10	8.30 × 10${}^{-6}$	<<1	1.23 × 10${}^{-7}$	0.08	3.22 × 10${}^{-9}$	7.92
11	8.30 × 10${}^{-6}$	<<1	1.22 × 10${}^{-7}$	<<1	2.33 × 10${}^{-9}$	3.39
12	8.30 × 10${}^{-6}$	<<1	1.22 × 10${}^{-7}$	<<1	1.96 × 10${}^{-9}$	0.63

displays the MAE and convergence order with respect to the parameter M for the fixed values N = 32, 128, 512. Here the convergence order is defined as:

$$C{O}_M=\frac{log\left(E_{max}\left(M_1\right)/E_{max}\left(M_2\right)\right)}{log\left(M_2/M_1\right)}.$$

(161)

When M is small enough, it defines the accuracy of the approximate solution for all values of N. For $N=32$, the accuracy remains the same for the growth of $M>8$. On the other hand, when $N=512$, the calculated error decreases with the increase of M in the whole range $2\le M\le 12$. This means that 512 Fourier’s harmonics provide an accurate approximation and the main error here is caused by the solving of FODEs. From the last row of this table, it can be seen evidently that the proposed method has high rates of convergence for the MAE, which provides reasonably accurate approximations for the unknown variables. It should be noted here that, with the increasing of M and N, the CO becomes flat, which means that the errors do not change for very large M or N. The proposed algorithm converges to the stable results.

In Figure 1, the observed behavior of the error is shown in more detail. Let us consider the left-hand side of Figure 1. The graphics $log\left(Emax\right(N\left)\right)$ have the same origin for all fixed M. With the growth of N, the curves $log\left(Emax\right(N\left)\right)$ change shape depending on the fixed value of M.

6.1.3. Example 3

The third example solved here is the high order TFPDE

$$\begin{array}{c}{D}_t^{\left(\alpha \right)}\left[u\left(x,t\right)\right]+\frac{sin\left(t\right)}{1+t}{D}_t^{\left({\alpha }_1\right)}\left[u\left(x,t\right)\right]+\frac{cos\left(t\right)}{1+t}{D}_t^{\left({\alpha }_2\right)}\left[u\left(x,t\right)\right]+log\left(1+t\right)u\left(x,t\right)\hfill \\ \hfill =\frac{sinh\left(t\right)}{1+t}{D}_t^{\left({\alpha }_3\right)}\left[\frac{{\partial }^{2}u}{\partial x^2}\left(x,t\right)\right]+\frac{cosh\left(t\right)}{1+t}{D}_t^{\left({\alpha }_4\right)}\left[\frac{{\partial }^{2}u}{\partial x^2}\left(x,t\right)\right]+\left(1+t^2\right)\frac{{\partial }^{2}u}{\partial x^2}\left(x,t\right)+f\left(x,t\right),\end{array}$$

(162)

where $\alpha =\sqrt{21}$, ${\alpha }_1=\pi$, ${\alpha }_2=\sqrt{5}$, ${\alpha }_3=\sqrt{3}$, ${\alpha }_4$=$\sqrt{2}$. As far as $4<\alpha =\sqrt{21}<5$, the TFPDE needs the following five ICs:

$$\frac{{\partial }^{i}u}{\partial t^i}\left(x,t=0\right)=u_i\left(x\right),i=0,1,2,3,4.$$

(163)

The equation is subjected to the BC

$$\frac{\partial u}{\partial x}\left(x=0,t\right)-{\pi }^{2}u(x=0,t)=g_W\left(t\right),\frac{\partial u}{\partial x}\left(x=1,t\right)+e^{2}u(x=1,t)=g_E\left(t\right),$$

(164)

where the functions $f\left(x,t\right)$, $u_i\left(x\right)$, $g_W\left(t\right)$, $g_E\left(t\right)$ can be easily computed from the following exact solution

$$u\left(x,t\right)=cos\left(t\right)\left[cos\left(x\right)+cosh\left(x\right)\right].$$

(165)

Table 4. The MAE of the u and $\partial u/\partial t$ versus N using M = 16 (left) and 24 (right).

	$\mathit{M}=16$			$\mathit{M}=24$
$\mathit{N}$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}\left(\mathit{u}\right)$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}\left(\partial \mathit{u}/\partial \mathit{t}\right)$	${\mathit{E}}_{{\mathit{H}}^{\mathit{1}}}$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}\left(\mathit{u}\right)$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}\left(\partial \mathit{u}/\partial \mathit{t}\right)$	${\mathit{E}}_{{\mathit{H}}^1}$
100	3.53 × 10${}^{-8}$	5.51 × 10${}^{-8}$	8.43 × 10${}^{-9}$	3.53 × 10${}^{-8}$	5.51 × 10${}^{-8}$	8.42 × 10${}^{-9}$
200	4.36 × 10${}^{-9}$	6.92 × 10${}^{-9}$	8.16 × 10${}^{-10}$	4.39 × 10${}^{-9}$	6.84 × 10${}^{-9}$	7.33 × 10${}^{-10}$
300	1.27 × 10${}^{-9}$	2.10 × 10${}^{-9}$	3.93 × 10${}^{-10}$	1.30 × 10${}^{-9}$	2.02 × 10${}^{-9}$	1.76 × 10${}^{-10}$
400	5.21 × 10${}^{-10}$	9.31 × 10${}^{-10}$	3.56 × 10${}^{-10}$	5.47 × 10${}^{-10}$	8.51 × 10${}^{-10}$	6.42 × 10${}^{-11}$
500	2.54 × 10${}^{-10}$	5.15 × 10${}^{-10}$	3.51 × 10${}^{-10}$	2.80 × 10${}^{-10}$	4.36 × 10${}^{-10}$	2.94 × 10${}^{-11}$
1000	5.80 × 10${}^{-11}$	1.94 × 10${}^{-10}$	3.50 × 10${}^{-10}$	3.49 × 10${}^{-11}$	5.44 × 10${}^{-11}$	2.59 × 10${}^{-12}$
2000	5.80 × 10${}^{-11}$	1.94 × 10${}^{-10}$	3.50 × 10${}^{-10}$	4.35 × 10${}^{-12}$	6.81 × 10${}^{-12}$	2.65 × 10${}^{-13}$

presents the maximum absolute errors of the solution and its first derivative in time with the increase in N for $M=16$ and $M=24$. In the case when $M=24$, the errors decrease sharply for $100\le N\le 2000$. For $M=16$, the errors are the same for $N>500$. In Figure 2, this behavior of the error is shown in more detail.

6.2. $(2+1)$-Dimensional Problems

6.2.1. Example 4

Let us consider the following two-dimensional multi-term time-fractional mixed sub-diffusion and diffusion-wave equation defined in the unit square

$$\begin{array}{c}{D}_t^{\left(\alpha \right)}\left[u\left(x,y,t\right)\right]+\frac{\partial u}{\partial t}\left(x,y,t\right)+D_t^{\left({\alpha }_1\right)}\left[u\left(x,y,t\right)\right]+u\left(x,y,t\right)\hfill \\ ={\nabla }^{2}u\left(x,y,t\right)+D_t^{\left({\alpha }_2\right)}\left[{\nabla }^{2}u\left(x,y,t\right)\right]+f\left(x,y,t\right),\left(x,y\right)\in {\left[0,1\right]}^2,\hfill \end{array}$$

(166)

with the exact solution

$$u\left(x,y,t\right)=\left(1+t^3\right)sin\left(\pi x\right)sin\left(\pi y\right)\equiv w\left(t\right){\psi }_{1,1}\left(x,y\right),$$

(167)

for the case $\alpha =1.6$, ${\alpha }_1=0.7$, and ${\alpha }_2=0.3$.

The initial and Dirichlet boundary conditions of function $f\left(x,y,t\right)$ conform the exact solution. Thus, we get the TFPDE with a single spatial harmonic corresponding to the eigenfunction ${\psi }_{1,1}\left(x,y\right)=2sin\left(\pi x\right)sin\left(\pi y\right)$. As it is shown above, the problem can be reduced to the single FODE

$$\begin{array}{c}{D}_t^{\left(\alpha \right)}\left[w\left(t\right)\right]+\frac{dw\left(t\right)}{dt}+D_t^{\left({\alpha }_1\right)}\left[w\left(t\right)\right]+w\left(t\right)\hfill \\ =-2{\pi }^{2}w\left(t\right)-2{\pi }^2{D}_t^{\left({\alpha }_2\right)}\left[w\left(t\right)\right]+F\left(t\right),\hfill \end{array}$$

(168)

with initial conditions $w\left(0\right)=w_0,$$dw\left(0\right)/dt=w_1.$

Table 5. The MAE and the computational time versus the M at $t=1$.

M	1	2	3	4	5	6	10	15
$\sigma$ = 0.25	4.3 × 10${}^{-1}$	4.8 × 10${}^{-2}$	1.1 × 10${}^{-2}$	1.7 × 10${}^{-3}$	8.5 × 10${}^{-5}$	2.2 × 10${}^{-6}$	8.1 × 10${}^{-11}$	5.4 × 10${}^{-13}$
$\sigma$ = 0.35	4.3 × 10${}^{-1}$	3.6 × 10${}^{-2}$	5.2 × 10${}^{-3}$	4.9 × 10${}^{-4}$	1.8 × 10${}^{-15}$	1.1 × 10${}^{-15}$	-	-
$\sigma$ = 0.70	4.3 × 10${}^{-1}$	1.8 × 10${}^{-2}$	1.1 × 10${}^{-15}$	2.0 × 10${}^{-15}$	-	-	-	-

shows the behavior of the MAE versus the M in ${\phi }_m\left(t\right)=t^{{\delta }_m},$${\delta }_m=\sigma \left(m-1\right)$. As shown in the previous section (see Equations (97), (98) and (107)), the approximate solution $w\left(t,M\right)$ belongs to the linear span $S\left(\alpha ,\sigma ,M\right)=Span{\left(1,t,t^{\alpha +\sigma \left(m-1\right)}\right)}_{m=1}^M$. It is easy to check that for $\alpha =1.6$, $\sigma =0.25$, there is no such M that the exact $w\left(t\right)$ (see Equation (167)) belongs to $S\left(\alpha ,\sigma ,M\right)$. On the other hand, $w\left(t\right)\in S\left(1.6,0.35,M\ge 5\right)$ and $w\left(t\right)\in S\left(1.6,0.7,M\ge 3\right)$. The data placed in Table 5 illustrate this situation. For $\sigma =0.25$, the error decreases step-by-step with the growth of M, while it decreases sharply for $\sigma =0.35$ and $\sigma =0.7$ when the exact solution belongs to the corresponding manifold $S\left(\alpha ,\sigma ,M\right)$. Feng, Liu, and Turner have considered this problem [32] and constructed two finite element schemes for its numerical solution. Ezz-Eldien et al. [55] have studied this problem by the use of the combination of the shifted Legendre polynomials with the time-space spectral collocation method. The comparison of the two methods presented in Table 4 of [55] shows that most accurate result that has been achieved there has the error $E_{max}=1.407\times {10}^{-3}$ for the first technique and $E_{max}=1.807\times {10}^{-6}$ for the second one.

6.2.2. Example 5

In this example, the problem solved here is to show the applicability of the proposed algorithm for the multi-term time-fractional diffusion-wave equation in the unit square ${\left[0,1\right]}^2$ [56]

$$D_t^{\left(\alpha \right)}\left[u\right]+\frac{\partial u}{\partial t}+D_t^{\left({\alpha }_1\right)}\left[u\right]+u={\nabla }^{2}u+f\left(x,y,t\right).$$

(169)

The initial conditions, the Dirichlet BC, and the source term f correspond to the exact solution

$$u\left(x,y,t\right)=t^{2}sin\left(1-x\right)\left(e^x-1\right)sin\left(1-y\right)\left(e^y-1\right).$$

(170)

The data shown in

Table 6. The MAE and the elapsed time versus the N at $t=1$ using $M=5$ (left) and 10 (right).

	$\mathit{M}=5$				$\mathit{M}=10$
$\mathit{N}$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}\left(\mathit{u}\right)$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}\left(\partial \mathit{u}/\partial \mathit{t}\right)$		CPU, s	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}\left(\mathit{u}\right)$	${\mathit{E}}_{\mathbf{m}\mathbf{a}\mathbf{x}}\left(\partial \mathit{u}/\partial \mathit{t}\right)$		CPU, s
100	1.443 × 10${}^{-6}$	3.105 × 10${}^{-6}$		9	1.452 × 10${}^{-6}$	2.903 × 10${}^{-6}$		14
200	2.018 × 10${}^{-7}$	1.553 × 10${}^{-6}$		23	2.100 × 10${}^{-7}$	4.200 × 10${}^{-7}$		35
300	1.203 × 10${}^{-7}$	1.522 × 10${}^{-6}$		56	6.445 × 10${}^{-8}$	1.289 × 10${}^{-7}$		97
400	1.231 × 10${}^{-7}$	1.514 × 10${}^{-6}$		89	2.756 × 10${}^{-8}$	5.510 × 10${}^{-8}$		154
500	1.241 × 10${}^{-7}$	1.511 × 10${}^{-6}$		153	1.420 × 10${}^{-8}$	2.839 × 10${}^{-8}$		260
600	1.245 × 10${}^{-7}$	1.510 × 10${}^{-6}$		204	8.247 × 10${}^{-9}$	1.649 × 10${}^{-8}$		352
[56], Table 1, $E_{max}\left(u\right)$ = 6.145893072032060 × 10${}^{-6}$

are obtained using $\alpha =1.3$, ${\alpha }_1=0.3$, $\sigma =0.25$, and the numbers of the Müntz polynomials are $M=5$ and $M=10$. The same problem was considered by Shen, Liu, and Anh in [56] using an implicit difference method. From this table, it is clearly stated that our new approach is generally more accurate than others, even with a small number of N and M.

6.3. $(3+1)$-Dimensional Problems

Example 6

Let us consider the following time-fractional telegraph equation in three dimensions

$$D_t^{\left(\alpha +1\right)}\left[u\right]+D_t^{\left(\alpha \right)}\left[u\right]+u={\nabla }^{2}u+f\left(x,y,z,t\right),0\le t\le 1,$$

(171)

in the domain $\left(x,y,z\right)\in {\left[0,\pi \right]}^3$ with zero Dirichlet boundary conditions and the source term corresponding to the exact solution

$$u\left(x,y,z,t\right)=\frac{t^{\alpha +2}}{\mathsf{\Gamma }\left(\alpha +3\right)}sin2xsin2ysin2z.$$

(172)

Using the transform $\left(x,y,z\right)\to \left(\pi x,\pi y,\pi z\right)$, the equation is transformed into

$$D_t^{\left(\alpha +1\right)}\left[u\right]+D_t^{\left(\alpha \right)}\left[u\right]+u=\frac{1}{{\pi }^2}{\nabla }^{2}u+f\left(x,y,z,t\right),$$

(173)

with the solution

$$u\left(x,y,z,t\right)=\frac{t^{\alpha +2}}{\mathsf{\Gamma }\left(\alpha +3\right)}sin2\pi xsin2\pi ysin2\pi z.$$

(174)

Thus, we get a single spatial harmonic TFPDE. Substituting

$$u\left(x,y,z,t\right)=w\left(t\right)sin2\pi xsin2\pi ysin2\pi z,$$

(175)

we get the FODE

$$D_t^{\left(\alpha +1\right)}\left[w\right]+D_t^{\left(\alpha \right)}\left[w\right]+13w=F\left(t\right),0\le t\le 1,$$

(176)

with the source term and ICs corresponding to the exact solution

$$w\left(t\right)=\frac{t^{\alpha +2}}{\mathsf{\Gamma }\left(\alpha +3\right)}.$$

(177)

Table 7. The MAE versus M at $t=1$ for different $\sigma$.

M	1	2	3	4	5	6	10	15
$\sigma$ = 0.23	4.6 × 10${}^{-2}$	2.4 × 10${}^{-2}$	3.4 × 10${}^{-3}$	1.2 × 10${}^{-4}$	1.6 × 10${}^{-6}$	1.2 × 10${}^{-7}$	2.4 × 10${}^{-10}$	3.9 × 10${}^{-12}$
$\sigma$ = 0.25	4.6 × 10${}^{-2}$	2.3 × 10${}^{-2}$	3.0 × 10${}^{-3}$	8.2 × 10${}^{-5}$	1.1 × 10${}^{-17}$	2.8 × 10${}^{-17}$	-	-

shows the behavior of the absolute errors with the growth of M for two cases: $\sigma =0.23$ and $\sigma =0.25$. The results tabulated in the table are obtained by using $\alpha =0.9.$ As shown above (see Equations (97), (98) and (107)), the approximate solution $w\left(t,M\right)$ belongs to the linear span $S\left(\alpha ,\sigma ,M\right)=Span{\left(1,t,t^{\alpha +\sigma \left(m-1\right)}\right)}_{m=1}^M$. It is easy to check that for $\sigma =0.23$ there is no such M that the exact solution $w\left(t\right)$ belongs to $S\left(\alpha ,\sigma ,M\right)$. On the other hand, $w\left(t\right)\in S\left(\alpha ,0.25,M\ge 5\right)$. Indeed, $\alpha +0.25\times \left(5-1\right)=\alpha +2.$ The data placed in the table illustrate this situation. For $\sigma =0.23$, the error decreases step-by-step with the growth of M, while it decreases sharply for $\sigma =0.25$ when the exact solution belongs to the corresponding manifold $S\left(\alpha ,\sigma ,M\right)$. Yang et al. have considered this problem in [57] using an ADI finite difference scheme. The most accurate result shown in Table 2 of the reference has the error $1.1243\times {10}^{-3}$.

7. Conclusions

In the present work, an accurate method that can reach the computer rounding errors has been proposed for solving multi-term time-fractional equations. Let us note that the proposed analytical–numerical method collides with two key issues. The first one is the method to handle non-homogeneous time-dependent boundary conditions, which is critical to the derivation of $v_g$. The second problem is the method to handle the non-homogeneous time-dependent source term of the equation. The derived function $v_g$ solves the first problem cardinally. This article presents the analytical function $v_g$ exactly satisfying the boundary conditions. Let us note that this function is not unique because one can locally change it inside the cube. The algorithm can give the function in the explicit analytical form with the help of math software packages like Maple or Mathematica if needed. The method of the Laplace transform and the Green function method are traditionally used for solving the second one. The BSM technique can also handle the non-homogeneous time-dependent source term. As mentioned above, in the solution procedure of the present method, one gets the system of the equation for each Fourier harmonic. The FODEs were solved independently with the help of Müntz polynomial bases. Additionally, the Fourier expansion of the source term over the eigenfunction basis should be performed at several fixed time moments $t_j$ only. The number of these points is proportional to M. The value of M is not too large. As the numerical experiment shows, even $M=10$ for Müntz’s basis functions are enough for a rather precise approximate solution of the $(3+1)$ multi-term FPDE of the high order. And the convergence order with respect to the M and N is larger than 3. Generally speaking, $M>10$ and $N>100$ are sufficient to provide accurate results for the tested problems.

In this paper, the method is demonstrated by solving an important class of fractional problems described in the Introduction. This technique can also be extended onto iterative quasi-linear PDEs and onto the problems in other orthogonal coordinate systems. The main drawback of this work is that the derivation of $v_g$ is only applicable for regular domains. Actually, numerical methods can be used for this purpose. This is the subject of further study.