images
Computer Modeling in Engineering & Sciences
images

DOI: 10.32604/cmes.2023.030449

ARTICLE

A Novel Accurate Method for Multi-Term Time-Fractional Nonlinear Diffusion Equations in Arbitrary Domains

Tao Hu1, Cheng Huang2, Sergiy Reutskiy3,*, Jun Lu4 and Ji Lin5,*

1Ningbo Yongxin Engineering Consulting Co., Ltd., Ningbo, 315042, China
2Taizhou Huanfeng Water Conservancy Engineering Co., Ltd., Taizhou, 318050, China
3A. Pidhornyi Institute of Mechanical Engineering Problems, National Academy of Sciences of Ukraine, Kharkiv, 61046, Ukraine
4Material and Structure Department, Nanjing Hydraulic Research Institute, Nanjing, 210029, China
5College of Mechanics and Materials, Hohai University, Nanjing, 211100, China
*Corresponding Authors: Sergiy Reutskiy. Email: sergiyreutskiy@gmail.com; Ji Lin. Email: linji@hhu.edu.cn
Received: 06 April 2023; Accepted: 20 June 2023; Published: 17 November 2023

Abstract: A novel accurate method is proposed to solve a broad variety of linear and nonlinear (1+1)-dimensional and (2+1)-dimensional multi-term time-fractional partial differential equations with spatial operators of anisotropic diffusivity. For (1+1)-dimensional problems, analytical solutions that satisfy the boundary requirements are derived. Such solutions are numerically calculated using the trigonometric basis approximation for (2+1)-dimensional problems. With the aid of these analytical or numerical approximations, the original problems can be converted into the fractional ordinary differential equations, and solutions to the fractional ordinary differential equations are approximated by modified radial basis functions with time-dependent coefficients. An efficient backward substitution strategy that was previously provided for a single fractional ordinary differential equation is then used to solve the corresponding systems. The straightforward quasilinearization technique is applied to handle nonlinear issues. Numerical experiments demonstrate the suggested algorithm’s superior accuracy and efficiency.

Keywords: Müntz polynomial basis; backward substitution method; collocation method; meshless method; fractional equation

1  Introduction

In this study, we consider an efficient and accurate numerical technique for the problem of the following form:

ℒt[v]=ℳt[L(x)[v]]+N(v,vx,vy)+h(x,t), t≥0, x∈Ω⊂R2,(1)

where

L(x)[v]=div(D^(x)∇v)=d1,1(x)∂2v∂x12+2d1,2(x)∂2v∂x1∂x2+d2,2(x)∂2v∂x22+(∂d1,1(x)∂x1+∂d1,2(x)∂x2)∂v∂x1+(∂d1,2(x)∂x1+∂d2,2(x)∂x2)∂v∂x2,(2)

is the operator of the anisotropic diffusivity. d1,1(x), d1,2(x), and d2,2(x) are the elements of the anisotropic matrix D^. ℒt and ℳt are given as follows:

ℒt=Dt(μ)+∑k=1Iak(t)Dt(μk), ℳt=∑k=I+1Kak(t)Dt(μk),(3)

which are time differential operators of integer or fractional orders. And N(u,ux,uy) is a nonlinear function in the unknown solution and its derivatives. The values 0≤μk<μ≤1 are fractional constant values and ak(t), k=1,...,K are real smooth time-varying functions available in advance. The operator Dt(ν) is the Caputo-type fractional derivative

Dt(ν)[f(x,t)]={1Γ(n−ν)∫0t∂t(n)f(x,τ)dτ(t−τ)ν−n+1,n−1<ν<n,∂t(n)f(x,t),ν=n.(4)

Eq. (1) can be used to describe a wide variety of mathematical models in engineering and science as particular cases. First of all, these are the models of the heat equations obeying the classical Fourier law. At the same time, there are many models of heat exchange based on the non-Fourier laws which also can be represented by this equation. Besides, it can also be used for simulating the unsteady flow of the non-Newtonian fractional Maxwell fluid [1] and the Oldroyd-B fluid [2]. The general equation given in Eq. (1) also includes the Benjamin-Bona-Mahony and the Benjamin-Bona-Mahony-Burgers equations [3,4], which have been used for modeling long waves of small amplitudes. Besides, the Sobolev equation also falls into this group [5,6].

So, the multi-term systems of equations are introduced for modeling many real-life applications which involved some complicated processes that cannot be accurately described by systems of the single-term. It is evident that only a small number of problems can be analyzed analytically under some idealized conditions which are useful for parametric research. However, the derivation of analytical solutions for the nonlinear fractional equations with general spatial differential operators is not a trivial task. Applying numerical methods to deal with these equations would be more desirable. The finite difference method (FDM) [7–9] and the finite element method (FEM) [10–12] have been proposed widely as available techniques. Jin et al. [13] investigated the multi-term time-fractional diffusion equations using the Galerkin FEM. Recently an alternating direction implicit (ADI) scheme has been presented by Huang et al. in [14]. The weighted meshless spectral method was proposed by Hussain and Haq for problems arising in mass and heat transfer [15]. A method based on the Haar wavelets and the finite difference scheme was proposed by Oruç, Esen, and Bulut for two-dimensional time fractional reaction-sub-diffusion equation in [16]. An accurate computational method based on the linear barycentric interpolation method has been proposed by Oruç [17] recently for solving fractional Rayleigh-Stokes problems. The numerical method which combines the Laplace transform with the local radial basis functions method was presented by Li et al. in [18] for fractional diffusion models. A Galerkin method based on the second kind Chebyshev wavelets was presented by Soltani Sarvestani et al. in [19] for solving the multi-term time fractional diffusion-wave equation. Very recently an alternating direction implicit Legendre spectral method has been developed by Liu et al. in [20]. A mixed meshless method was proposed by Nikan et al. in [21]. The weak Galerkin finite element method for the space discretization combined with the discretization of the Caputo time derivative by L1 method was proposed by Zhou et al. in [22] for time-fractional quasi-linear diffusion equation. A method that combines the L1 interpolation of the fractional derivative with a high-order compact finite volume scheme was applied by Su et al. in [23] for the 2D multi-term time fractional sub-diffusion equation. The technique which combines approximation of the time-fractional derivative via the L1 formula and approximation of the integer order space derivatives by truncated one and two-dimensional wavelet series was presented by Ghafoor et al. for one and two dimensional higher order multi-term time-fractional partial differential equations in [24].

Inspired by the advantages of meshless methods which do not need the requirement of domain mesh, a variety of meshless methods have been proposed in the literature for solving fractional equations in both theoretical and numerical aspects. There exist several types of meshless methods, from which the radial basis function (RBF)-based meshless methods are the most popular type. With the merits of Euclidian distances as variables, RBF-based methods are flexible and usual for high-dimension problems under irregular complicated domains. Piret et al. [25] made the first attempt to use the RBF discretization for the fractional diffusion. Hosseini et al. [26] solved the time-fractional telegraph equations using the RBF discretization with the finite difference method for temporal discretization. Liu et al. [27] provided the implicit RBF approach with error analysis. The radial basis function-finite difference method (RBF-FD) [28] was also used to solve the fractional equations. In this paper, an accurate and efficient method based on the key idea of the backward substitution method has been proposed [29–34]. In this method, the general approximation based on the pure radial basis functions and their correcting terms are formed. Applying the collocation method, we can reduce the considered problems into the system of fractional ordinary differential equations. Then the corresponding problems are solved efficiently with the help of the Müntz polynomial basis since the fractional derivative of a Müntz polynomial is again a Müntz polynomial.

The remainder is organized as follows. These preliminaries including the definition of fractional derivative, the backward substitution method for the linear system of fractional ordinary differential equations, and the quasilinearization techniques are described in Section 2. The main method of application to (1+1)-dimensional and (2+1)-dimensional problems are described in Section 3 where the description of the method in application to the problems with the spatial operator of the fourth order has also been discussed. The numerical examples which illustrate the method presented are placed with the comparisons by some well-known numerical methods in Section 4. Finally, in Section 5, a short conclusion and discussion is provided.

2  Preliminaries

2.1 Solutions for Linear System of Fractional Ordinary Differential Equations

The method presented in this study is based on the effective algorithm for the solution of linear systems of the fractional ordinary differential equations (FODEs). This algorithm is described detailed in this subsection. Let us consider the linear system of fractional ordinary differential equation

A^Dt(μ)[P(t)]=∑k=1KA^k(t)Dt(μk)[P(t)]+B^(t)P(t)+F(t), 0≤t≤T,(5)

with zero initial conditions

P(0)=0,(6)

where P(t)=[P1(t),P2(t),...,PN(t)]T,A^ is a constant non-singular matrix of N×N dimension, A^k(t) and B^(t) are time dependent matrices of the same size, F(t)=[f1(t),f2(t),...,fN(t)]T. In order to solve the fractional equation easily, we introduce the Müntz polynomial basis (MPB)

φm(t)=tδm, δm=σ(m−1), 0<σ≤1,m=1,2,3,...,(7)

as the basis function where the right hand side of Eq. (5) can be approximated by the linear combination of these functions

∑k=1KA^k(t)Dt(μk)[P(t)]+B^(t)P(t)+F(t)=A^∑m=1∞qmφm(t), 0≤t≤T,(8)

where qm are unknown coefficients to be determined which will be discussed in the following section.

Some properties of Müntz polynomials should be discussed and noted here. Let {δ0,δ1,δ2,...} be a sequence of distinct real numbers such that 0≤δ0<δ1<δ2<...→∞. The Müntz polynomials of the form ∑k=0ncktδk, with real coefficients c0,c1,...,cn are dense in L2[0,1] if and only if

∑k=1∞1δk=+∞.(9)

Also if δ0=0, then the Müntz polynomials are dense in C[0,1], with the uniform norm, if and only if Eq. (9) is established [35,36]. With the help of Eqs. (8) and (5), we have

A^Dt(μ)[P(t)]=A^∑m=1∞qmφm(t).(10)

Followed from Eq. (10), we have the following equation:

Dt(μ)[P(t)]=∑m=1∞qmφm(t),(11)

if the matrix A^ in Eq. (10) is a non-singular matrix. As it is easily to prove that the following analytical expression:

Φm(t)=Γ(δm+1)Γ(δm+μ+1)tδm+μ,(12)

satisfies

Dt(μ)[Φm(t)]=φm(t),Φm(0)=0.(13)

Taking into account of Eq. (6), the following series:

P(t)=∑m=1∞qmΦm(t),(14)

is the analytical solution to Eq. (11) and the unknown vectors qm can be solved by

∑k=1KA^k(t)Dt(μk)[∑m=1∞qmΦm(t)]+B^(t)∑m=1∞qmΦm(t)+F(t)=A^∑m=1∞qmφm(t),(15)

or

∑m=1∞[A^φm(t)−B^(t)Φm(t)−∑k=1KA^k(t)Φm(μk)(t)]qm=F(t), t∈[0,T],(16)

here

Φm(μk)(t)=Dt(μk)[Φm(t)]=Dt(μk)[Γ(δm+1)Γ(δm+μ+1)tδm+μ]=Γ(δm+1)tδm+μ−μkΓ(δm+μ+1−μk).(17)

It is worth noting that if Eq. (16) is fulfilled for any t in the time sequence, the infinite sequence given in Eq. (14) is the desired solution. For the sake of calculation, we consider the truncated series

PM(t)=∑m=1MqmΦm(t),(18)

which satisfies the corresponding equation of truncated series as follows:

Dt(α)[PM(t)]=∑m=1Mqmφm(t),(19)

where the unknown parameters can be obtained by enforcing the Eq. (16) at several selected time steps tj. To ensure the solvability of the corresponding system, the number of equations (i.e., the number of time steps) should be taken larger than the number of unknown variables.

2.2 Quasilinearization Procedure for Nonlinear Problem

The nonlinear term N(v,vx,vy) in the Eq. (1) can be transformed into linear terms using the quasilinearization technique [37] which is briefly illustrated as follows.

Let us denote ξi=∂v∂xi and assume that v0(x) and ξi,0(x) are known functions which are used as initial approximations of the exact values. So, we can write

v=v0+(v−v0)=v0+Δv, ξi=ξi,0+(ξi−ξi,0)=ξi,0+Δξi.(20)

Assuming that Δv and Δξi are small enough, then we can neglect the squares of their values, in this way

N(v,vx,vy,x,t)≈N(v0,ξ1,0,ξ2,0,x,t)+∂N(v0,ξ1,0,ξ2,0,x,t)∂v(v−v0)+∑i=12∂N(v0,ξ1,0,ξ2,0,x,t)∂ξi(ξi−ξi,0).(21)

As a result the original equation is transformed to a sequence of linear equations with the initial approximations of the form

ℒt[v]=ℳt[L(x)[v]]+b1(x,t)vx+b2(x,t)vy+c(x,t)v+f(x,t),(22)

and the system of iteration will stop under the control of the error given by two successive evaluations maxx,t|v(x,t)−v0(x,t)| which is smaller than the selected tolerance. In practical calculations, we can simply fix the number of iterations and from the numerical examples, the iteration will converge in only several steps.

3  The Solution Scheme for the Considered Problem

In this section, the solution scheme for the Eq. (22) with the known boundary conditions (BC) and initial conditions (IC) will be discussed in details, as follows:

ℬ[v(x,t)]=g(x,t), x∈∂Ω,(23)

v(x,0)=v0(x),(24)

in which ℬ is the boundary condition operator. In this case, we consider the first kind and the third kind boundary conditions where the second kind boundary condition can be easily obtained, as follows:

u=g1(x,t), x∈Γ1,(25)

r(x)u+∂u∂n=r(x)u+n1∂u∂x1+n2∂u∂x2=g2(x,t), x∈Γ2,(26)

in which n=(n1,n2) represents the unit outward normal vector to the physical boundary and r(x), g1(x,t), and g2(x,t) given in Eqs. (25) and (26) are known in advances with sufficient smoothness. Note that the boundary condition of the second kind (Neumann BC) is given if r(x)≡0 in Eq. (26).

3.1 Algorithm for (1+1)-Dimensional Problem

In this case the governing equation takes the form

ℒt[v]=ℳt[L(x)[v]]+b(x,t)vx+c(x,t)v+f(x,t),(27)

or in the explicit form

ℒt[v]=Dt(μ)[v]+∑k=1Iak(t)Dt(μk)[v]=∑k=I+1Kak(t)[∂x(d(x)∂xDt(μk)[v])]+b(x,t)∂xv+c(x,t)v+f(x,t)=ℳt[L(x)[v]]+b(x,t)∂xv+c(x,t)v+f(x,t),(28)

with the IC

v(x,0)=v0(x),(29)

and the BCs

LW[v(x,t)]≡αW∂xv(0,t)+βWv(0,t)=gW(t),(30)

LE[v(x,t)]≡αE∂xv(1,t)+βEv(1,t)=gE(t),(31)

at the endpoints of the interval Ω=[0,1]. Let us define the function

u(x,t)=v(x,t)−v0(x),(32)

which satisfies the equation

ℒt[u]=ℳt[L(x)[u]]+b(x,t)∂xu+c(x,t)u+f1(x,t),(33)

the BCs

LW[u(x,t)]=gW(t)−LW[v0(x)]=hW(t),(34)

LE[u(x,t)]=gE(t)−LE[v0(x)]=hE(t),(35)

and zero IC

u(x,0)=0.(36)

Here

f1(x,t)=f(x,t)+b(x,t)∂xv0(x)+c(x,t)v0(x).(37)

Let us define the following two functions where the parameters are determined from the initial condition and boundary condition:

θE(x)=αW−βWxαWβE−βW(αE+βE),(38)

θW(x)=βEx−(αE+βE)αWβE−βW(αE+βE),(39)

satisfying the conditions

LW[θW(x)]=1, LE[θW(x)]=0,(40)

LW[θE(x)]=0, LE[θE(x)]=1.(41)

Then, it is easy to prove that the following combination:

ug(x,t)=θE(x)hE(t)+θW(x)hW(t),(42)

satisfies the boundary conditions

LW[ug(x,t)]=hW(t), LE[ug(x,t)]=hE(t).(43)

Suppose that the ug(x,t) provided in Eq. (42) is the initial approximation to the solution, then by using the decomposition method, we have the equation for the correction item w(x,t), as follows:

u(x,t)=ug(x,t)+w(x,t),(44)

where w(x,t) should satisfy the following governing equation:

ℒt[w]=ℳt[L(x)[w]]+b(x,t)∂xw+c(x,t)w+f2(x,t),(45)

where

f2(x,t)=f1(x,t)−ℒt[ug]+ℳt[L(x)[ug]]+b(x,t)∂xug+c(x,t)ug,(46)

and the following initial condition and boundary condition:

w(x,0)=0,(47)

LW[w(x,t)]=0, LE[w(x,t)]=0.(48)

In order to approximate w(x,t), we consider the radial basis function ψi(x) which has been widely used in meshless methods

ψi(x)=(x−ζi)2+c2,(49)

where c is the artificial selection shape parameter and ζi are the centers of the considered function distributed in the solution domain. We define the modified basis functions as follows:

ϕi(x)=ψi(x)+ci,0+ci,1x,(50)

where coefficients ci,0 and ci,1 are chosen in the way that the following equations are fulfilled:

LW[ϕi(x)]=LE[ϕi(x)]=0.(51)

As a result we get the linear system

αWci,1+βWci,0=−αW∂xψi(x)x=0−βWψi(0),(52)

αEci,1+βE(ci,0+ci,1)=−αE∂xψi(x)x=1−βEψi(1),(53)

for each pair of the coefficients ci,0 and ci,1. The system can be solved easily in the analytical form. By considering the properties of the modified functions ϕi(x) which satisfying the homogeneous boundary conditions, the linear combination of ϕi(x) is used as the approximation to wN(x,t), as follows:

wN(x,t)=∑i=1Nϕi(x)Pi(t),(54)

in which the unknown Pi(t) have to be determined. Substituting Eq. (54) into Eq. (45), we can yield the following system

∑i=1Nϕi(x)ℒt[Pi(t)]=∑i=1NL(x)[ϕi(x)]ℳt[Pi(t)]+∑i=1N[b(x,t)∂xϕi(x)+c(x,t)ϕi(x)]Pi(t)+f2(x,t).(55)

Let xj∈[0,1], j=1,…,N be selected nodes in the solution domain, the following system of FODEs will be obtained after considering collocation procedure:

∑i=1Nϕi(xj)ℒt[Pi(t)]=∑i=1NL(x)[ϕi(xj)]ℳt[Pi(t)]+∑i=1N[b(xj,t)∂xϕi(xj)+c(xj,t)ϕi(xj)]Pi(t)+f2(xj,t),(56)

or simplified in the matrix form as shown below:

A^Dt(μ)[P(t)]+∑k=1Iak(t)A^Dt(μk)[P(t)]=∑k=I+1Kak(t)DkDt(μk)[P(t)]+B^(t)P(t)+F(t), 0≤t≤T.(57)

Here the derivatives can be obtained in the analytical form

∂xϕi(x)=x−ζiψi(x)+ci,1, ∂xxϕi(x)=1ψi(x)−(x−ζi)2(ψi(x))3.(58)

A^, D^k, and B^(t) are N×N matrices with the components

A^=(aj,i)j,i=1N=(ϕi(xj))j,i=1N, (59)

B^(t)=(b(xj,t)∂xϕi(xj)+c(xj,t)ϕi(xj))j,i=1N,(60)

D^k=(∂∂x(d(xj)∂xϕi(xj)))j,i=1N(61)

and P(t)=[P1(t),P2(t),...,PN(t)]T, F(t)=[f2(x1,t),f2(x2,t),...,f2(xN,t)]T are N-vectors. If we denote

ak(t)A^=−A^k(t), k=1,...,I; ak(t)D^k=A^k(t), k=I+1,...,K,(62)

then the system Eq. (57) coincides with the linear system of FODEs shown in Eq. (5) and can be solved using the algorithm described there. Using M Müntz polynomials in solving of the system of the FODEs shown in Eq. (5), we get the approximate solution

wN,M(x,t)=∑i=1Nϕi(x)PM,i(t),(63)

with the vector PM(t) given in Eq. (18). Then, the approximate solution of the original problem Eqs. (27), (29)–(31) can be given in the following way:

vN,M(x,t)=uN,M(x,t)+v0(x)=wN,M(x,t)+ug(x,t)+v0(x)=∑i=1Nϕi(x)PM,i(t)+ug(x,t)+v0(x).(64)

The same algorithm can be applied for solving problems of Eq. (27) with the spatial operator of the fourth order

L(x)[v]=−∂xx(r(x)∂xxv)+∂x(q(x)∂xv)−p(x)v.(65)

In this case, the modified RBF takes the form

ϕi(x)=ψi(x)+ci,0+ci,1x+ci,2x2+ci,3x3,(66)

where the coefficients ci,0, ci,1, ci,2, ci,3 are chosen to satisfy the homogeneous BCs. For example, let us consider the BCs:

ϕi(0)=ϕi(1)=∂xxϕi(0)=∂xxϕi(1)=0,(67)

which leads to the four conditions

1) ci,0=−ψi(0), 2) ci,0+ci,1+ci,2+ci,3=−ψi(1), 3) 2ci,2=−∂xxϕi(0),4) 2ci,2+6ci,3=−∂xxϕi(1).(68)

So, from 3) ⇒ci,2=−1/2∂xxϕi(0), then from 4) ⇒ci,3=(∂xxϕi(0)−∂xxϕi(1))/6 and finally, from 2) ci,1=−ψi(1)−ci,0−ci,2−ci,3. The derivatives of the modified RBF can be calculated in explicit analytical form:

∂xϕi(x)=x−ζiψi(x)+ci,1+2ci,2x+3ci,3x2, ∂xxϕi(x)=1ψi(x)−(x−ζi)2ψi3(x)+2ci,2+6ci,3x,∂xxxϕi(x)=−3(x−ζi)ψi3(x)+3(x−ζi)3ψi5(x)+6ci,3, ∂xxxxϕi(x)=−3ψi3(x)+18(x−ζi)2ψi3(x)−15(x−ζi)4ψi7(x).(69)

3.2 Algorithm for (2+1)-Dimensional Problem

Following the decomposition for the (1+1)-dimensional problem, we have

u(x,t)=v(x,t)−v(x,0)≡v(x,t)−v0(x),(70)

which can transform problems of Eq. (22) into the form

ℒt[u]=ℳt[L(x)[u]]+b1(x,t)∂xu(x,t)+b2(x,t)∂yu(x,t)+c(x,t)u(x,t)+f1(x,t),(71)

where

 f1(x,t)=f(x,t)+b1(x,t)∂xv0(x)+b2(x,t)∂yv0(x)+c(x,t)v0(x).(72)

The solution u(x,t) satisfies homogeneous initial condition:

u(x,0)=0,(73)

and the updated boundary conditions, as follows:

ℬ[u(x,t)]=g(x,t)−ℬ[v(x,0)], x∈∂Ω,(74)

where the operator in the last equation is described by the formulae of Eqs. (25) and (26). To transform the time-dependent fractional equations into a system of fractional ordinary equations and apply the algorithm described in the previous subsection for (1+1)-dimensional problems, we should perform the following steps for this goal:

1) We transform Eq. (71) to the homogeneous BCs by introducing the primary approximation ug:

u(x,t)=ug(x,t)+w(x,t),(75)

where ug(x,t) is a man-made function which has to satisfy

ℬ[ug(x,t)]=g(x,t), x∈∂Ω.(76)

As discussed in the (1+1)-dimensional problems, by considering the Eq. (75) and the governing equation, we have the equation to determine the correcting solution w, as follows:

ℒt[w]=ℳt[L(x)[w]]+b1(x,t)∂xw(x,t)+b2(x,t)∂yw(x,t)+c(x,t)w(x,t)+f2(x,t),(77)

where the new source term is

f2(x,t)=f1(x,t)−ℒt[ug(x,t)]+ℳt[L(x)[ug(x,t)]]+b1(x,t)∂xug(x,t)+b2(x,t)∂yug(x,t)+c(x,t)ug(x,t),(78)

which can now be solved using the method as described in the (1+1)-dimensional problems. The main difference for the solution process of the (1+1)-dimensional problems and the (2+1)-dimensional problems is that the initial approximation of the solution ug and the correcting function ωi(x) for the (1+1)-dimensional problem can be obtained analytically. However, for the (2+1)-dimensional problems, these functions can only be obtained numerically. In the following subsection, the goal for obtaining such functions are detailed.

3.3 The Procedure for Obtaining ug(x,t) and ωi(x)

Since the artificially designed ug(x,t) and ωi(x) do not have the requirements to satisfy the governing equations, we apply the trigonometric basis, as follows:

θk(β,x)=sin⁡(kπx+β2β),(79)

θk1,k2(β,x)=sin⁡(k1πx1+β2β)sin⁡(k2πx2+β2β),(80)

for the (1+1)-dimensional problems and the (2+1)-dimensional problems, respectively, where the parameter β is selected in the way that Ω⊂Ωβ. Then we use the following approximations of the linear combination of θ with weighted parameters

ωi(x)=∑k=1𝒦pi,kθk(β,x), i=1,...,N, (81)

ug(x,t)=∑k=1𝒦sk(t)θk(β,x),(82)

where the weighted parameters are determined using the collocation approach

∑k=1𝒦pi,kℬ[θk(β,ζl)]=Hi, ζl∈∂Ω, i=1,…,N, l=1,…,𝒦c,(83)

∑k=1𝒦sk(t)ℬ[θk(β,ζl)]=G(t), ζl∈∂Ω, l=1,…,𝒦c,(84)

and the Hi and G(t) are given as follows:

Hi=−[B[ψi(ζ1)],B[ψi(ζ2)],...,B[ψi(ζ𝒦c)]]T, G(t)=[g(ζ1,t),g(ζ2,t),...,g(ζ𝒦c,t)]T.(85)

It is important to note that we need approximation of the function ug in the way which admits calculation the fractional derivatives (see Eq. (78)). When the sought solution has the particular form u(x,t)=u0(x)ς(t), then a simplified procedure can be applied: we find ug,0(x) as a boundary approximation of u0(x) and define ug(x,t)=ug,0(x)ς(t). Then the derivatives can be found as Dt(α)[ug(x,t)]=ug,0(x)Dt(α)[ς(t)] at each time moment t. In general case, this approach is impossible.

To make this possible, we use the point interpolation method (PIM) [38,39]. Let 0≤t1<t2<...<tI≤T be the collocation points in the time interval. Let Sk=[sk(t1),sk(t2),...,sk(tI)] be the vector of the values of the smooth function sk(t) calculated at these points. We seek approximation of sk(t) at a point of interest t in the form of

s~k(t,I)=∑i=1Iti−1ak,i,=[1,t,...,tI−1]⏟pT{ak,1...ak,I}=pTak,(86)

and the coefficients ak,1 in Eq. (86) can be determined by enforcing s~k(t,I) to pass through the values at the collocation points tj. This yields I equations for each collocation point, i.e.,

sk(t1)=ak,1+ak,2t1+...+ak,It1I−1sk(t2)=ak,1+ak,2t2+...+ak,It2I−1...sk(tI)=ak,1+ak,2tMp+...+ak,ItII−1}⇒Sk=𝒫ak,(87)

where

𝒫^=[1,t1,...,t1I−11,t2,...,t2I−1...1,tMp,...,tII−1],(88)

is the so-called square moment matrix. Assuming that 𝒫 is a non-singular matrix, we get

ak=𝒫^−1Sk.(89)

Using vectors ak, k=1,...,𝒦, we obtain ug(x,t) in the form of the smooth function differentiable in time.

ug(x,t)=∑k=1𝒦∑i=1Iti−1ak,iθk(β,x)=∑i=1I(∑k=1𝒦ak,iθk(β,x))ti−1=∑i=1IΘi(β,x)ti−1.(90)

As a result, we get the approximate solution vN,M(x,t) in the form similar to the one obtained in the case of the (1+1) dimensional problem

vN,M(x,t)=∑i=1Nϕi(x)PM,i(t)+ug(x,t)+v0(x),(91)

where the approximate solution also depends on the parameters 𝒦 and I of the function ug(x,t).

4  Numerical Examples

In this section, the accuracy and efficiency of the proposed method is verified. The maximum absolute error (MAE, Emax(t)), the relative root mean square error (Erel(t)), and also the error in discrete H1(Ω) norm (EH1(t)) are used as criteria. The detailed computational strategies are given:

Emax(t)=max1≤i≤Nt|uN,M(xi,t)−uex(xi,t)|,(92)

Erel(t)=∑i=1Nt[uex( xi,t)−uN,M(xi,t)]2/∑i=1Ntuex2(xi,t),(93)

EH1(t)=1Nt∑i=1Nt[(uex(xi,t)−uN,M(xi,t))2+∑k=12(∂uex∂xk(xi,t)−∂uN,M∂xk(xi,t))2],(94)

where uex(xi,t) and uN,M(xi,t) are analytical and numerical solutions, respectively, Nt is the number of the test points. To investigate the influence of N (the number of RBFs) and M (the number of Müntz polynomials) on the accuracy of the proposed scheme, we introduce the convergence orders (COs)

CON=log⁡(Emax(N1)/Emax(N2))log⁡(N1/N2), M is fixed,COM=log⁡(Emax(M1)/Emax(M2))log⁡(M1/M2), N is fixed.(95)

Furthermore, for the 1D problems, we let c=0.5 and for the 2D problems, we take c=2. The data provided in Examples 1 to 4 demonstrate that in the whole range σ∈[0.1,0.3], the Müntz basis provides a quite precise approximation in time. The calculation results are very close for all values of the parameter σ from this range. As the parameter σ increases, the error grows rapidly. It reaches a maximum at σ=1 when the usual polynomials of the integer order tm, m=0, 1, 2,... are used in approximation of time. This demonstrates the advantage of using the Müntz basis. Therefore, for the parameters in the Müntz polynomials, we have σ=0.3 for the rest examples and illustrations.

4.1 (1+1)-Dimensional Problems

4.1.1 Example 1

Let us consider the multi-term problem discussed in [13]

Dt(0.95)[u]+Dt(0.2)[u]=∂xxu+f(x,t), 0≤x,t≤1,(96)

with the first kind boundary condition and initial condition conform to u(x,t)=(1+t2)(x2−x).

Fig. 1 and Table 1 show the behaviour of the errors with respect to the number of the modified RBFs ϕi(x), N in the approximation of the solution (see Eq. (54)) with fixed M. For M=10, the predicted errors can be shapely reduced by increasing of the N for all the computed 2≤N≤26. This means that the error occurring in spatial approximation is the main error source in this case. While for M=5 the computed errors cannot be decreased for the selected N>10. Therefore, in this case, the dominant error is the error in the solution of the system. With the increasing of N, the minimum error of this example is 5.80E−10 and the convergence order vs. N is larger than 2. It is noticeable that for the M=10 and N=26, we can obtain the solution in 2.9 s which is sufficiently efficient. Fig. 1 shows that all the curves Emax(N), EH1(N) originally lie on the same curve. The predicted errors move away from this curve depending on the value of M, when the error occurring in the approximation in time plays an essential role.

Figure 1: The Emax (left) and EH1 (right) errors with the respect to the number of the RBFs used in the approximate solution

Fig. 2 and Table 2 show the computed Emax and EH1 errors vs. the number of Müntz polynomials M with fixed N. The same problem has already been considered by Jin et al. in [13] using the Galerkin finite element method and finite difference discretization of the time-fractional derivatives. In their practical computations, the size of the mesh is h=2−10 and the time step for the temporal approximation is τ=1/160. They obtained Emax=4.20E−4 which was placed in the Table 1 of the reference [13]. It is evident that the proposed method is far more accurate than the reference. The domain errors with N=26 and M=10 is given in Fig. 3 from which we can see that in the whole domain, the proposed method maintains the same error level for the test nodes which are close to the boundary.

Figure 2: The Emax (left) and EH1 (right) errors vs. the number of Müntz’s polynomials M

Figure 3: The domain absolute errors with N=26 and M=10

Finally, we show the behaviour of the errors Erel with the growth of N for different σ in Table 3. This table demonstrates that in the whole range σ∈[0.1,0.3], the Müntz basis provides a quite precise approximation in time.

4.1.2 Example 2

Let us consider the multi-term problem of the following form:

Dt(π/8)[u]+(1+t)Dt(0.25)[u]+(1+t2)Dt(0.125)[u]=sin⁡(t)Dt(0.2)[L(x)[u]]+cos⁡(t)Dt(0.1)[L(x)[u]]+etL(x)[u]+f(x,t), 0≤x,t≤1,(97)

with the boundary condition of third kind (Robin BCs)

∂xu(0,t)−πu(0,t)=g0(t), ∂xu(1,t)−eu(1,t)=g1(t).(98)

The BCs and IC conform the exact solution u(x,t)=exp⁡(x+t). Here the spatial operator is L(x)[u(x,t)]=∂x(cosh⁡(x)∂xu(x,t)).

Some results of the calculations are presented in Tables 4 and 5 and displayed in the Figs. 4 and 5. From Tables 4 and 5, we can see that for the problem of the Robin boundary conditions, the proposed method can yield accurate results with small number of basis functions and collocation nodes. And the overall convergence order of the propose scheme is large than two. The same conclusions can be made from Figs. 4 and 5, where the computed maximum absolute errors vs. the M and N are displayed. From these figures, we can see that with the increasing of N for the case N≥24, the errors adjust within a small range which may due to properties of the coefficient matrix such as the ill-conditioning. The following three ways can be used to solve this problem. In the first strategy, we can increase the number of M, for N≥24, then we can obtain more accurate and stable results from the over-determined system. The second attempt can be done to use the high-precision matrix solver to overcome the ill-conditioning of the matrix. Furthermore, the recently proposed multiple-scale method by Oruç [40] can also be used to reduce the condition number of the coefficient matrix. For the fact that with the increasing of M and the fixed N, the approximate solutions become more accurate and the errors only adjust within a small range, we will not discuss in details of these techniques for improving the solution stability.

Figure 4: The Emax (left) and EH1 (right) errors vs. N

Figure 5: The Emax (left) and EH1 (right) errors vs. M

In Table 6, we show the effect of the σ on the accuracy. The same conclusion can be made that in the whole range σ∈[0.1,0.3], the Müntz basis provides a quite precise approximation in time.

4.1.3 Example 3

Consider the nonlinear time-fractional Huxley equation of the following form:

Dt(0.5)[u(x,t)]=∂xxu(x,t)+u(x,t)(1−u(x,t))(u(x,t)−γ)+f(x,t), 0≤x,t≤1,(99)

with u(x,t)=x3(t+1).

Table 7 and Fig. 6 show the trends of the computed errors with the growth of N with fixed M. The predicted solutions are obtained after 3 iteration of the quasilinearization procedure. The initial approximations of the solution and its derivatives are provided randomly. From the results, we can verify that for the randomly provided initial approximations, we can obtained the accurate solutions in three iterations. And with the increasing of the N, the proposed algorithm converges fast. The same problem was considered by Hadhoud et al. in [41] using a numerical technique based on the cubic B-spline collocation method and the mean value theorem for integrals. The maximal absolute errors obtained there are near ∼10−6 for the mesh size Δx=0.01 and the time step Δt=0.01 (see Tables 4–6 of the original paper).

Figure 6: The maximal absolute Emax (left) and EH1 (right) errors as functions of the number of the RBFs used in the approximate solution

Now let us move onto the multi-term nonlinear time-fractional equation with the Huxley nonlinear term,

Dt(0.5)[u(x,t)]+sin⁡(t)Dt(0.3)[u(x,t)]+cos⁡(t)Dt(0.2)[u(x,t)]=sin⁡h(t)Dt(0.29)[L(x)[u(x,t)]]+cosh⁡(t)Dt(0.19)[L(x)[u(x,t)]]+etL(x)[u(x,t)]+u(x,t)(1−u(x,t))(u(x,t)−γ)+f(x,t), 0≤x,t≤1,(100)

with the same analytical solution. Here L(x)[u(x,t)]=∂x((1+x2)∂xu(x,t)).

Fig. 7 demonstrates the behaviour of the errors with the grows of the total amount of basis function ϕi(x) for the approximation of the solution (see Eq. (54)). The data correspond to the 3th iterations of the quasilinearization procedure with randomly provided initial approximations of the solution and its derivatives. If we apply much more steps of iterations, the solution accuracy do not change which means that three steps of iterations are sufficient for the proposed method.

Figure 7: The Emax (left) and EH1 (right) errors vs. the number of N

Then, we discuss the solution accuracy on the σ in Table 8. It is evident that the σ=0.3 is sufficient to do the approximation.

4.1.4 Example 4

Let us consider the nonlinear problem with the spatial differential operator of the fourth order

Dt(0.9)[u(x,t)]+tDt(0.1)[u(x,t)]+t2Dt(0.2)[u(x,t)]=11+tDt(0.29)[L(x)[u(x,t)]]+11+t2Dt(0.19)[L(x)[u(x,t)]]+L(x)[u(x,t)]+N(u(x,t),∂xu(x,t))+f(x,t), 0≤x,t≤1,(101)

where the spatial operator

L(x)[u(x,t)]=−∂xx[cos⁡(x)∂xxu(x,t)]+∂x[cosh⁡(x)∂xu(x,t)]−sinh⁡(x)u(x,t),(102)

and the nonlinear term is

N(u(x,t),∂xu(x,t))=sin⁡(u(x,t))+(∂xu(x,t))2.(103)

The boundary condition and initial condition conform to the exact u(x,t)=sin⁡(t)sin⁡(πx). The predicted results shown in Table 9 and Fig. 8 are obtained after 3th iterations of the quasilinearization procedure with randomly provided initial guess of solution and its derivatives. From Table 9, we can conclude that the proposed approach can produce accurate solution for the problems of high-order with a small number of basis functions both in terms of the spatial and temporal approximations. For the M=6, the method converges with N=14. And if we increase the M=12, the proposed method converges with N=20. And even with the small number of N=10 and M=6, the computational accuracy of the problem with the spatial differential operator of the fourth order is high around 10−5.

Figure 8: The maximal absolute Emax (left) and EH1 (right) errors as functions of the number of the RBFs

In Table 10, we display the effect of the parameter σ in the calculations. From this table, we can see that σ=0.3 provides a quite precise approximation in time. Therefore, in the following examples, we take σ=0.3 without further discussion.

4.2 (2+1)-Dimensional Problems

4.2.1 Example 5

Let us consider the linear fractional equation of the following form:

D(π/4)[u]+tD(0.3)[u]+t2D(0.2)[u]+t3D(0.1)[u]=(1+t2)Dt(0.33)[L(x,y)[u]]+(1+t3)Dt(0.22)[L(x,y)[u]]+(1+t4)Dt(0.11)[L(x,y)[u]]+(1+t)L(x,y)[u]+f(x,y,t),(104)

in the domain governed by

x=ρ(θ)cos⁡θ, y=ρ(θ)sin⁡θ, ρ(θ)=(cos⁡4θ+3.6+sin2⁡4θ)1/3, 0≤θ≤2π,(105)

as shown in Fig. 9.

Figure 9: Profile of the irregular domain

Here

L(x,y)[u]=div(D^(x,y)∇u),(106)

where the matrix of diffusivity is

D^=(cosh⁡(x+y),0.1(x2+y2)0.1(x2+y2),cosh⁡(x−y)).(107)

The boundary condition of the Robin type:

∂u∂n+(x2+y2)u=g(x,y,t), (x,y)∈∂Ω,(108)

and IC conform the exact solution uex(x,t)=exp⁡(x+y+t).

Table 11 shows the errors as the functions of N where M is fixed at 25, and the maximum error is about 10−7 for M=25 and N=10. From this table, it is clear that with the increasing of N, the errors decrease slowly indicating that N=10 is sufficient for the approximation. Additionally, for the small number of basis function and collocation nodes, the computational cost is low. For larger parameters, more computational resources should be deployed with less efficiency.

4.2.2 Example 6

In the last example, we consider a nonlinear problem

D(0.75)[u]+cosh⁡(t)D(0.3)[u]+t2D(0.2)[u]=(1+t2)Dt(0.33)[L(x,y)[u]]+L(x,y)[u]+u∂xu+u∂yu−u2+f(x,y,t),(109)

in the star-like domain shown in Fig. 10. Here

L(x,y)[u]=div(D^(x,y)∇u),(110)

Figure 10: Profile of irregular star-like solution domain

where the matrix of diffusivity is

D^=(cos⁡h(x+y),0.1(x+y)20.1(x+y)2,cosh⁡(x−y)).(111)

The Dirichlet boundary conditions and IC can be obtained from uex(x,t)=exp⁡(t)sinh⁡(x)sinh⁡(y). The results displayed in Table 12 are obtained after 3th iterations of the quasilinearization procedure where M=10 and t=1 and the initial guesses are randomly provided. From the table, we can see that the parameters M=10 and N=10 are sufficient to obtain approximate solutions and its derivatives where the three computed errors Emax, Erel, and EH1 are about 10−9. From this table, we notice that for the N=50 and M=10, the computational time is about 505 s for the fact that the full dense matrix system should be solved at each time step of the iterations. The computational cost can be reduced if the localized scheme is absorbed where the sparse matrix system can be solved efficiently. We refer the readers to [42]. The accuracy of the proposed method can also be verified from Figs. 11 to 13 where the domain absolute errors of u, ∂u/∂x, and ∂u/∂y are plotted. The maximum errors are occurred at the boundary.

Figure 11: Absolute errors of u

Figure 12: Absolute errors of ∂u/∂x

Figure 13: Absolute errors of ∂u/∂y

5  Conclusion

In order to solve (1+1)-dimensional and (2+1)-dimensional time-fractional nonlinear diffusion equations of multi-term accurately, we propose a simple technique in this paper. First, we convert the original problem into a new one under homogeneous initial conditions by subtracting the initial conditions from the unknown desired solutions. Then the key issue is to solve the problem with zero initial conditions. After that, we try to seek the function ug(x,t) analytical or numerically which satisfies the specify boundary conditions of the transformed problem. By the decomposition with the computed ug(x,t) which makes it possible that the original problem of homogeneous conditions can be solved efficiently. By this conversion, we can seek the function in the form of the truncated series over the modified radial basis functions. The modified radial basis functions and the linear combinations pose the properties satisfying the homogeneous boundary conditions. The collocation at the centers of the RBF transforms the problems into a system of ordinary equations which are solved by the recently proposed backward substitution technique. The numerical examples including linear and nonlinear problems prove the accuracy and efficiency of the proposed method. From the numerical results, we can see that with a small number of collocations nodes, the problems can be solved efficiently with high accuracy and converging order compared with the results in the literature. It should be noted here that this paper provides a general algorithm for solving time-dependent fractional equations which can be considered as an alternative tool for engineering.

Acknowledgement: The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.

Funding Statement: This research was funded by the National Key Research and Development Program of China (No. 2021YFB2600704), the National Natural Science Foundation of China (No. 52171272), and the Significant Science and Technology Project of the Ministry of Water Resources of China (No. SKS-2022112).

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: T. Hu, C. Huang, and S. Reutskiy; data collection: J. Lu; analysis and interpretation of results: J. Lin and S. Reutskiy; draft manuscript preparation: J. Lu and J. Lin. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: None.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Liu, L., Feng, L., Xu, Q., Zheng, L., Liu, F. (2020). Flow and heat transfer of generalized Maxwell fluid over a moving plate with distributed order time fractional constitutive models. International Communications in Heat and Mass Transfer, 116(6), 104679.

2. Zheng, L., Liu, Y., Zhang, X. (2012). Slip effects on MHD flow of a generalized Oldroyd-B fluid with fractional derivative. Nonlinear Analysis: Real World Applications, 13(2), 513–523.

3. Lyu, P., Vong, S. (2019). A high-order method with a temporal nonuniform mesh for a time-fractional benjamin-bona-mahony equation. Journal of Scientific Computing, 80(3), 1607–1628.

4. Bayarassou, K. (2021). Fourth-order accurate difference schemes for solving benjamin-bona-mahony-burgers (BBMB) equation. Engineering with Computers, 37(1), 123–138.

5. Dehghan, M., Shafieeabyaneh, N., Abbaszadeh, M. (2020). Application of spectral element method for solving sobolev equations with error estimation. Applied Numerical Mathematics, 158, 439–462.

6. Zhao, J., Fang, Z., Li, H., Liu, Y. (2020). A crank-nicolson finite volume element method for time fractional sobolev equations on triangular grids. Mathematics, 8(9), 1591.

7. Ren, J., Sun, Z. (2015). Efficient numerical solution of the multi-term time fractional diffusion-wave equation. East Asian Journal on Applied Mathematics, 5(1), 1–28.

8. Salehi, R. (2017). A meshless point collocation method for 2-D multi-term time fractional diffusion-wave equation. Numerical Algorithms, 74, 1145–1168.

9. Yang, X., Qiu, W., Zhang, H., Tang, L. (2021). An efficient alternating direction implicit finite difference scheme for the three-dimensional time-fractional telegraph equation. Computers & Mathematics with Applications, 102, 233–247.

10. Shen, S., Liu, F., Anh, V. (2019). The analytical solution and numerical solutions for a two-dimensional multi-term time fractional diffusion and diffusion-wave equation. Journal of Computational and Applied Mathematics, 345, 515–534.

11. Stynes, M., O’Riordan, E., Gracia, J. (2017). Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM Journal on Numerical Analysis, 55(2), 1057–1079.

12. Soori, Z., Aminataei, A. (2018). Sixth-order non-uniform combined compact difference scheme for multi-term time fractional diffusion-wave equation. Applied Numerical Mathematics, 131, 72–94.

13. Jin, B., Lazarov, R., Liu, Y., Zhou, Z. (2015). The Galerkin finite element method for a multi-term time-fractional diffusion equation. Journal of Computational Physics, 281, 825–843.

14. Huang, J., Zhang, J., Arshad, S., Tang, Y. (2021). A numerical method for two-dimensional multi-term time-space fractional nonlinear diffusion-wave equations. Applied Numerical Mathematics, 159, 159–173.

15. Hussain, M., Haq, S. (2019). Weighted meshless spectral method for the solutions of multi-term time fractional advection-diffusion problems arising in heat and mass transfer. International Journal of Heat and Mass Transfer, 129, 1305–1316.

16. Oruç, O., Esen, A., Bulut, F. (2019). A haar wavelet approximation for two-dimensional time fractional reaction-subdiffusion equation. Engineering with Computers, 35(1), 75–86.

17. Oruç, O. (2022). An accurate computational method for two-dimensional (2D) fractional Rayleigh-Stokes problem for a heated generalized second grade fluid via linear barycentric interpolation method. Computers & Mathematics with Applications, 118, 130–131.

18. Li, J., Dai, L., Nazeer, W. (2020). Numerical solution of multi-term time fractional wave diffusion equation using transform based local meshless method and quadrature. AIMS Mathematics, 5(6), 5813–5839.

19. Soltani Sarvestani, F., Heydari, M., Niknam, A., Avazzadeh, Z. (2019). A wavelet approach for the multi-term time fractional diffusion-wave equation. International Journal of Computer Mathematics, 96(3), 640–661.

20. Liu, Y., Yin, X., Liu, F., Xin, X., Shen, Y. et al. (2022). An alternating direction implicit legendre spectral method for simulating a 2D multi-term time-fractional Oldroyd-B fluid type diffusion equation. Computers & Mathematics with Applications, 113, 160–173.

21. Nikan, O., Tenreiro Machado, J., Golbabai, A. (2021). Numerical solution of time-fractional fourth-order reaction-diffusion model arising in composite environments. Applied Mathematical Modelling, 89, 819–836.

22. Zhou, J., Xu, D., Qiu, W., Qiao, L. (2022). An accurate, robust, and efficient weak Galerkin finite element scheme with graded meshes for the time-fractional quasi-linear diffusion equation. Computers & Mathematics with Applications, 124, 188–195.

23. Su, B., Jiang, Z. (2020). High-order compact finite volume scheme for the 2D multi-term time fractional sub-diffusion equation. Advances in Difference Equations, 2020, 1–22.

24. Ghafoor, A., Khan, N., Hussain, M., Ullah, R. (2022). A hybrid collocation method for the computational study of multi-term time fractional partial differential equations. Computers & Mathematics with Applications, 128, 130–144.

25. Piret, C., Hanert, E. (2013). A radial basis functions method for fractional diffusion equations. Journal of Computational Physics, 238, 71–81.

26. Hosseini, V., Chen, W., Avazzadeh, Z. (2014). Numerical solution of fractional telegraph equation by using radial basis functions. Engineering Analysis with Boundary Elements, 38, 31–39.

27. Liu, Q., Gu, Y., Zhuang, P., Nie, Y. (2011). An implicit RBF meshless approach for time fractional diffusion equations. Computational Mechanics, 48, 1–12.

28. Qiao, Y., Zhai, S., Feng, X. (2017). RBF-FD method for the high dimensional time fractional convection-diffusion equation. International Communications in Heat and Mass Transfer, 89, 230–240.

29. Reutskiy, S. (2017). A new semi-analytical collocation method for solving multi-term fractional partial differential equations with time variable coefficients. Applied Mathematical Modelling, 45, 238–254.

30. Reutskiy, S., Fu, Z. (2018). A semi-analytic method for fractional-order ordinary differential equations: Testing results. Fractional Calculus and Applied Analysis, 21(6), 1598–1618.

31. Reutskiy, S., Lin, J. (2018). A semi-analytic collocation method for space fractional parabolic PDE. International Journal of Computer Mathematics, 95(6–7), 1326–1339.

32. Lin, J., Zhang, Y., Reutskiy, S. (2021). A semi-analytical method for 1D, 2D and 3D time fractional second order dual-phase-lag model of the heat transfer. Alexandria Engineering Journal, 60(6), 5879–5896.

33. Tian, X., Reutskiy, S., Fu, Z. (2022). A novel meshless collocation solver for solving multi-term variable-order time fractional PDEs. Engineering with Computers, 38, 1527–1538.

34. Lin, J., Bai, J., Reutskiy, S., Lu, J. (2023). A novel RBF-based meshless method for solving time-fractional transport equations in 2D and 3D arbitrary domains. Engineering with Computers, 39, 1905–1922.

35. Pourbabaee, M., Saadatmandi, A. (2022). A new operational matrix based on Müntz-Legendre polynomials for solving distributed order fractional differential equations. Mathematics and Computers in Simulation, 194, 210–235.

36. Borwein, P., Erdelyi, T., Zhang, J. (1994). Müntz systems and orthogonal Müntz-Legendre polynomials. Transactions of the American Mathematical Society, 342(2), 523–542.

37. Bellman, R., Kalaba, R. (1965). Quasilinearization and nonlinear boundary-value problems. New York, USA: American Elsevier Publishing Company.

38. Liu, G., Gu, Y. (2001). A point interpolation method for two-dimensional solids. International Journal for Numerical Methods in Engineering, 50(4), 937–951.

39. Liu, G., Gu, Y. (2005). An introduction to meshfree methods and their programming. Berlin, Germany: Springer Science & Business Media.

40. Oruç, O. (2023). A strong-form meshfree computational method for plane elastostatic equations of anisotropic functionally graded materials via multiple-scale Pascal polynomials. Engineering Analysis with Boundary Elements, 146, 132–145.

41. Hadhoud, A., Abd Alaal, F., Abdelaziz, A., Radwan, T. (2021). Numerical treatment of the generalized time-fractional Huxley-Burgers’ equation and its stability examination. Demonstratio Mathematica, 54(1), 436–451.

42. Lin, J., Xu, Y., Zhang, Y. (2020). Simulation of linear and nonlinear advection-diffusion-reaction problems by a novel localized scheme. Applied Mathematics Letters, 99, 106005.

images
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.