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

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)