# Numerical Solution for a Non-Linear Functional Fractional Differential Equation (FFDE)

I tried solve non-linear Functional-Fractional Differential Equation (FFDE) with this method, but it works on only for range: $$x\in \{0,1\}$$.

I what extend the solution range for example for general a and b.

a =-3;
b = 2;
DSolve[{CaputoD[y[x], {x, 3/2}] - Sqrt[y[x]] - y[2 x- 1] == 0, y[x /; x <= 0] == 2, y'[x /; x <= 0] == 1}, y[x], {x, a, b}]


How to change the code to work for range $$x\in \{-3,2\}$$ and plot solution ?

• This is fractional delayed differential equation FDDE since we have y[2 x - 1]. In a theory we need initial data for x<0 as well to define function for x>0. Feb 11 at 12:24
• @AlexTrounev .Ok I update the code. Feb 11 at 12:30
• Initial data your given are not satisfier to equation. Feb 11 at 12:40
• @AlexTrounev Ok putt some random Initial data that satisfier to equation. Feb 11 at 12:44
• It means that you are looking for solution at 2>=x>0 only since at x<0 we define function with first derivative. Also please pay attention that in a new redaction you put y[x /; x <= 0] == 2, y'[x /; x <= 0] == 1 that is apparently wrong. :) Feb 11 at 12:53

This is fractional delayed differential equation (FDDE), therefore we can't use method described in my answer here. Instead of the Haar wavelets method we try to use implicit difference scheme described in this paper. Note that for this type of equations we need to define $$y(x)$$ at $$x\le 0$$ and at $$2x-1\ge a$$ as well. Without this data we can't compute an unique solution. While for $$x\le 0$$ we have $$y(x)=2, y'(x)=0$$ (it could not be y'(x)=1 as given above), for $$2x-1\ge 2$$ we have no data. So, we put for example $$f(2 x-1)=0, 2 x-1\ge 2$$. To reduce the Caputo derivative order we introduce a new variable $$z=y'(x)$$, then our code is follows

L = 2; dx = L/100; xgrid = Range[0, L, dx]; nn = Length[xgrid];
M1 = NDSolveFiniteDifferenceDerivative[Derivative[1], xgrid,
DifferenceOrder -> 2]@"DifferentiationMatrix";
Z = Array[zz, {nn}]; Y = Array[yy, {nn}]; y0 = 2; dy0 = 0;
eq1 = Drop[Z - M1 . Y, 1]; bc = {Z[[1]] - dy0, Y[[1]] - y0};

a[j_][q_] := (j + 1)^(1 - q) - j^(1 - q);
b[j_][q_] := ((j + 1)^(2 - q) - j^(2 - q))/(2 -
q) - ((j + 1)^(1 - q) + j^(1 - q))/2;

eq2 = Table[0, {nn - 1}]; q = 1/2; Do[
c[0] = If[n == 1, 1, a[0][q] + b[0][q]];
If[n >= 2,
Do[c[j] = a[j][q] + b[j][q] - b[j - 1][q];, {j, 1, n - 2}];
c[n - 1] = a[n - 1][q] - b[n - 2][q]];
eq2[[n - 1]] =
dx^(-q)/Gamma[2 - q] (c[0] Z[[n]] - c[n - 1] Z[[1]] -
Sum[(c[n - j - 1] - c[n - j]) Z[[j]], {j, 1, n - 1}]) -
Sqrt[Y[[n]]] -
If[2 xgrid[[n]] - 1 <= 0, y0,
If[Round[(2 xgrid[[n]] - 1)/dx] >= nn + 1, 0,
Y[[Round[(2 xgrid[[n]] - 1)/dx]]]]];, {n, 2, nn}];

varALL = Join[Y, Z]; eqAll = Join[eq1, eq2, bc];
{res, {evx}} =
Reap[FindRoot[Table[eqAll[[i]] == 0, {i, Length[eqAll]}],
Table[{varALL[[i]], .1}, {i, Length[varALL]}],
Method -> {"Newton", "StepControl" -> "TrustRegion"},
MaxIterations -> 1000, EvaluationMonitor :> Sow[varALL]]];


Visualization

ListLinePlot[Transpose[{xgrid, Y /. res}], AxesLabel -> {"x", "y"}]

`

• This wasn't what I expected. If I had a different equation, I would have to redo the entire code each time. Maybe it is possible to modify the code to be for general equation. and general initial conditions This The paper costs 39.95 euros quite expensive for me .Thanks positive answer (+1) Feb 12 at 9:50
• @MariuszIwaniuk Probably it is possible to modify code so that it could be applicable for general equation, but even for this simpler case we need some analysis to understand the function behavior. Feb 12 at 11:20