# Numerical solution for a non-linear Fractional Differential Equation (FDE)

As shown below, a neat explicit expression is obtained for F=2, however an exact solution is not present for 1< F < 2.

How do we obtain numerical values for F = 1.5 (for instance)?

There have been suggestions of approximate methods (Adomian decomposition method (ADM), the variational iteration method (VIM)), these could involve time-consuming data crunching, and not sure if such techniques are already utilized in Mathematica.

The exact solution for the F=2 case does allow the testing of these possible alternative methods, hopefully one where rapid convergence is noted. An algorithm based on the Adomian decomposition method and the Padé approximants technique for nonlinear fractional differential equations with initial or boundary conditions is described here, but I am unable to access this code.

F = 2;

sol = DSolve[{CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 0,
y[0] == 1, y'[0] == 2}, y[x], x]

Table[{x, Evaluate[y[x] /. sol[[1]]]}, {x, 0, 1, 0.1}]

Plot[Evaluate[y[x]  . sol[[1]]], {x, 0, 1}]

> (* {{y[x] ->     1/16 E^(-2 x) (E^x + 54 E^(2 x) + 81 E^(3 x) -
>       12 Sqrt[E^(3 x) (1 + 9 E^x)^2])}} *)
>
> (* {{0., 1.}, {0.1, 1.21698}, {0.2, 1.47099}, {0.3, 1.76705}, {0.4,
>   2.11074}, {0.5, 2.50829}, {0.6, 2.96662}, {0.7, 3.49344}, {0.8,
>   4.09733}, {0.9, 4.78782}, {1., 5.57553}} *)


• NDSolve can't deal with CaputoD yet ... Aug 5 at 13:00
• Fractional differential equation is a relatively new topic (according to my limited knowledge). Are you aware of any standard method for numerically solving the problem? If the answer is yes, mentioning it in the question or even changing the question to something like "how to implement ×× method for solving this FDE" will make the question more attractive. Aug 5 at 15:10
• "but I am unable to access this code" - you didn't look for it hard enough, I would say: data.mendeley.com/datasets/dn82gkvf4f/1 yesterday
• @J.M. many thanks, I got into this package, but its Maple based, will be absolutely bliss if this can be coded in Mathematica, and then compared with the awesome answers provided below, not just for this problem but for other generalized equations as well.. yesterday

## Update

Though my old answer is wrong, I manage to think out a correct one.

The idea is again simple: given that usually the solution of differential equation can be approximated (without Runge's phenomenon) with n-th order polynomial whose coefficient is determined by function values at Chebyshev-Gauss-Lobatto grid (notice this is essentially the spirit of pseudospectral method), and CaputoD is capable of calculating derivative of polynomial, we assume the solution of the initial value problem (IVP) in the whole interval $$[0,1]$$ can be approximated by

$$y=c_0+c_1 x+…+c_{n-1} x^{n-1}$$

n = 16;
range = Range[0, n - 1];
coef = c /@ range;
poly[x_] = coef . x^range;


Then substitute it into the differential equation and i.c.s, and calculate the specific algebraic equations at Chebyshev-Gauss-Lobatto grid:

CGLGrid[xl_, xr_, n_Integer /; n > 1] :=
1/2 (xl + xr + (xl - xr) Cos[(π Range[0, n - 1])/(n - 1)])
grid = CGLGrid[0, 1, n];

F = 3/2;
{eq, ic} = {CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 0,
{y[0] == 1, y'[0] == 2}};

{approxeq, approxic} = {eq, ic} /. y -> poly; // AbsoluteTiming
cguess = 1;
crule = FindRoot[
Table[approxeq, {x, grid[[3 ;;]]}]~Join~approxic, {#, cguess} & /@
coef]; // AbsoluteTiming

pseudo = Plot[poly[x] /. crule // Evaluate, {x, 0, 1},
PlotRange -> All, GridLines -> Automatic]


Now the solution is consistent with those shown by Alex.

BTW, the following is a compiled version of the predictor-corrector method shown in Alex's answer:

g[t_, y_] := 9/4 Sqrt[y] + y;

alexsolver = Hold@Compile[{{n, _Integer}, α}, Module[{
h = 1./n,
y = Table[0., {n + 1}],
y0 = Table[0., {n + 1}],
a = Table[0., {n + 1}],
b = Table[0., {n + 1}],
p = Table[0., {n + 1}], k, j},
y[0] = y0[0] = 1;
y0[1] = 2;
For[k = 1, k <= n, k++, b[k] = k^α - (k - 1)^α;
a[
k] = -(2*k^(α + 1)) + (k - 1)^(α + 1) + (k + 1)^(α +
1);];
For[j = 1, j <= n, j++,
p[j] = (h^α*Sum[b[j - A]*g[A*h, y[A]], {A, 0, j - 1}])/
Gamma[α + 1] + y0[0] + j h y0[1];
y[
j] = (h^α*(Sum[a[j - f]*g[f*h, y[f]], {f, 1, j - 1}] +
g[h*j, p[j]] + ((j - 1)^(α + 1) - (-α + j - 1)*
j^α)*g[0, y[0]]))/Gamma[α + 2] + y0[0] + j h y0[1]];
y], CompilationTarget -> "C", RuntimeOptions -> "Speed"] /.
DownValues@g /. (head : y | y0 | a | b | p)[i_] :> head[[i + 1]] // ReleaseHold;

n = 10^4; h = 1/n; F = 3/2; ylst = alexsolver[n, F]; // AbsoluteTiming

alexsol = ListInterpolation[ylst, {0, 1}];
pseudo~Show~ListLinePlot[alexsol, PlotStyle -> {Red, Dashed}]


The following answer is wrong, because definition of CaputoD involves integration from 0 to x, so we can not approximate CaputoD only with function value in a small interval $$[x_0,x_0+dx]$$ (except for $$x_0=0$$).

Since currently NDSolve cannot handle fractional differential equation, I'd like to add a primary but working home-made solution.

The idea is simple: given that usually the solution of differential equation can be approximated with piecewise polynomial (notice this is essentially the spirit of finite difference method (FDM)), and CaputoD is capable of calculating derivative of polynomial, we assume that the solution of the initial value problem (IVP) in a small interval $$[x_0,x_0+dx]$$ can be approximated by

$$y\approx a+bx+cx^2$$

yfunc = x |-> a + b x + c x^2;


Substitue this into the IVP and solve:

ivp = {CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 0 /. x -> x0 + dx,
y[x0] == y0, y'[x0] == dy0};
solrule = Solve[ivp /. y -> yfunc, {a, b, c}]; // AbsoluteTiming


Solve::nongen There may be values of the parameters for which some or all solutions are not valid.

(* {0.0378036, Null} *)


Hmm… something improper seems to happen in Solve. Further check via Dimensions@solrule suggests 2 solutions are found, which is a bit counterinituive. I don't bother to looking into what happens in Solve, so I simply check the solutions numerically:

(Subtract @@@ ivp /. y -> yfunc /. # /. {x0 -> 0, y0 -> 1, dy0 -> 2,
dx -> 0.1, a -> 1, b -> 2} // N) & /@ solrule
(* {{-4.89298, 0., 0.}, {-4.44089*10^-16, 0., 0.}} *)


OK, seems that only the last solution is valid. Extract it, and substitute it into yfunc, we obtain the approximate solution at $$x_0+dx$$. Needless to say, the process can be done repeatively, so we can obtain the numeric solution in arbitrary interval. The programming is straightforward:

cf = Compile[{x0, y0, dy0, dx},
Evaluate[{x0 + dx, yfunc[x0 + dx], yfunc'[x0 + dx], dx} /. solrule[[-1]]]]

solver[x0_, y0_, dy0_, dx_, tend_] :=
MapAt[List,
NestList[cf @@ # &, {x0, y0, dy0, dx}, tend/dx // Round]\[Transpose] // Most,
{1, All}] // Transpose // Interpolation


Here, 4th syntax in document of Interpolation is used to take the derivative into account when building InterpolatingFunction.

Usage of solver:

tend = 1;
nsol = solver[0, 1, 2, 0.1, tend]
nsol // ListPlot


I've used an undocumented syntax of ListPlot here. See this post for more info.

Finally let's check the convergency of this method:

With[{lst = 10.^-Range[5]},
Table[solver[0, 1, 2, dx, tend], {dx, lst}] //
ListLinePlot[#, PlotLegends -> lst] &]


Um… not too bad (at least in my view). As mentioned in the beginning, this is merely a primary solution. Seeing it working, I'm satisfied enough.

• Nice approach (+1). Note, we can test your algorithm for CaputoD[y[x], {x, F}] - y[x] == 0 with F={1/2, 1}. Also there are wavelets solution on mathematica.stackexchange.com/questions/221609/… and predictor-corrector solution on mathematica.stackexchange.com/questions/239669/… Aug 6 at 0:01
• I have compared 3 methods in my answer, and your method seems not so good. Aug 6 at 17:21
• @xzczd Yes correct about the residual error being insignificant (Trounev's method) , also issues with smaller F ( less than 1.5) in your method where the simple series is taken, perhaps the coefficients to be used needs further examination. In any case both suggested answers are good steps forward. 2 days ago
• @xzczd Do you mean predictor-corrector method? Actually it is not mine method, it taken from the NASA report published on ntrs.nasa.gov/api/citations/20020024453/downloads/… . Also I don't understand what do you try to compare with? 2 days ago
• @AlexTrounev I add a new method that agrees with yours. 2 days ago

There are several numerical methods to solve fractional DE and PDE, including predictor-corrector method described here, and wavelets method discussed in our paper and here. First we solve this problem for $$1\le F\le 2$$ with Haar wavelets. Note, that in this case we have the second order problem, and it is why we need two boundary conditions. In a case $$F=2$$ we have CaputoD[y[x], {x, F}]=D[y[x], {x, 2}], therefore problem turns to second order ODE with known solution. This solution we gonna use as a test. With Haar wavelets we can compute CaputoD[y[x], {x, q}] as follows

h[x_, k_, m_] :=
WaveletPsi[HaarWavelet[], m x - k, WorkingPrecision -> Infinity];

p[x_, k_, m_] :=
Piecewise[{{(1 + k - m*x)/m, k >= 0 && 1/m + (2*k)/m - 2*x < 0 &&
1/m + k/m - x >= 0 && m > 0}, {(-k + m*x)/m,
k >= 0 && 1/m + (2*k)/m - 2*x >= 0 &&
k/m - x < 0 && 1/m + k/m - x >= 0 && m > 0}}, 0];
h1[x_] := WaveletPhi[HaarWavelet[], x, WorkingPrecision -> Infinity];

p1[x_] := Piecewise[{{1, x > 1}}, x];

pc[t_, k_, m_, q_] :=
Piecewise[{{-(t^(1 - q)/(-1 + q)), k == 0 && 1/m - 2*t >= 0 &&
m > 0 && t > 0 && 1/m - t >= 0},
{-((m^(-1 + q)*(1/(-k + m*t))^(-1 + q))/(-1 + q)),
k > 0 && 1/m + (2*k)/m - 2*t > 0 && k/m - t < 0 && m > 0 &&
1/m + k/m - t > 0},
{(-t^q + 2*m*t^(1 + q) - m*t*(-(1/(2*m)) + t)^q)/
(t^q*(-(1/(2*m)) + t)^q*(m*(-1 + q))),
k == 0 && m > 0 && 1/m - 2*t < 0 && 1/m - t >= 0},
{(1/(-1 + q))*((2^(-1 + q)*m^(-1 + 2*q)*(-(-(k/m) + t)^q -
2*k*(-(k/m) + t)^q + 2*m*t*(-(k/m) + t)^q +
2*k*(-((1/2 + k)/m) + t)^q -
2*m*t*(-((1/2 + k)/m) + t)^
q))/((1 + 2*k - 2*m*t)*(k - m*t))^q),
k > 0 && 1/m + (2*k)/m - 2*t == 0 && m > 0 &&
1/m + k/m - t > 0},
{-((1/(-1 + q))*((2^(-1 + q)*m^(-1 + 2*q)*
(-2*(-((1/2 + k)/m) + t)^
q*((1 + 2*k - 2*m*t)*(k - m*t))^
q - 2*k*(-((1/2 + k)/m) + t)^q*
((1 + 2*k - 2*m*t)*(k - m*t))^q +

2*m*t*(-((1/2 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*
(k - m*t))^q + (-((1 + k)/m) + t)^q*
((1 + 2*k - 2*m*t)*(k - m*t))^q +

2*k*(-((1 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*(k - m*t))^
q - 2*m*t*(-((1 + k)/m) + t)^q*
((1 + 2*k - 2*m*t)*(k - m*t))^
q + (-(k/m) + t)^q*
((1 + 2*k - 2*m*t)*(1 + k - m*t))^q +

2*k*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
q -

2*m*t*(-(k/m) + t)^q*((1 + 2*k - 2*m*t)*(1 + k - m*t))^
q - 2*k*(-((1/2 + k)/m) + t)^q*
((1 + 2*k - 2*m*t)*(1 + k - m*t))^q +

2*m*t*(-((1/2 + k)/m) + t)^q*((1 + 2*k - 2*m*t)*
(1 + k - m*t))^
q))/(((1 + 2*k - 2*m*t)*(k - m*t))^q*
((1 + 2*k - 2*m*t)*(1 + k - m*t))^q))),
k > 0 && m > 0 && 1/m + (2*k)/m - 2*t <= 0 &&
1/m + k/m - t <= 0},
{-((1/(2*m*(-1 + q)))*((2^q*m^(2*q)*t^q*(-(1/m) + t)^q*
(-(1/(2*m)) + t)^q -
2^(1 + q)*m^(1 + 2*q)*t^(1 + q)*
(-(1/m) + t)^q*(-(1/(2*m)) + t)^q -
2^(1 + q)*m^(2*q)*
t^q*(-(1/(2*m)) + t)^(2*q) +
2^(1 + q)*m^(1 + 2*q)*
t^(1 + q)*(-(1/(2*m)) + t)^(2*q) +
t^q*((-1 + m*t)*(-1 + 2*m*t))^q - 2*m*t^(1 + q)*
((-1 + m*t)*(-1 + 2*m*t))^q +
2*m*t*(-(1/(2*m)) + t)^q*
((-1 + m*t)*(-1 + 2*m*t))^q)/(t^
q*(-(1/(2*m)) + t)^q*
((-1 + m*t)*(-1 + 2*m*t))^q))),
k == 0 && 1/m - 2*t < 0 && 1/m - t < 0 && m > 0},
{(1/(-1 + q))*((2^(-1 + q)*m^(-1 + q)*((-m^q)*(-(k/m) + t)^q -
2*k*m^q*(-(k/m) + t)^q +
2*m^(1 + q)*t*(-(k/m) + t)^q +
2*k*m^q*(-((1/2 + k)/m) + t)^q - 2*m^(1 + q)*t*
(-((1/2 + k)/m) + t)^
q - ((1 + 2*k - 2*m*t)*(k - m*t))^q*
(1/(-1 - 2*k + 2*m*t))^q -
2*k*((1 + 2*k - 2*m*t)*(k - m*t))^q*
(1/(-1 - 2*k + 2*m*t))^q +
2*m*t*((1 + 2*k - 2*m*t)*(k - m*t))^q*
(1/(-1 - 2*k + 2*m*t))^q))/((1 + 2*k -
2*m*t)*(k - m*t))^
q), 1/m + (2*k)/m - 2*t < 0 && k > 0 && m > 0 &&
1/m + k/m - t > 0}}, 0]

pc1[t_, q_] := Piecewise[{{-(t^(1 - q)/(-1 + q)), t <= 1}},
-(((-1 + t)^q*t + t^q - t^(1 + q))/((-1 + t)^q*t^q*(-1 + q)))];


Here p, p1 are integrals of h,h1, and pc, pc1 are Caputo derivatives of h, h1. With these functions we can compute numerical solution for $$0\le q\le 1$$ by reducing order FDE as

eqs={CaputoD[z[x],{x,q}] - 9/4 Sqrt[y[x]] - y[x] == 0, z[x]==y'[x],
y[0] == 1, z[0] == 2}


Finally we define numerical functions y[x],z[x] and first derivatives y1[x], z1[x], and compute eqs in collocation points xcol

J = 4; M = 2^J; dx = 1/(2*M); xl = Table[l dx, {l, 0, 2 M}]; xcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
z1[x_, q_] :=
Sum[v[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
v1 pc1[x, q];
z[x_] :=
Sum[v[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
v1 p1[x] + v0;
y1[x_] :=
Sum[u[i, j] h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
u1 h1[x];
y[x_] :=
Sum[u[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
u1 p1[x] + u0;
varM = Join[{v0, v1, u0, u1},
Flatten[Table[v[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]],
Flatten[Table[u[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]];
eq[q_] :=
Flatten[Table[{1/Gamma[1 - q] z1[x, q] - 9/4 Sqrt[Abs[y[x]]] - y[x] ==
0, y1[x] - z[x] == 0}, {x, xcol}]]; bc = {y[0] == 1, z[0] == 2};


The system of algebraic equations can be solved with FindRoot for several q

qlist = {99/100, 1/2, 1/3};

Do[
eqM[i] = Join[eq[i], bc];
solv[i] =
FindRoot[eqM[i], Table[{varM[[j]], 1/10}, {j, Length[varM]}],
MaxIterations -> 1000, WorkingPrecision -> 30];, {i,
qlist}]


Visualization together with test data

Do[lst[i] =
Table[{x , Evaluate[y[x] /. solv[i]]}, {x, 0, 1, .01}];, {i, qlist}]

test = {{0., 1.}, {0.1, 1.21698}, {0.2, 1.47099}, {0.3,
1.76705}, {0.4, 2.11074}, {0.5, 2.50829}, {0.6, 2.96662}, {0.7,
3.49344}, {0.8, 4.09733}, {0.9, 4.78782}, {1., 5.57553}};

Show[ListLinePlot[Table[lst[i], {i, qlist}], Frame -> True,
FrameLabel -> {"x", "y"},
PlotRange -> All,
PlotLegends ->
Table[Row[{"F = ", qlist[[i]] + 1}], {i, Length[qlist]}]],
ListPlot[test, PlotStyle -> Red]]


Predictor-corrector code is based on the pseudocode published in the end of paper

g[t_, y_] := 9/4 Sqrt[y] + y;

\[Alpha] = 199./100;
h = 1./100;
y[0] = y0[0] = 1;
y0[1] = 2; n = 100;
For[k = 1, k <= n, k++, b[k] = k^\[Alpha] - (k - 1)^\[Alpha];
a[k] = -(2*k^(\[Alpha] + 1)) + (k - 1)^(\[Alpha] + 1) + (k +
1)^(\[Alpha] + 1);];
For[j = 1, j <= n, j++,
p[j] = (h^\[Alpha]*Sum[b[j - A]*g[A*h, y[A]], {A, 0, j - 1}])/
Gamma[\[Alpha] + 1] + y0[0] + j h y0[1];

y[j] = (h^\[Alpha]*(Sum[a[j - f]*g[f*h, y[f]], {f, 1, j - 1}] +
g[h*j, p[
j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
j^\[Alpha])*g[0, y[0]]))/Gamma[\[Alpha] + 2] + y0[0] +
j h y0[1];
;];


Visualization together with test data

lst = Table[{j h, y[j]}, {j, n}];
test = {{0., 1.}, {0.1, 1.21698}, {0.2, 1.47099}, {0.3,
1.76705}, {0.4, 2.11074}, {0.5, 2.50829}, {0.6, 2.96662}, {0.7,
3.49344}, {0.8, 4.09733}, {0.9, 4.78782}, {1., 5.57553}}; Show[
ListLinePlot[lst], ListPlot[test, PlotStyle -> Red]]


Finally note, that solver proposed by xzczd is not so good compare to wavelets and predictor-corrector method and needs to be improved. Code is given by

F = 3/2; yfunc =
x |-> a + b x +
c x^2; ivp = {CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] == 0 /.
x -> x0 + dx, y[x0] == y0, y'[x0] == dy0};
solrule =
Solve[ivp /. y -> yfunc, {a, b,
c}]; (Subtract @@@ ivp /. y -> yfunc /. # /. {x0 -> 0, y0 -> 1,
dy0 -> 2, dx -> 0.1, a -> 1, b -> 2} // N) & /@ solrule; cf =
Compile[{x0, y0, dy0, dx},
Evaluate[{x0 + dx, yfunc[x0 + dx], yfunc'[x0 + dx], dx} /.
solrule[[-1]]]];

solver[x0_, y0_, dy0_, dx_, tend_] :=
MapAt[List,
NestList[cf @@ # &, {x0, y0, dy0, dx},
tend/dx // Round]\[Transpose] // Most, {1, All}] // Transpose //
Interpolation
tend = 1;
nsol = solver[0, 1, 2, 0.01, tend]
nsol // ListPlot


In the picture below shown Figure 1 with wavelets method results, predictor-corrector method for F=3/2 - dashed red line and nsol - doted line. We see that two methods are in a good agreement while nsol is out of range.

Updated code by xzczd is in a good agreement with Haar wavelets method and predictor-corrector method as well. Code

n = 16;
range = Range[0, n - 1];
coef = c /@ range;
poly[x_] = coef . x^range;
CGLGrid[xl_, xr_, n_Integer /; n > 1] :=
1/2 (xl + xr + (xl - xr) Cos[(\[Pi] Range[0, n - 1])/(n - 1)])
grid = CGLGrid[0, 1, n];

F = 3/2;
{eq, ic} = {CaputoD[y[x], {x, F}] - 9/4 Sqrt[y[x]] - y[x] ==
0, {y[0] == 1, y'[0] == 2}};

approx = eq /. y -> poly; // AbsoluteTiming
cguess = 1;
crule = FindRoot[
Table[approx, {x, grid[[3 ;;]]}]~
Join~(ic /. y -> poly), {#, cguess} & /@ coef]; // AbsoluteTiming

pseudo =
Plot[poly[x] /. crule // Evaluate, {x, 0, 1}, PlotRange -> All,
GridLines -> Automatic, PlotStyle -> Dashed]


In the picture below the dashed line is the numerical solution computed with xzczd code compare to wavelets (left) and predictor-corrector method (right):