Courtesy of  Santos Bravo Yuste  ☺   
Dpt. Fisica. UEx. Badajoz. Spain  
email:  santos@unex.es
http://www1.unex.es/eweb/fisteor/santos/sby.html

An explicit finite difference method for fractional diffusion equations

Here it is provided a Mathematica notebook that uses the explicit difference numerical algorithm discussed in 
Ref [1] = S. B. Yuste and L. Acedo, On an explicit finite difference method for fractional diffusion equations,  SIAM Journal of Numerical Analysis, 52, 1862-74 (2005)

The fractional diffusion equation (FDE) we are going to solve is (see Eq. (1.2) of Ref[1]) : 

			((∂u)/(∂t)) = K (((Subscript[, o])D)^(1-γ)) (((∂^2)u)/(∂(x^2)))	
			
In particular, we are going to consider two initial conditions:
	  □ (Case 1)   Dirac' delta initial condition,  u(x,0)= δ(x,0),  and the boundary condition u(-∞,t)=u(∞,t)=0 
	  	
	  □ (Case 2)   u(x,0)=x(1-x)  and Dirichlet boundary condition u(0,t)=u(1,t)=0

The solution u evaluated at position (Subscript[x, j])  at time (Subscript[t, m]),   u((Subscript[x, j]),(Subscript[t, m])), is denoted here as u[j,m,γ], where the third index γ  is the anomalous diffusion exponent of the FDE.

The explicit difference scheme/algorithm for solving the FDE (1) [see Eq. (3.3) of Ref[1]) :

BDF1 Coefficients:  w[k_, alfa_]≡(Subscript[ω, k]^((α)))  (see formula  (2.6) of Ref[1]) :

w[0,alfa_]:=1;w[k_,alfa_]:=w[k,alfa]=(1-(alfa+1)/k)*w[k-1,alfa]

In what follows, 
	lambda≡λ  (weight factor; see Eq. (17) of [2]) ,  
	 sfrac≡S  (parameter defined in Eq. (19b) of  [2]) and 
	 Nx determines the space discretization of the interval in which the equation is going to be solve (the total number of points (Subscript[x, j]) is   2Nx+1)

uME[j, m, gamma]  provides the value of the numeral solution at the point (Subscript[x, j])  at time  (Subscript[t, m])  when the  anomalous diffusion exponent γ is gamma.

uME[j_,m_,gamma_]:=uME[j,m,gamma]=Module[{kk},
kk=Sum[w[k,1-gamma]*{uME[j-1,m-1-k,gamma],uME[j,m-1-k,gamma],uME[j+1,m-1-k,gamma]},{k,0,m-1}];
uME[j,m-1,gamma]+sfrac*(kk[[1]]-2*kk[[2]]+kk[[3]])]

Case 1:   
Initial condition:  Dirac’ delta initial condition,  u(x,0)= δ(x,0)

An example of case 1.  
Parameters of the example:   gamma=1/2, ht=0.0005,Kg=1, sfrac=0.33 (⇒the algorithm is stable), deltax=0.260307 
(Notation:   gamma=γ,  ht=Δt  ,   Kg=K,   sfrac=S,  deltax=Δx )

The exact solution  (propagator)  for  γ=1/2  and  (Subscript[K, γ])=1:

uPropagador[x_,t_,1/2]:=(1/(4 (t^(3/4))))((1/(Gamma[(3/4)]))2 (Sqrt[t]) HypergeometricPFQ[{},{(1/2),(3/4)},-(((Abs[x])^4)/(256 t))]+Abs[x] (-(1/(Sqrt[π]))2 (t^(1/4)) HypergeometricPFQ[{},{(3/4),(5/4)},-(((Abs[x])^4)/(256 t))]+(1/(Gamma[(1/4)]))Abs[x] HypergeometricPFQ[{},{(5/4),(3/2)},-(((Abs[x])^4)/(256 t))]))

Here we clear the saved values of uME   (the evaluation of this cell is required if you are  going to work with several examples in a Mathematica session)

Clear[uME];
uME[j_,m_,gamma_]:=uME[j,m,gamma]=Module[{kk},
kk=Sum[w[k,1-gamma]*{uME[j-1,m-1-k,gamma],uME[j,m-1-k,gamma],uME[j+1,m-1-k,gamma]},{k,0,m-1}];
uME[j,m-1,gamma]+sfrac*(kk[[1]]-2*kk[[2]]+kk[[3]])]

{gamma=1/2,ht=0.0005,Kg=1.;sfrac=0.33,deltax=(Kg*ht^gamma/sfrac)^(1/2)}

2*Nx+1=number of spatial points

Nx=25;

Parameters of the example

{ sfracEfect=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2,Sx=2^(gamma-2)//N}

criterion=If [sfracEfect<=Sx,{"stable"," is not"},{"unstable"," is "}];Print["The algorithm is "<>criterion[[1]]<> " because sfracEfect="<>ToString[sfracEfect] <>criterion[[2]]<>" larger than (Subscript[S, \[Cross]])=2^(gamma-2)="<>ToString[Sx]<>" (see formula (3.12) of Ref[1])"  ]

Here we set (we approximate)  the initial condition  [u(x,0)= δ(x,0)]

uME[nx_,0,alfa_]:=1./deltax/;nx==0;

uME[nx_,0,alfa_]:=0/;nx≠ 0;

Table[{nx,uME[nx,0,1]},{nx,-Nx,Nx}]

Here we impose a Dirichlet boundary condition at the two ends of the interval (it is assumed that for the times considered the true solution satisfy, to a large extent, this boundary condition)

uL[t_]:=0.;
uR[t_]:=0;

uME[-Nx,m_,alfa_]:=uL[m*ht];
uME[Nx,m_,alfa_]:=uR[m*ht];

A small test to check that everything goes well:

{uME[1,9,0.5],uME[2,4,0.5]}

soluME[t, gamma]   is a function that generates a table  with the solution  u  evaluated at time t. The table has the format {x, u (x, t)}

soluME[t_,gamma_]:=Module[{mfin,kk,uMEfinalTab},
mfin=Floor[t/ht];kk=Timing[Do[uMEfinalTab=Table[{j*deltax,uME[j,m,gamma]},{j,-Nx,Nx}],{m,0,mfin}]][[1]];Print["mfin=",mfin,", t=",mfin*ht,",  ht=",ht,",  Delta x=",deltax, ",   S_γ=",sfrac,",   Computation time=",kk," :: Explicit method: Yuste-Acedo SIAM-JNA(2005)   "];uMEfinalTab]

Solution for t=0.005

tv=0.005;soluME[tv,gamma]

Solution for t=0.1

tv=0.1;soluME[tv,gamma]

Plot of the numerical solution for t=0.1

tv=0.1;figNum=ListPlot[soluME[tv,gamma],PlotRange->{{-Nx deltax,Nx deltax},All},PlotStyle->PointSize[0.02]]

Plot of the exact solution for t=0.1

tv=0.1;figExact=Plot[uPropagador[x,tv,1/2],{x,-6,6},PlotStyle->{Red,Thick},PlotRange->All]

Comparison of the exact and numerical solution for t=0.1

Show[figNum,figExact,PlotRange->All]

Case 2:   
Initial condition  u(x,0)=x(1-x). Boundary condition: u(0,t)=u(1,t)=0.

An example of case 2.  
Parameters of the example:  gamma=1/2, ht=0.0005, Kg=1, sfrac=0.33 (⇒the algorithm is stable), deltax=0.260307  
(Notation:   gamma=γ,  ht=Δt  ,   Kg=K,   sfrac=S,  deltax=Δx )

Exact solution   for  γ=1/2  and  (Subscript[K, γ])=1

uExactx1mx[x,t,γ,Kg,numterm] is the  exact solution  of the FDE (1) written as a (finite) series with numterm terms (the exact solution corresponds to numterm→∞)

uExactx1mx[x_,t_,γ_,Kg_,numterm_]:=8/Pi^3*Sum[1/(2*n+1)^3*Sin[(2*n+1)*Pi*x]*MittagLefflerE[γ,-Kg*(2*n+1)^2*Pi^2*t^γ ],{n,0,numterm}]

(The  Mittag-Leffler function for γ = 1/2  can be written in a simpler way:   MittagLefflerE[1/2,z]=(e^(z^2)) Erfc[-z])

Here we clear the saved values of uME   (the evaluation of this cell is required if you are  going to work with several examples in a Mathematica session)

Clear[uME];
uME[j_,m_,gamma_]:=uME[j,m,gamma]=Module[{kk},
kk=Sum[w[k,1-gamma]*{uME[j-1,m-1-k,gamma],uME[j,m-1-k,gamma],uME[j+1,m-1-k,gamma]},{k,0,m-1}];
uME[j,m-1,gamma]+sfrac*(kk[[1]]-2*kk[[2]]+kk[[3]])]

2*Nx+1=number of spatial points

Nx=10;

Parameters of the example

{gamma=0.5,deltax=1./(2*Nx),Kg=1.;sfrac=.33,ht=(sfrac*deltax^2/Kg)^(1/gamma)}

{ sfracEfect=sfrac*(Sin[(2*Nx-1)*Pi/(4*Nx)])^2,Sx=2^(gamma-2)//N}

criterion=If [sfracEfect<=Sx,{"stable"," is not"},{"unstable"," is "}];Print["The algorithm is "<>criterion[[1]]<> " because sfracEfect="<>ToString[sfracEfect] <>criterion[[2]]<>" larger than (Subscript[S, \[Cross]])=2^(gamma-2)="<>ToString[Sx]<>" (see formula (3.12) of Ref[1])"  ]

Initial condition : u(x,0)=x(1-x)

uME[nx_,0,alfa_]:=(nx+Nx)*deltax*(1-(nx+Nx)*deltax)

Table[{nx,uME[nx,0,1]},{nx,-Nx,Nx}]//N

Boundary condition: u(0,t)=u(1,t)=0

uL[t_]:=0.;
uR[t_]:=0;

uME[-Nx,m_,alfa_]:=uL[m*ht];
uME[Nx,m_,alfa_]:=uR[m*ht];

soluME[t, gamma]   is a function that generates a table  with the solution  u  evaluated at time t. The table has the format {x, u (x, t)}

soluME[t_,gamma_]:=Module[{mfin,kk,uMEfinalTab,eps=$MachineEpsilon},
mfin=Floor[(t+eps)/ht];kk=Timing[Do[uMEfinalTab=Table[{(j+Nx)*deltax//N,uME[j,m,gamma]},{j,-Nx,Nx}],{m,0,mfin}]][[1]];Print["mfin=",mfin," : t=",mfin*ht," : ht=",ht," : Delta x=",deltax, " : S_γ=",sfrac," : Computation time=",kk," : Explicit method: Yuste-Acedo SIAM-JNA(2005)   "];uMEfinalTab]

mfin=8;soluME[mfin*ht,gamma]

Plot of the numerical solution for  t=200*Δt=0.000136125

mfin=200;figNum=ListPlot[soluME[mfin*ht,gamma],PlotRange->All,PlotStyle->PointSize[0.02]]

Plot of the exact solution for t=200*Δt=0.000136125

tv=200*ht;figExact=Plot[uExactx1mx[x,tv,1/2,1,20],{x,0,1},PlotStyle->{Red,Thick},PlotRange->All]

Comparison of the exact and numerical solution for t=200*Δt=0.000136125

Show[figNum,figExact,PlotRange->All]