# Nonpolynomial Spline Interpolation for Solving Fractional Subdiffusion Equations

Authors: Homan Emadifar; Faraidun K. Hamasalh
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

The nonpolynomial spline interpolation is proposed to distinguish numerical analysis from the senes boundary conditions, accurance error estimations. The idea used in this article is readily applicable to obtain numerical solution of nonpolynomial spline interpolation. These analyze the methods that are suitable for the numerical solution of subdiffusion equation. The method has been shown to be stable by using von Neumann technique. The accuracy and efficiency of the scheme are checked by several examples to obtain numerical tests.


---

## Body

## 1. Introduction

In recent years, the derivation of solutions to fractional differential equations has become a hot topic in many areas of applied science and engineering. There are still not enough extremely precise numerical methods, although there have been many works on numerical methods for fractional differential equations. It is also very important to conduct research on fractional-order differential equations in unconstrained fields. A popular approach to solving these problems is to use artificial boundary conditions (ABCs) [1]. There are several applications of spline methods to the numerical solution of differential equations and, in particular, fractional differential equations [2–6]. The derivation of solutions to fractional differential equations has become a first-rate topic in many areas of applied science and engineering. Since fractional derivatives provide more accurate models than integer derivatives for most real-world problems, fractional differential equations are used to formulate a considerable number of applied problems. Applications are diverse and include viscoelastic and viscoelastic flows, control theory, transport problems, tumor evolution, random walks, continuum mechanics, and turbulence [7–13]. Finite differences, Fourier method, orthogonal spline collocation, implicit schemes with alternating direction, and compact finite difference methods have been successfully used in [14–22] to solve fractional subdiffusion equations. The author of [23] has presented an implicit numerical method for the fractional diffusion equation, where the fractional derivative is discretized by spline and the Crank–Nicolson discretization is used for the time variable. The authors of [24] have developed a novel nonpolynomial spline method for solving second-order hyperbolic equations, and their results are numerically more accurate than some finite difference methods, see [25–28]. The authors of [29] presented a quadratic nonpolynomial spline approach to approximate the solution of a system of second-order boundary value problems associated with a one-sided obstacle, and compared to collocation, finite difference, and certain common polynomial spline approaches, we used contact problems and produced approximations that were more accurate.We consider the fractional differential equation on a bounded domain expressed by [1](1)Dtαux,t−uxxx,t=fx,t,x,t∈Ω,(2)ux,0=ψx,x∈0,1,(3)u0,t=ϕt,t>0,(4)u1,t=φt,t>0,where Ω=x,t|0≤x≤1,0≤t<+∞ and 0<α<1.


## 2. Basic Definition

Definition 1.
Letu(x) be a function defined on (a, b), then the Riemann–Liouville fractional derivative is of the following form [12]:(5)Dαux=1Γm−αdmdxm∫0xx−tm−α−1utdt,α>0,m−1<α<m.Definition 2.
Letu(x) be a function defined on (a, b), then the Caputo fractional derivative is of the following form [12]:(6)Dαux=1Γm−α∫0xx−tm−α−1umtdt,α>0,m−1<α<m.


## 3. Nonpolynomial Spline Method

In this section, we present nonpolynomial spline method for solving the problem equations (1)–(4) in the finite interval 0,T. For the positive integers n and k, we take h=1/n, τ=T/k as the spatial stepsize and temporal stepsize, respectively. We use the notations xi=ih, 0≤i≤n, tj=jτ, 0≤j≤k, and uij=uxi,tj. Let us consider a mesh with nodal points xi=ih on 0,1 such that(7)Δ:0=x0<x1<x2<…<xn−1<xn=1.For each segmentxi,xi+1,i=0,1,2,…,n−1, the nonpolynomial spline Qix,tj has the following form:(8)Qix,tj=aitjcoskx−xi+bitjsinkx−xi+citjx−xi5+ditjx−xi4+eitjx−xi3+κitjx−xi2+gitjx−xi+ritj,i=0,…,n,where k is the frequency of the trigonometric part of the spline functions which can be real or pure imaginary and which will be used to raise the accuracy of the method [30]. Let ux,t be the exact solution, and let Sij be an approximation to uij obtained by the segment Qix,tj passing through the points xi,Sij and xi+1,Si+1j.To derive the coefficients of the equation (8), we first define(9)Qixi,tj=Sij,Qi′xi,tj=mij,Qi2xi,tj=Mij,Qi6xi,tj=LijQixi+1,tj=Si+1j,Qi′xi+1,tj=mi+1j,Qi2xi+1,tj=Mi+1j,Qi6xi+1,tj=Li+1j.By algebraic manipulation, we getaitj,bitj,citj,ditj,eitj,κitj,gitj, and ritj. Using the continuity of the third, fourth, and fifth derivatives, that is, Qi−1ρxi,tj=Qiρxi,tj,ρ=3,4, and 5, we obtain the following useful consistency relation in terms of the second derivative of spline Mi and ui,(10)α1Mi−3j+Mi+3j+α2Mi−2j+Mi+2j+α3Mi−1j+Mi+1j+α4Mij=1h2ui+3j+ui−3j+β1ui+2j+ui−2j+β2ui+1j+ui−1j+β3uij,i=3,…,n−3.By expanding equation (10) in Taylor series about xi at each fixed t=tj, 0≤j≤M, we obtain the following local truncation error:(11)tij=−2+2β1+2β2+β3uij+h2d2udx2ij−9+2α1+2α2+2α3+α4−4β1−β2+1/12h4d4udx4ij−81+108α1+48α2+12α3−16β1−β2+1360h6d6udx6ij−729+2430α1+480α2+30α3−64β1−β2+120160h8d8udx8ij−6561+40824α1+3584α2+56α3−256β1−β211814400h10d10udx10ij−59049+590490α1+23040α2+90α3−1024β1−β2+1239500800h12d12udx12ij−531441+7794468α1+135168α2−132α3−4096β1+143589145600h14d14udx14ij−4782969+96722262α1+745472α2−182α3−16384β1−β2+…,where α1,α2,α3,α4,β1,β2, and β3 are given in [1]. In this paper, we use fourth- and sixth-order method. From [30], α1,α2,α3,α4,β1,β2, and β3 are α1=50,α2=10,α3=−450,α4=900,β1=24,β2=15, and β3=−80 for fourth-order method and α1=1/42,α2=20/7,α3=397/14,α4=1208/21,β1=24,β2=15, and β3=−80 for sixth-order method.Lemma 1.
LetQ is an estimation of the smoothness function at class C2, then the error estimation obtained as(12)ui″−Qi″≤h212Qi4+Oh4,Qi″−ui″≤h212ui4+Oh4,ui′−Qi′≤h4180Qi5+Oh5,Qi′−ui′≤h4180ui5+Oh5.Proof.
The above relations are due to the expansion of Taylor sentencesMi+1,Mi−1,mi+1,mi−1,ui+1, ui−1, at point xi can be obtained as follows:(13)ui″=Qi″+h212Qi4+h4240Qi6−h66048Qi8+Oh7,i=1,…,n−1,I,Qi″=ui″−h212ui4+h4360ui6+17h660480ui8+Oh7,i=1,…,n−1.II..
Also, we have the same way.(14)ui′=Qi′+h4180Qi5+h61512Qi7−h814400Qi9+Oh9,i=1,…,n−1,III,Qi′=ui′−h4180ui5+h61512ui7+h825920ui8+h9,i=1,…,n−1.IV..
And this is the end of the proof.The local truncation error for the first- and second-order spline method is obtained as follows:(15)IISecond−ordermethod:ti=307h6ui6+Oh7,IIISixth−ordermethod:ti=h8ui8252+Oh9,IVEighth−ordermethod:ti=65629h10ui1027936090+Oh11.To obtain unique solution for the nonlinear system equation (10), we need four more equations. We define the following identities:(16)∑k=04γkukj+h2∑k=112ηkuk′′j+ν1h14u140j=0,ati=1,∑k=05μkukj+h2∑k=112σkuk′′j+ν2h14u140j=0,ati=2,∑k=05μkun−kj+h2∑k=112σkun−k′′j+νn−2h14u14nj=0,ati=n−2,∑k=04γkun−kj+h2∑k=112ηkun−k′′j+νn−1h14u14nj=0,ati=n−1.by using Taylor’s expansion, we obtain the unknown coefficients in equation (16) as follows:(17)γ0,γ1,γ2,γ3,γ4=65,−104,14,24,1,ν1=νn−1=1210210269217/326918592000,η1,η2,η3,η4,η5,η6,η7,η8,η9,η10,η11,η12=−2248215317/19958400,4539175/17600,−6055918291/6652800,9918918899/4989600,892256279/285120,18019157507/4989600,−30650022317/9979200,2380569353/1247400,−16851712321/19958400,503717713/1995840,−130447723/2851200,18976637/4989600μ0,μ1,μ2,μ3,μ4,μ5=26,14,−80,15,24,1,ν2=νn−2=521085679991373621248000,σ1,σ2,σ3,σ4,σ5,σ6,σ7,σ8,σ9,σ10,σ11,σ12,=−7712745923159667200,2598859335702400,−1922274760153222400,835179191911404800,−37121331073193344,1588304771280,−26143703812280960,2846619240739916800,−720403331322809600,16761737177408,−434354892534400,1137958737933600,where from equation (1), we can write Mij in the form(18)Mij=∂2uij∂x2=0cDtαuxi,tj−fxi,tj.Now, we need to use the definition of discrete approximation of fractional derivative in time of any appropriate functionz2, 10 as(19)Dtαzxi,tj=a∑q=0j−1Rj,qziq+1−ziq+Oτ2−α,where a=1/Γ2−ατα and Rj,q=j−q1−α−j−q−11−α.Then, from equations (18) and (19), we can express Mij in the following form:(20)Mij=a∑q=0j−1Rj,quiq+1−uiq−fij.By substituting equations (20) in (10), we have following systems:(21)α1a−1h2ui−31+ui+31+α2a−β1h2ui−21+ui+21+α3a−β2h2ui−11+ui+11+α4a−β3h2ui1=aα1ui−30+ui+30+α2ui−20+ui+20+α3ui−10+ui+10+α4ui0+α1fi−31+fi+31+α2fi−21+fi+21+α3fi−11+fi+11+α4fi1,fori=3,4,…,n−3, if j=1, and(22)α1a−1h2ui−3j+ui+3j+α2a−β1h2ui−2j+ui+2j+α3a−β2h2ui−1j+ui+1j+α4a−β3h2uij=aα1ui−3j−1+ui+3j−1+α2ui−2j−1+ui+2j−1+α3ui−1j−1+ui+1j−1+α4uij−1+α1fi−3j+fi+3j+α2fi−2j+fi+2j+α3fi−1j+fi+1j+α4fij+μij,for i=3,4,…,n−3, if j=2,3,…,k, with(23)μij=a∑q=0j−2Rj,qα1ui−3q+ui+3q+α2ui−2q+ui+2q+α3ui−1q+ui+1q+α4uiq−a∑q=0j−2Rj,qα1ui−3q+1+ui+3q+1+α2ui−2q+1+ui+2q+1+α3ui−1q+1+ui+1q+1+α4uiq+1.Substituting equations (20) in (16), we get(24)∑k=04γkukj+h2∑k=112ηka∑q=0j−1Rj,qukq+1−ukq−fkj+ν1h14u140j=0,ati=1,∑k=05μkukj+h2∑k=112σka∑q=0j−1Rj,qukq+1−ukq−fkjν2h14u140j=0,ati=2,∑k=05μkun−kj+h2∑k=112σka∑q=0j−1Rj,qun−kq+1−un−kq−fn−kj+νn−2h14u14nj=0,ati=n−2,∑k=04γkun−kj+h2∑k=112ηka∑q=0j−1Rj,qun−kq+1−un−kq−fn−kj+νn−1h14u14nj=0,ati=n−1.In the above systems, for eachj=1,2,…,k, clearly, we have n linear equations in terms of n unknowns u1j,u2j,…,unj.


## 4. Stability

In this section, we analyze the stability of equation (21), by means of the von Neumann method. We consider(25)uij=ξjexpIiΦh,where I=−1 and Φ is a real mode number and ξ is a complex number. We substitute equations (25) into (21) to obtain(26)ξ1α1a−1h2expIi−3θ+expIi+3θ+α2a−β1h2expIi−2θ+expIi+2θ+α3a−β2h2expIi−1θ+expIi+1θ+α4a−β3h2expIiθ=aξ0α1expIi−3θ+expIi+3θ+α2expIi−2θ+expIi+2θ+α3expIi−1θ+expIi+1θ+α4expIiθ,where θ=Φh, and this equation reduces to ξ1=ρξ0 in which ρ=x/y with(27)x=a2α1cos3θ+2α2cos2θ+2α3cosθ+α4,y=a2α1cos3θ+2α2cos2θ+2α3cosθ+α4−1h22cos3θ+2β1cos2θ+2β2cosθ+β3.Then, we obtain(28)ρ=xx−1/h22cos3θ+2β1cos2θ+2β2cosθ+β3,in whichx≥0 and(29)2cos3θ+2β1cos2θ+2β2cosθ+β3≤0,hold for all realθ and α1,α2,α3,α4,β1,β2, and β3 in fourth- and sixth-order method. Hence, the condition ρ≤1 is satisfied which implies ξ1≤ξ0. Similarly, by substituting equations (25) into (22), we obtain(30)ξjα1a−1h2expIi−3θ+expIi+3θ+α2a−β1h2expIi−2θ+expIi+2θ+α3a−β2h2expIi−1θ+expIi+1θ+α4a−β3h2expIiθ=aξj−1α1expIi−3θ+expIi+3θ+α2expIi−2θ+expIi+2θ+α3expIi−1θ+expIi+1θ+α4expIiθ+μij,in which(31)μij=a2α1cos3θ+2α2cos2θ+2α3cosθ+α4expIθi∑q=0j−2Rj,qξq−ξq+1.From the denoted values ofx and y, the equation (30) can be simplified as(32)yξj=xξj−1+x∑q=0j−2Rj,qξq−ξq+1.Then, we obtain(33)ξj=ρξj−1+ρ∑q=0j−2Rj,qξq−ξq+1,j=2,3,…,k.Now, from this equation, it is obtained thatξj≥ξ0 for all j=2,3,…,k (see [3]). Hence, the nonpolynomial spline method is unconditionally stable.


## 5. Numerical Results

In this section, the wealth of the presented method nonpolynomial spline method will be tested on three different examples, which are widely used in the literature. All problems were solved with best error estimations, and their results were compared with analytical or spline solutions given in the literature. In the following section, detailed information about this comparison will be presented.Example 1.
Consider problem equations (1)–(4) with(34)fx,t=x−12exp−xt3x−12Γ4+α/6−21−10x+x2tα−4x2x−1214x2−14x+3,x,t∈Ωψx=x4x−14,x∈0,1,ϕt=t3+α,φt=0,t∈0,T.
The exact solution of the problem is as follows [7]:(35)ux,t=x−14exp−xt3+α+x4,x,t∈Ω.
This problem is solved using the methods described in Section2 for T=1,2,3. The results are reported in Tables 1 and 2.Table 1
Maximum absolute error forh = 1/20, τ=1/2000.
TFourth-order methodSixth-order method17.7006e−43.1388e−728e−39.1704e−733.38e−21.7069e−6Table 2
2-norm error forh = 1/20, τ=1/2000.
TFourth-order methodSixth-order method12.5e−39.4389e−722.59e−22.7779e−631.086e−15.1856e−6Example 2.
Consider problem equations (1)–(4) with(36)fx,t=Γ3Γ3−αt2−αsin2πx+2π2t2sin2πx,x,t∈Ωψx=0,x∈0,1,ϕt=0,φt=0,t∈0,T.
It can be checked that the exact solution isux,t=t2sin2πx [10].
This problem is solved using the methods described in Section2 for T=1/4,1/2,1. The results are reported in Table 3.Table 3
2-norm error forh = 1/16, τ=1/4096.
TFourth-order methodSixth-order methodThe method in [10]¼8.3e−39.8444e−81.0666e−6½3.39e−22.5638e−83.8858e−611.3e−13.0239e−71.4463e−5


## 6. Conclusion

In this paper, we constructed the application of the nonpolynomial spline method to obtain the numerical point solution of fractional partial differential equations. Numerical results showed that these techniques achieve accuracy, in numerical results illustrative graphs, the Figure1 the comparison of the present results with the result of [10] shows the efficiency of the suggested technique. Also, the stability and error estimates of the methods are investigated and shown some figures of the results. Moreover, illustrations for the stability region of the schemes are derived.Figure 1
Comparison of the solution u(x, t) with respect values of obtained by scheme spline of Example 5.2 in [10].



---

*Source: 7354121-2022-09-05.xml*