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 ... Commented Aug 5, 2022 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. Commented Aug 5, 2022 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 Commented Aug 7, 2022 at 22:53

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):

• Can you extend the solution range? For example and for general: $x\in \{a,b\}$ maybe for : $x\in \{-2,3\}$. Thanks for excellent answer. (+1) Commented Feb 10 at 10:59
• Yes we can solve it on (a, b) using map (a, b)->(0,1) in a form $x=a+s (b-a)$ with $s\in \{0,1\}$ :) Commented Feb 10 at 11:52
• Can you update the code of yours last findings on the comment? because for me it's dark magic. Commented Feb 10 at 12:53
• @MariuszIwaniuk This answer is too long. Could you ask your question as a separate post? Commented Feb 10 at 17:07

Update 2: A general-purpose IVP solver for FDE

The following is a general-purpose implementation for the predictor-corrector method in this paper linked by Alex:

BeginPackage["fdeivp"];
fdesolve; fdesolvecore;
Begin["private"];
Clear[fdesolve, fdesolvecore];
Options[fdesolve] = {Compiled -> True};
fdesolve[α_, f_, y0_?NumericQ, rest__] := fdesolve[α, f, {y0}, rest]
fdesolve[α_, f_, y0_?VectorQ, rest__] /; Length@y0 != Ceiling@α :=
fdesolve[α, f, {y0}, rest]
fdesolve[α_, f_, y0_, T_, n_Integer : 25, OptionsPattern[]] :=
With[{size = Rest@Dimensions@y0},
fdesolvecore[α, f, size, OptionValue@Compiled][y0, T, n]]

fdesolvecore[α_, f_, size_, tf : True | False] :=
fdesolvecore[α, f, size, tf] =
With[{m = Ceiling[α], dim = 1 + Length@size},
Hold@function[argu, with[{h = T/n,
b = Table[If[k == 0, 0, k^α - (k - 1)^α], {k, 0, n}],
a = Table[If[k == 0, 0, (k + 1)^(α + 1) -
2 k^(α + 1) + (k - 1)^(α + 1)], {k, 0, n}]},
module[{y = Table[0., {n + 1}, ##], p},
y[0] = y0[0];
Do[p = Sum[(j h)^k/k! y0[k], {k, 0, m - 1}] +
h^α/Gamma[α + 1] Sum[b[j - k] f[k h, y[k]], {k, 0, j - 1}];
y[j] = Sum[(j h)^k/k! y0[k], {k, 0, m - 1}] +
h^α/Gamma[α + 2] (f[j h, p] + ((j - 1)^(α + 1) - (j - 1 - α) j^α) f[
0, y[0]] + Sum[a[j - k] f[k h, y[k]], {k, j - 1}]),
{j, n}]; y]]] & @@ (List /@ size) /. (head : y | y0 | a | b)[i_] :>
head[[i + 1]] /. {module -> Module, with -> With} /. {function ->
If[tf, Compile, Function],
argu -> If[tf, {{y0, _Real, dim}, {T, _Real}, {n, _Integer}}, {y0, T, n}]} /.
Factorial[n_] :> Gamma[1 + n] // ReleaseHold];
End[];
EndPackage[];


With this package, OP's problem can be solved as follows:

order = 1.5; rhsfunc = {t, y} |-> 9/4 Sqrt[y] + y; iclst = {1, 2}; tend = 1;

sollst = fdesolve[order, rhsfunc, iclst, tend];


The output is a List, it can be visualized with e.g.

ListPlot[sollst, DataRange -> {0, tend}]


If you prefer a function relationship as the output, build a interpolating function with e.g.

solfunc = ListInterpolation[sollst, {0, tend}]


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/… Commented Aug 6, 2022 at 0:01
• I have compared 3 methods in my answer, and your method seems not so good. Commented Aug 6, 2022 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. Commented Aug 7, 2022 at 3:35
• @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? Commented Aug 7, 2022 at 6:56
• @AlexTrounev I add a new method that agrees with yours. Commented Aug 7, 2022 at 9:04

Good reference is provided by "A review of the Adomian decomposition method and its applications to fractional differential equations ", Commun.Frac.Calc.3 (2) (2012) 73-99

Note that the non-linear operator is taken to be semi-analytic, and to simplify the computation we consider just the first four terms starting with y0 which is set to 1, in this code we identify the Adomian polynomials by A_0, A_1, etc and an image of the code is attached below for easy viewing of Subscripted terms

ClearAll;
Jf[u_, v_, x_] := Gamma[v + 1]/Gamma[u + v + 1] x^(u + v)
Subscript[y, 0] = 1;
Subscript[A, 0] = 1;
Subscript[y, 1][u_, x_] := 2 x + 9/4 Jf[u, 0, x] + Jf[u, 0, x]
Subscript[A, 1][u_, x_] := Subscript[y, 1][u, x]/(
2 Subscript[y, 0]^0.5)

Subscript[y, 2][u_, x_] := 17/4 Jf[u, 1, x] + (221/32)/Gamma[u + 1] Jf[u, u, x]

Subscript[A, 2][u_, x_] := Subscript[y, 2][u, x]/(2 Subscript[y, 0]^0.5) -  Subscript[y, 1][u, x]^2/(4 Subscript[y, 0]^(3/2))

Subscript[y, 3][u_,x_] := -(9/8) Jf[u, 2, x] + (17^2/(32*Gamma[u + 2]) - 117/(32*Gamma[u + 2])) Jf[u, 1 + u, x] + (3757/(256*Gamma[2 u + ]) - 1521/(512*Gamma[u + 2]^2)) Jf[
u, 2 u, x]

Subscript[A, 3][u_, x_] :=
Subscript[y, 3][u, x]/(2 Subscript[y, 0]^0.5) -
Subscript[y, 2][u, x]^2/(4 Subscript[y, 0]^(3/2)) - (
Subscript[y, 1][u, x] Subscript[y, 2][u, x])/(
2 Subscript[y, 0]^(3/2)) + (3 Subscript[y, 1][u, x]^2)/(
8 Subscript[y, 0]^(5/2))

Subscript[y, 4][u_, x_] :=
9/8 Jf[u, 3, x] +
1/16 ((-9*34)/(4*Gamma[u + 3]) - (18*17)/(4*Gamma[u + 2]) + (
27*13)/(4*Gamma[u + 1])) Jf[u, u + 2, x]

Subscript[y, T][u_, x_] :=
Subscript[y, 0] + Subscript[y, 1][u, x] + Subscript[y, 2][u, x] +
Subscript[y, 3][u, x]

Subscript[y, exact][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])

Plot[{Subscript[y, T][1.25, x], Subscript[y, T][1.5, x],
Subscript[y, T][1.75, x], Subscript[y, T][2, x]}, {x, 0, 1},
PlotTheme -> "Detailed", GridLines -> Automatic,
PlotLegends -> {1.25, 1.5, 1.75, 2, Exact}]

Plot[Subscript[y, T][2, x] - Subscript[y, exact][x], {x, 0, 1},
PlotTheme -> "Detailed", GridLines -> Automatic]

Table[{x, Subscript[y, T][2, x] - Subscript[y, exact][x]}, {x, 0, 1,
0.1}] // TableForm

Plot[{Subscript[y, T][1.5, x], poly[x] /. crule // Evaluate}, {x, 0,
1}, PlotTheme -> "Detailed", GridLines -> Automatic,

• ClearAll; doesn't have the effect in your mind. Please read the document of it carefully. Commented Aug 25, 2022 at 6:18