# A New Variant of B-Spline for the Solution of Modified Fractional Anomalous Subdiffusion Equation

Authors: M. S. Hashmi; Zainab Shehzad; Asifa Ashraf; Zhiyue Zhang; Yu-Pei Lv; Abdul Ghaffar; Mustafa Inc; Ayman A. Aly
Journal: Journal of Function Spaces
 (2021)

---

## Abstract

The objective of this paper is to present an efficient numerical technique for solving time fractional modified anomalous subdiffusion equation. Anomalous diffusion equation has its role in various branches of biological sciences. B-spline is a piecewise function to draw curves and surfaces, which maintain its degree of smoothness at the connecting points. B-spline provides an active process of approximation to the limit curve. In current attempt, B-spline curve is used to approximate the solution curve of time fractional modified anomalous subdiffusion equation. The process is kept simple involving collocation procedure to the data points. The time fractional derivative is approximated with the discretized form of the Riemann-Liouville derivative. The process results in the form of system of algebraic equations, which is solved using a variant of Thomas algorithm. In order to ensure the convergence of the procedure, a valid method named Von Neumann stability analysis is attempted. The graphical and tabular display of results for the illustrated examples is presented, which stamped the efficiency of the proposed algorithm.


---

## Body

## 1. Introduction

The current study deals with the investigation of modified fractional anomalous subdiffusion equation using hybrid B-spline-based collocation method. The general form of the equation is(1)∂ψξ,τ∂τ=∑i=1nDτ1−αi∂2ψξ,τ∂ξ2+gξ,τa≤ξ≤b,0<τ≤τf,with BC(2)ψξ0,τ=p0τ,ψξN,τ=p1τ,and IC(3)ψξ,0=q0ξ,where Dτ1−αi is the Riemann-Liouville version of the fractional partial derivatives of order (1−αi), αi′s are the fractional parameters with 1≥αi>0 for i=1,2,⋯,n are fractional parameters appears in the equation. Here, g0τ, g1τ, h0ξ, and h1ξ are smooth functions, and ξ,τ are space and time variables.Fractional PDEs have implementation in numerous fields of engineering and science. Some of them are placed in the great book of Podlubny [1] and in some recent work can be found in [2–6]; their applications can be found in the field of electromagnetic and mechanical engineering.Liu et al. [7] proposed a finite difference technique for the modified anomalous subdiffusion equation. Liu et al. [8] presented two versions of finite element approximation with semi- and full discretizations in a finite domain and compact difference scheme by Wang et al. [9]. An unconditionally stable high-order scheme was proposed by Mohebbi et al. [10]. The proposed scheme has a nonlinear source term. In another paper by Mohebbi et al. [11], the solution of 2D modified anomalous subdiffusion equation was proposed, and this method was based on radial basis functions.Dag et al. [12] presented a numerical solution of Burger’s equation using B-spline collocation method. Zahra and Elkholy [13] use cubic splines to represent a numerical solution of fractional differential equations. Nonclassical diffusion problems by Ismail et al. [14], advection-diffusion problems by Nazir et al. [15], and fractional subdiffusion equation by Zhu et al. [16] are solved by using B-spline collocation methods. In recent times, Hashmi et al. [17, 18] has solved Hunter Saxton equation and space fractional PDE by cubic trigonometric and hybrid B-spline method. The numerical solution of time-space fractional PDEs using B-spline wavelet method is presented by Kargar and Saeedi [19].In present research, numerical formulation of the Riemann-Liouville derivative for anomalous subdiffusion equation described by Dehghan et al. [20] is used as time fractional derivative which is given as(4)Dτ1−αψξ,τ=∂∂τI0+αψξ,τ,where(5)I0+αψξ,τn+1−I0+αψξ,τndτ=dταΓ1+αψξ,τn+1+∑p=0n−1wp+1−wpψξ,τn−p+Rα,in which(6)Rα≤Cwpdτ+α,1=w0>w1>⋯>wp>0,and wp=p+1α−pα. Modified fractional anomalous subdiffusion equation have not been solved by hybrid B-spline collocation method yet. So, we solved modified fractional equation using hybrid B-spline collocation method. The proposed method is evaluated by Von Neumann stability analysis for convergence, which stamped that iterative method is stable without any condition over the domain. Numerical examples have been tested to highlight the performance of the numerical method. Numerical results are presented through tables and graphs. It shows that result very much agree with closed form of solution.The paper is organized as follows: the description of proposed method along with methodology is presented in Section2. In Section 3, initial state is calculated. Section 4 presents the convergence analysis of the iterative method. The accuracy of the iterative scheme is shown in Section 5, and Section 6 contains the conclusion.


## 2. Description of Proposed Method

Hybrid B-spline collocation method is utilized to resolve the modified anomalous subdiffusion equation. The approximate solutionψξ,τ to the analytical solution u¯ξ,τ is considered as(7)ψξ,τ=∑j=−1N+1ςjτH3,jξ.Here, time-dependent unknowns are denoted byςjτ, and H3,jξ is a hybrid cubic B-spline basis function of third order and given as(8)H3,jξ=γB3,jξ+1−γT3,jξ.The value ofγ plays a very significant role in the hybrid cubic basis function. The hybrid nature of the proposed method occurred in 0<γ<1, and at boundaries of this hybrid parameter, it produces trigonometric and polynomial B-spline. The nodal values are placed in Table 1.(9)Q1=γ6+1−γa1,Q2=2γ3+1−γa2,Q3=γ2h+1−γa3,Q4=γh2+1−γa5,Q5=−2γh2+1−γa6.Table 1
Values of hybrid functionH3,jξ.
ξξj−2ξj−1ξjξj+1ξj+2H3,j0Q1Q2Q10H3,j′0Q30−Q30H3,j′′0Q4Q5Q40Here,H3,j−1ξj,H3,jξj, and H3,j+1ξj are only three nonzero basis functions that are included over subinterval ξj,ξj+1 owing to the local support property of B-spline basis function. Hence, the approximate solution and its derivatives with respect to ξ at ξj,τn, by using (8), are(10)ψjn=Q1ςj−1n+Q2ςjn+Q1ςj+1n,ψxjn=Q3ςj−1n−Q3ςj+1n,ψxxjn=Q4ςj−1n+Q5ςjn+Q4ςj+1n.


### 2.1. Numerical Formulation

In order to apply the hybrid B-spline collocation method, we take (1) at ξj,tn and can be written in a form as(11)∂ψξ,τ∂τξ=ξj,τ=τn=Dτ1−α1ψξj,τn+Dτ1−α2ψξj,τn∂2ψξj,τn∂ξ2+gξj,τn.Now, using the above discretized form of forward difference scheme and Riemann-Liouville derivative as a time fractional derivative will have the form(12)ψjn+1−ψjndτ=Dτ1−α1ψξξjn+Dτ1−α2ψξξjn+gξj,τn.Simplifying the equation, we get(13)Q1ςj−1n+1+Q2ςjn+1+Q1ςj+1n+1+Q1ςj−1n+Q2ςjn+Q1ςj+1n=dτκ+λQ4ςj−1n+1+Q5ςjn+1+Q4ςj+1n+1+dτκ+λ∑p=0n−1dpα1+dpα2Q4ςj−1n−p+Q5ςjn−p+Q4ςj+1n−p,where wp+1−wp=dpαi,κ=dτα1/Γ1+α1andλ=dτα2/Γ1+α2.


## 2.1. Numerical Formulation

In order to apply the hybrid B-spline collocation method, we take (1) at ξj,tn and can be written in a form as(11)∂ψξ,τ∂τξ=ξj,τ=τn=Dτ1−α1ψξj,τn+Dτ1−α2ψξj,τn∂2ψξj,τn∂ξ2+gξj,τn.Now, using the above discretized form of forward difference scheme and Riemann-Liouville derivative as a time fractional derivative will have the form(12)ψjn+1−ψjndτ=Dτ1−α1ψξξjn+Dτ1−α2ψξξjn+gξj,τn.Simplifying the equation, we get(13)Q1ςj−1n+1+Q2ςjn+1+Q1ςj+1n+1+Q1ςj−1n+Q2ςjn+Q1ςj+1n=dτκ+λQ4ςj−1n+1+Q5ςjn+1+Q4ςj+1n+1+dτκ+λ∑p=0n−1dpα1+dpα2Q4ςj−1n−p+Q5ςjn−p+Q4ςj+1n−p,where wp+1−wp=dpαi,κ=dτα1/Γ1+α1andλ=dτα2/Γ1+α2.


## 3. Initial State

The numerical procedure of the iterative process can be initiated using initial vectorς0=ς−10,ς00,⋯⋯,ςN+10T. Initial conditions have been utilized to attain the initial vector by using the following procedure:(14)ψNξ,0=∑j=−1N+1ςj0B3,jξ,where ςj0′s are unknown parameters. We need the initial approximation ψNξ,0 to satisfy the following conditions:(15)ψξj,0=q0ξj=0,1,⋯,N,ψxξ0,0=q0′ξj=0,ψxξN,0=q0′ξj=N.This produce a square matrix of orderN+3 of the form(16)Qς0=B,where ς0=ς−10,ς00,⋯⋯,ςN+10T(17)Q=Q30−Q3⋯⋯⋯⋯0Q1Q2Q1⋱⋮0Q1Q2Q1⋱⋮⋱⋮⋮Q1Q2Q10⋯⋯⋯⋯Q30−Q3,and B=q0′ξ0,q0ξ0,⋯,q0ξN,q0′ξNT. The penta-diagonal system can be solved by the Thomas algorithm.


## 4. Stability Analysis

The stability of numerical technique can be examined by performing Von Neumann stability analysis. The stability of numerical schemes is closely related with numerical errors. The Von Neumann stability analysis provides us growth of the error in terms of Fourier series. For this, we used(18)ςln=ηneilθh,in equation, where θ is a rational number and i is a imaginary symbol. Utilizing (13), we get(19)Q1ςl−1n+1+Q2ςln+1+Q1ςl+1n+1+Q1ςl−1n+Q2ςln+Q1ςl+1n=dτκ+λQ4ςl−1n+1+Q5ςln+1+Q4ςl+1n+1+dτκ+λ∑p=0n−1wp+1−wpQ4ςl−1n−p+Q5ςln−p+Q4ςl+1n−p.Simplifying the equation in the presence of (18), we have(20)ηn+12Q1cosθh+Q2−dτκ+λ2Q4cosθh+Q5=ηn2Q1cosθh+Q2+dτκ+λ2Q4cosθh+Q5∑p=1n−1dpα1+∑p=1n−1dpα2ηn−p+ηn2Q4cosθh+Q5dτκ+λ.After simplification, we get(21)ηn+1=νηn+μ∑p=1n−1dpα1+∑p=1n−1dpα2ηn−p.Thus,ηn+1≤ηn.Hence, our described method is stable unconditionally.


## 5. Numerical Experiments

In current section, numerical examples are computed to prove feasibility and accuracy of the scheme, and the results are displayed through graphs and tables.Example 1.
Consider the following problem(22)∂ψξ,τ∂τ=Dτ1−α1+Dτ1−α2∂2ψξ,τ∂ξ2+gξ,τ,in the domain 0≤ξ≤1 and for the interval 0<τ≤1, with the boundary conditions(23)ψ0,τ=τ1+α1+α2,ψ1,τ=eτ1+α1+α2,and the IC(24)ψξ,0=0,where(25)gξ,τ=eξ1+α1τα1−Γ2+α1+α2Γ1+2α1+α2τ2α1+α2−Γ2+α2+α2Γ1+α1+2α2τα+2α2.
The exact solution is(26)u¯ξ,τ=eξτ1+α1+α2.Table2 describes the numerical errors at distinct rational values of γ, and it is evident that γ=0 is an appropriate solution than other case. So, trigonometric spline is the obvious choice to proceed further. In Table 3, the numerical display of the solutions at different values of τf is shown for α1=0.2, α2=0.8, and N=50. Table 4 shows the numerical display of solution for different values of α1 at final time τf=0.01 and α2=0.3. Table 5 shows the numerical display of solutions at different values of N when α1=0.6, α2=0.9, and τf=0.001. Figure 1(a) shows the graphical error when α1 varies at final time τf=0.05 but α2 remains same. Figure 1(b) shows the graphical representation of errors at different terminal time τf when α1 and α2 remains same. Figure 2(a) represents the comparison of numerical and approximate solution in the entire domain; on the other hand, the comparison at different values of τf is presented in Figure 2(b).Example 2.
Consider the problem(27)∂ψξ,τ∂τ=12Dτ1−α1+Dτ1−α2∂2ψξ,τ∂ξ2+gξ,τ,in the domain 0≤ξ≤1 and for the time 0<τ≤1
with the BC(28)ψ0,τ=τ1+α1+τ1+α2,ψ1,τ=eτ1+α1+τ1+α2,and the IC(29)ψξ,0=0,where(30)gξ,τ=eξ1+α1τα1−Γ2+α1Γ1+2α1τ2α1+1+α2τα2−Γ2+α2Γ1+2α2τ2α2.
The exact solution is(31)u¯ξ,τ=eξτ1+α1+τ1+α2.Table 2
The numerical display of error of Example1 when N=50, α1=0.3, α2=0.7, and τf=0.1 for different values of γ.
ξγ=0γ=0.5γ=10.106.71631×10−46.71643×10−46.71645×10−40.301.66881×10−31.66884×10−41.66885×10−40.502.12015×10−32.12019×10−42.12020×10−30.701.90480×10−31.90484×10−31.90484×10−30.908.75119×10−38.75133×10−38.75136×10−3Table 3
The numerical display of the solutions for Example1 at different values of τf with N=50 and α1=0.2, α2=0.8.
τf=0.1τf=0.01ξExactNumericalErrorExactNumericalError0.100.01105170.01172056.68809×10−40.00011050.00011726.6956×10−60.300.01349860.01516041.66181×10−40.00013490.00015161.6637×10−50.500.01648720.01859852.11125×10−30.00016490.00018602.1136×10−50.700.02013750.02203431.89679×10−30.00020140.00022041.8989×10−50.900.02459600.02546758.71428×10−40.00024600.00025478.7240×10−6Table 4
The numerical display of errors for Example1 at different values of α1 with α2=0.3, N=70, and τf=0.01.
α1=0.1α1=0.5ξExactNumericalErrorExactNumericalError0.100.00175160.00185811.06557×10−40.00027760.00029491.72452×10−50.300.00213940.00240412.64734×10−40.00033910.00038194.28495×10−50.500.00261300.00294943.36343×10−40.00041410.00046865.44384×10−50.700.00319160.00349383.02230×10−40.00050580.00055474.89091×10−50.900.00389820.00403711.38893×10−40.00061780.00064032.24703×10−5Table 5
Numerical display of solutions of Example1 of the TFPDE at different values of N with α1=0.6, α2=0.9, and τf=0.001.
N=70N=140ξExactNumericalErrorExactNumericalError0.103.495×10−83.706×10−82.11095×10−93.495×10−83.706×10−82.10877×10−90.304.269×10−84.793×10−85.24514×10−94.269×10−84.793×10−85.23974×10−90.505.214×10−85.880×10−86.66371×10−95.214×10−85.879×10−86.65684×10−90.706.368×10−86.967×10−85.98679×10−196.368×10−86.966×10−85.91222×10−90.907.778×10−88.053×10−82.75044×10−97.778×10−88.053×10−82.74761×10−9Figure 1
Numerical display of deviations from exact solution of Example1 for various values of α1 (a)and for distinct values of τf (b).
(a)(b)Figure 2
Graphical representation of exact and approximate solution for different values ofξ in entire domain (a) and various values of τf i.e., τf=0.1,0.2,0.3 in downward direction (b) of Example 1.
(a)(b)In Table6, the absolute error at different values of γ is displayed, and it shows that γ=0 gives better result. Table 7 shows the numerical solution for different values of α2 with final time τf=0.05 and α1=0.3. In Table 8, the numerical display of solutions for different values of τf when α1=0.5, α2=0.7, and N=70 is presented. Figure 3(a) represents absolute error Example 2 of various values of α2 with fixed τf. Figure 3(b) describes the graphical representation of errors at various values of terminal time τf when α1 and α2 remains same. Figure 4(a) represents the comparison of numerical and approximate solution in entire domain; on the other hand, the comparison at different values of α1 is presented in Figure 4(b).Table 6
The absolute error of Example2 when N=70, α1=0.5, α2=0.2, and τf=0.5 for different values of γ.
ξγ=0γ=0.3γ=10.101.79881×10−41.78163×10−41.77495×10−40.304.51450×10−44.47218×10−44.45572×10−40.505.72311×10−45.66922×10−45.64827×10−40.705.07258×10−45.02359×10−45.00454×10−40.902.27224×10−42.24926×10−42.24033×10−4Table 7
The numerical display of the solutions at different values ofα with N=40 and τf=0.05 for Example 2.
α1=0.3, α2=0.5α1=0.3, α2=0.9ξExactNumericalErrorExactNumericalError0.100.03485140.03538435.32893×10−40.02622320.02567545.47808×10−40.300.04256760.04386951.29287×10−30.03202910.03070431.32481×10−30.500.05199220.05364361.65143×10−30.03912040.03742721.69323×10−30.700.06350340.06503571.53223×10−30.04778180.04620441.57735×10−30.900.07756330.07830787.44550×10−40.05836080.05758907.71809×10−4Table 8
The numerical display of solutions of Example2 at different time level with N=70, α1=0.5, and α2=0.7.
τf=0.01τf=0.001ξExactNumericalErrorExactNumericalError0.100.00154510.00153876.4223×10−60.00004370.00004361.6572×10−70.300.00188720.00187281.4444×10−50.00005340.00005313.1578×10−70.500.00230510.00228651.8558×10−50.00006520.00006523.9949×10−70.700.00281540.00279681.8627×10−50.00007970.00007924.4739×10−70.900.00343880.00342841.0402×10−50.00009730.00009703.1454×10−7Figure 3
Graphical display of error for different values ofα2 (a) and τf (b) Example 2.
(a)(b)Figure 4
Graphical representation of exact and approximate solution for different values ofξ in entire domain (a) and various values of α1 i.e., α1=0.1,0.2,0.3 in upward direction (b) of example 2.
(a)(b)


## 6. Conclusion

A B-spline collocation method is utilized to find the numerical approximation to the solution curve of fractional form of anomalous subdiffusion equation in modified form. First, formulate the formula for the modified fractional anomalous subdiffusion equation using discretized form of Riemann-Liouville Derivative for the time fractional derivative. Then, we show the stability of B-spline method using the Fourier version of stability method, i.e., Von Neumann stability analysis. The iterative method is applied on numerical examples to show the applicability of described method. The graphs and tables demonstrate the numerical results. It shows that results are very much in agreement with exact solution. The technique can be further extended to higher dimension and nonlinear fractional PDEs in future.



---

*Source: 8047727-2021-07-05.xml*