Solving a Caputo fractional diffusion equation in cylindrical coordinates

Question

I would like to solve this following equation with d, the order of the factional derivative (in the sense of Caputo) :

So I tried the following code :

R1 = 0.001; 
R2 = 1; 
K = 1;
d = 0.9;

eq = 1/K  CaputoD[T[t, r, θ], {t, d}] - 
   1/ (r^2)*Derivative[0, 0, 2][T][t, r, θ] - 
   1/r*Derivative[0, 1, 0][T][t, r, θ] - 
   Derivative[0, 2, 0][T][t, r, θ];

(*Initial and boundary conditions*)
ic = T[0, r, θ ] == 0;
ic1 = Derivative[1, 0, 0][T][0, r, θ] == 0;

bc = DirichletCondition[T[t, r, θ] == 1, 
   0 <= θ <= 2 Pi && r==R1];


pbc1 = PeriodicBoundaryCondition[
   T[t, r, θ], θ == 2 Pi + 0.01 && R1 < r < R2, 
   TranslationTransform[{0, -2 Pi}]];
pbc2 = PeriodicBoundaryCondition[
   T[t, r, θ], θ == -0.01 && R1 < r < R2, 
   TranslationTransform[{0, 2 Pi}]];

(*solution w/ NDSolveValue*)

sol = NDSolveValue[{eq == 0, ic, bc, pbc1, pbc2}, 
  T, {r, R1, R2}, {θ, -0.01, 2 *Pi + 0.01}, {t, 0, 10}]

It naturally works for d = 1 but fails for other fractional values. A quick search for previous questions asked such as this one here led me to realize that NDSolve cannot handle Caputo yet and that you have to rely on DSolve with some Laplace transforms to get solutions to such equations. A numerical solution would alternatively require a manual discretization using tools like pdetoode.

However, from my humble understanding of the answers presented in this post, the numerical resolution with pdetoode relies on a n-th order polynomial approximation and I have no idea how to adapt this to my 2D cylindrical coordinates problem. Other posts showcase different methods like the Haar wavelet method here but again, I can't see how to adapt this one to solving my problem.

Surprisingly, when i try to compute this equation instead :

using this modified code :

R1 = 0.001; 
R2 = 1; 
K = 1;
d = 0.9;

eq = 1/K  (Derivative[2, 0, 0][T][t, r, θ] + 
      CaputoD[T[t, r, θ], {t, d}]) - 
   1/ (r^2)*Derivative[0, 0, 2][T][t, r, θ] - 
   1/r*Derivative[0, 1, 0][T][t, r, θ] - 
   Derivative[0, 2, 0][T][t, r, θ];

(*Initial and boundary conditions*)
ic = T[0, r, θ ] == 0;
ic1 = Derivative[1, 0, 0][T][0, r, θ] == 0;

bc = DirichletCondition[T[t, r, θ] == 1, 
   0 <= θ <= 2 Pi && R1 <= r <= R1 + 0.001];


pbc1 = PeriodicBoundaryCondition[
   T[t, r, θ], θ == 2 Pi + 0.01 && R1 < r < R2, 
   TranslationTransform[{0, -2 Pi}]];
pbc2 = PeriodicBoundaryCondition[
   T[t, r, θ], θ == -0.01 && R1 < r < R2, 
   TranslationTransform[{0, 2 Pi}]];

(*solution w/ NDSolveValue*)

sol = NDSolveValue[{eq == 0, ic, ic1, bc, pbc1, pbc2}, 
  T, {r, R1, R2}, {θ, -0.01, 2 *Pi + 0.01}, {t, 0, 10, 1}]

I do get an interpolating function, which I can then visualise with

frames = 
  Table[DensityPlot[
    sol1[t, Sqrt[x^2 + y^2], 
     Mod[ArcTan[x, y], 2 Pi]], {x, y} ∈ reg2D, 
    Axes -> {True, True, False}, AxesLabel -> {x, y}, 
    PlotRange -> All, AxesStyle -> Directive[Purple, 12], 
    ColorFunction -> "Rainbow", ColorFunctionScaling  -> False, 
    PlotLegends -> BarLegend[{"Rainbow", {0, 1}}], 
    PlotLabel -> Row[{"t=", t}]], {t, 0, 10, 1}];

Export["test.mp4", frames]

So here are my two questions :

1 - How can I solve numerically the first equation (without the 2nd order derivative wrt time)

2 - I have read that NDSolve cannot handle Caputo but if I add a second order derivative w.r.t. time, it solves the equation. I'd like to know if the solution I get truly is the solution to my second equation or if NDSolve somehow disregards Caputo to compute the whole thing. The reason why I ask this is I have tested several values of fractional order d and there doesn't seem to be much change, hence the question.

I apologize for the long post. I am more interested in solving the first equation but if you have any idea about my second question, I would also appreciate your insight on the matter.

Have a great day.

Answer 1

First of all, let me extend the useful part of my comments: as to 2nd question, your guess is right. NDSolve doesn't parse the CaputoD correctly in 2nd sample. A simpler example showing the issue:

{NDSolveValue[{y'[x] + aaadvddd CaputoD[y[x], {x, 0.5}] == Sin[x], y[0] == 1}, 
              y, {x, 0, 2 Pi}], 
 NDSolveValue[{y'[x] == Sin[x], y[0] == 1}, y, {x, 0, 2 Pi}]} // 
 ListPlot[#, PlotMarkers -> {o, x}] &

As we can see, the meaningless term aaadvddd CaputoD[y[x], {x, 0.5}] is ignored silently.

Then let's deal with 1st question. The solution in this post can be extended to 2D in a staightforward manner. As usual, I'll use pdetoode for spatial discretization:

tend = 1/40;
R1 = 1/1000;
R2 = 1;
d = 9/10;
With[{T = T[t, r, θ]},
 eq = CaputoD[T, {t, d}] - D[T, {θ, 2}]/r^2 - D[T, r]/r - D[T, {r, 2}] == 0;
 ic = T == 0 /. t -> 0;
 bc = {T == 1 /. r -> R1, D[T, r] == 0 /. r -> R2};]
 
{points@t, points@r, points@θ} = {15, 25, 25};    
{domain@t, domain@r, domain@θ} = {{0, tend}, {R1, R2}, {-Pi, Pi}};

CGLGrid[xl_, xr_, n_Integer /; n > 1] := 
 1/2 (xl + xr + (xl - xr) Cos[(π Range[0, n - 1])/(n - 1)])
grid@t = CGLGrid[##, points@t] & @@ domain@t;
(grid@# = Array[# &, points@#, domain@#]) & /@ {r, θ};

difforder = 2;
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = pdetoode[T[t, r, θ], t, grid /@ {r, θ}, difforder, {False, True}];

del = #[[2 ;; -2]] &;
SetAttributes[ca, Listable]    
ode = del@ptoofunc[eq /. CaputoD[a_, b_] :> ca[a]]; // AbsoluteTiming
(* {0.357973, Null} *)
odeic = del@ptoofunc@ic; // AbsoluteTiming
(* {0.0016925, Null} *)    
odebc = ptoofunc@bc; // AbsoluteTiming
(* {0.0028338, Null} *)

Remark

1. ca with Listable attribute is introduced to temporarily replace CaputoD because the new-introduced CaputoD is incredibly slow on expression that actually doesn't need to be calculated. Minimal example: Clear[f, t]; Table[CaputoD[f[i, j][t], {t, 1/2}], {i, 5}, {j, 5}]; // AbsoluteTiming.

2. The periodic b.c. in $\theta$ direction is set by 5th argument of pdetoode. Here {False, True} means 1st dimension i.e. $r$ direction is not periodic but 2nd dimension i.e. $\theta$ direction is periodic.

Then we use n-th order polynomial whose coefficient is determined by function values at Chebyshev-Gauss-Lobatto grid to approximate solution in $t$ direction. An InterpolatingFunction is generated as the output for convenience:

range = Range[0, points@t - 1];
coef[x_, y_] = c[x, y] /@ range;
poly[x_, y_][t_] = coef[x, y] . t^range;

{approxode, approxic, approxbc} = {ode, odeic, odebc} /. T -> poly /. 
    ca[a_] :> CaputoD[a, {t, d}]; // AbsoluteTiming
(* {2.11858, Null} *)

var = Outer[coef, grid@r, grid@θ] // Flatten;
sys = Flatten@{Table[approxode, {t, Rest@#}], Table[approxbc, {t, #}], approxic} &@
    grid@t; // AbsoluteTiming
{barray, marray} = CoefficientArrays[sys, var]; // AbsoluteTiming
(* {6.03323, Null} *)
(* {6.30408, Null} *)

csol = LinearSolve[marray, -N[barray]]; // AbsoluteTiming
(* {1.87547, Null} *)

nsol = Block[{c}, Evaluate@var = csol;
    ListInterpolation[
     Table[poly[x, y][t], {t, grid@t}, {x, grid@r}, {y, grid@θ}], 
     grid /@ {t, r, θ}]]; // AbsoluteTiming
(* {0.633282, Null} *)

Finally, visualization. I use ToElementMesh to discretize the region first to speed up Plot3D:

<< NDSolve`FEM`

reg = 
  ToElementMesh[ImplicitRegion[R1^2 < x^2 + y^2 < R2^2, {x, y}], 
   MaxCellMeasure -> Pi R2^2/(4 10^3)];

lst = 
   Table[Plot3D[nsol[t, Sqrt[x^2 + y^2], ArcTan[x, y]], {x, y} ∈ reg, 
     Exclusions -> None, PlotRange -> {-1, 1}], {t, 0, tend, 
     tend/points@t}]; // AbsoluteTiming
(* {8.53243, Null} *)

lst // ListAnimate

Oops, I forgot RevolutionPlot3D is a better choice in this case:

lst = 
   Table[RevolutionPlot3D[nsol[t, r, θ], {r, R1, R2}, {θ, -Pi, Pi}, 
     PlotRange -> {-1, 1}], {t, 0, tend, tend/points@t}]; // AbsoluteTiming
(* {1.87516, Null} *)
    
lst // ListAnimate

As we can see, the solution reaches equilibrium quite fast, that's the reason I modify tend to 1/40.