# Solving a Caputo fractional diffusion equation in cylindrical coordinates

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.

• 1. What b.c. are you trying to set at r==R1? (The definition of bc looks a bit suspicious. ) 2. You're aware that zero Nemann condition will be set at r==R2, right? 3. In this toy sample the solution is actually axisymmetric, is the $θ$ dependence necessary for you in real case? 4. You wrote {t, 0, 10, 1}, do you really know what it means?: mathematica.stackexchange.com/a/240406/1871 Dec 5, 2022 at 2:32
• 5. As to 2nd question, your guess is right, NDSolve doesn't parse the CaputoD correctly, just add a meaningless coefficient to the CaputoD term and try, or use NDSolveFEMGetInactivePDE. A much simpler example showing the issue: NDSolveValue[{y'[x] + aaadvddd CaputoD[y[x], {x, 0.5}] == Sin[x], y == 1}, y, {x, 0, 30}] Dec 5, 2022 at 2:40
• @xzczd 1- I guess I was going for a "hot inner wall" type of bc, if we are reasoning with heat transfer. I guess a bc = DirichletCondition[T[t, r, [Theta]] == 1, 0 <= [Theta] <= 2 Pi && r == R1]; would have been enough. Dec 5, 2022 at 13:29
• @xzczd 2- I was not aware of this. But coincidentally, I was working previously with an additional boundary condition cb = NeumannValue[0, r == R1] + NeumannValue[0, r == R2] which I ended up removing since whether it was there or not it was there didn't seem to change the solution. I guess I now know why. 3 - I do have a Theta dependence in my real case. 4- I wasn't aware of this either and made this mistake probably because of my habits in python. I'll correct this right now. 5 - Thanks for confirming it. I now know that results for my second equation are not reliable. Dec 5, 2022 at 13:40

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 == 1},
y, {x, 0, 2 Pi}],
NDSolveValue[{y'[x] == Sin[x], y == 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,
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:

<< NDSolveFEM

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.

• I thank you for this answer. By the way, when I used your code, I ended up with 3 error messages along the line of  InterpolatingFunction::dmval: Input value {0,1.41401,-2.35619} lies outside the range of data in the interpolating function. Extrapolation will be used. as well as General::stop: Further output of InterpolatingFunction::dmval will be suppressed during this calculation. I imagine the 1st error has to do with points in the plotrange not being in the interpolation domain, but what about the second error. Do you perhaps know the reason for it ? Dec 6, 2022 at 14:47
• @ConfuzzledStudent Which version and OS are you in? I can't reproduce it on Wolfram Cloud, but the version number therein is "13.2.0 for Linux x86 (64-bit) (November 2, 2022)"(!). If I have to guess, this is a bug that's fixed in version 13.2, changing nsol[t, Sqrt[x^2 + y^2], ArcTan[x, y]] to If[R1^2 < x^2 + y^2 < R2^2, nsol[t, Sqrt[x^2 + y^2], ArcTan[x, y]], 0.] should fix the problem in your version. Dec 6, 2022 at 15:13
• I am on the 13.1 version on windows 10. Implementing the changes you suggested did remove the errors, but now the peak from the Plot3D looks like this. Well I guess this is just a minor issue so in the end, thank you for the enlightening answer. Dec 6, 2022 at 15:35
• @ConfuzzledStudent Does e.g. adding PlotPoints->100 to Plot3D` help? Dec 6, 2022 at 15:42
• It definitely helped alright. Thanks. Dec 6, 2022 at 15:55