12

We are familiar that the Predictor-Corrector Algorithm is applied to solve the Fractional-Order Differential systems. I need to solve the following Fractional Order Delay Differential System. Please share any MATHEMATICA Code for the following Fractional Order Delay Differential System:

D^α x[t] == z[t] + (y[t - τ] - a)*x[t];
D^α y[t] == 1 - b*y[t] - (x[t - τ])^2;
D^α z[t] == -x[t - τ] - c*z[t];

Where,

 a = 3; b = 0.1; c = 1, τ = 0.35 , and 
x[0] = 0.1; y[0] = 4; z[0] = 0.5 and α = 0.90.
Mariusz Iwaniuk
  • 13,841
  • 1
  • 25
  • 41
vicky
  • 353
  • 1
  • 11
  • 1
    vicky, you are aware that what you posted is not Mathematica code. Can you edit your question to make this code valid? – user21 Feb 09 '21 at 08:37
  • For the time delay problem we need to define x[t], y[t], z[t] for $t<0$ as well. – Alex Trounev Feb 09 '21 at 13:20
  • 1
    @AlexTrounev yes here i have mentioned those conditions too. Considered fractional order delay differential system,
    $\begin{eqnarray} D^{\alpha} x(t)&=&z+(y (t - \tau) - a)x\ D^{\alpha} y (t) &=& 1 - b y - (x (t - \tau))^2\ D^{\alpha} z (t) &= &-x (t - \tau) - c z.\ \end{eqnarray}$

    Where, $ a = 3,~ b =0.1,~ c = 1,~ \tau = 0.35 $, and $x(0) = 0.1, ~y(0) =4,~ z(0) = 0.5$ and $\alpha = 0.90.$ and for, $t<0,~ x(t)=0.1, ~y(t)=4,~ z(t)=0.5$

     – vicky Feb 09 '21 at 13:43
  • I found some Mathematica code in this forum but all those code use NDSolve built-in function. So dear all if possible please share me any link that contains code which applies numerical analysis to solve a fractional order delay differential system or even an integer order delay differential system is fine. – vicky Feb 09 '21 at 13:53
  • 2
    @vicky Is time interval up to infinity or limited? Please, pay attention that here is Mathematica forum, not Python. So function notation is f[t] not f(t ) :) – Alex Trounev Feb 09 '21 at 13:55
  • @AlexTrounev yes time interval is finite. For example time can be considered till 2000 – vicky Feb 09 '21 at 13:59
  • @vicky Why interval is required to be so large? There are several definitions of the fractional derivatives by Caputo, Riemann-Liouville, Atangana–Baleanu. What $D^{\alpha}$ means in your code? – Alex Trounev Feb 09 '21 at 15:48
  • @AlexTrounev I planned to study the chaotic nature of the above fractional order delay system so the interval should be at least 500 because chaos exist at large time. I adopt Caputo derivative. – vicky Feb 09 '21 at 15:54
  • @AlexTrounev In Mathematica I have written code for Predictor-Corrector Algorithm, which works good for a fractional order differential system but i am not able to edit that code for a delay system. Could you please edit my code so that it could work for a delay system. – vicky Feb 09 '21 at 15:59
  • 2
    @vicky Maybe, I can, but where is this code? – Alex Trounev Feb 09 '21 at 16:02
  • Crossposted here. – Rohit Namjoshi Feb 09 '21 at 17:05
  • @AlexTrounev I am waiting for your code so could you please post the code. – vicky Feb 10 '21 at 17:45
  • @vicky Ah, sorry! See my answer. – Alex Trounev Feb 10 '21 at 21:33

1 Answers1

15

There is no theorem for the fractional delay system like for the fractional system discussed in the paper Analysis of Fractional Differential Equations and used to solve it with using the predictor-corrector method in this paper. But we can suggest that we can use this method at least for $0\le t \le \tau $, since functions $x(t-\tau),y(t-\tau)$ are given as initial condition. On this step we compute x[t],y[t],z[t], and hence we can make next step on $\tau \le t \le 2\tau $. Repeating this step we can solve problem. But final code very slow, and numerical solution slightly differ from that we can get with NDSolve for $\alpha=1$ (see code below and pictures after)

a1 = 3; b1 = 0.1; 
c1 = 1; 
f[t_, x_, y_, z_] := x*(y - a1) + z; 
g[t_, x_, y_, z_] := -(b1*y) - x^2 + 1; 
d[t_, x_, y_, z_] := -(c1*z) - x; 
\[Alpha] =1; 
h = 0.01; 
x0 = 0.1; 
y0 = 4; 
z0 = 0.5; tau = 0.15; nn = Round[tau/h]; n = 40 nn; Do[x[i] = x0; y[i] = y0; z[i] = z0, {i, -nn, 0}]; 
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); ]; 
Do[Do[x1[i] = x[i - nn]; y1[i] = y[i - nn]; , {i, 0, s*nn}]; For[j = 1, j <= s*nn, j++, p[j] = (h^\[Alpha]*Sum[b[j - A]*f[A*h, x[A], y[A], z[A]], {A, 0, j - 1}])/Gamma[\[Alpha] + 1] + x[0]; 
     l[j] = (h^\[Alpha]*Sum[b[j - B]*g[B*h, x1[B], y[B], z[B]], {B, 0, j - 1}])/Gamma[\[Alpha] + 1] + y[0]; r[j] = (h^\[Alpha]*Sum[b[j - C]*d[C*h, x1[C], y[C], z[C]], {C, 0, j - 1}])/Gamma[\[Alpha] + 1] + z[0]; 
     x[j] = (h^\[Alpha]*(Sum[a[j - K]*f[h*K, x[K], y1[K], z[K]], {K, 1, j - 1}] + f[h*j, p[j], l[j], r[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*f[0, x[0], y[0], z[0]]))/Gamma[\[Alpha] + 2] + x[0]; 
     y[j] = (h^\[Alpha]*(Sum[a[j - F]*g[F*h, x1[F], y[F], z[F]], {F, 1, j - 1}] + g[h*j, p[j], l[j], r[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*g[0, x[0], y[0], z[0]]))/Gamma[\[Alpha] + 2] + y[0]; 
     z[j] = (h^\[Alpha]*(Sum[a[j - G]*d[G*h, x1[G], y[G], z[G]], {G, 1, j - 1}] + d[h*j, p[j], l[j], r[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*d[0, x[0], y[0], z[0]]))/Gamma[\[Alpha] + 2] + z[0]; ]; , 
  {s, 1, 40}]
lst1 = Table[{x[j], y[j], z[j]}, {j, n}]; lst = 
 Table[{x[j], y[j], z[j]}, {j, 10, n, 10}]; lst2 = 
 Table[{x[j], y[j]}, {j, 10, n, 10}];

Now we can compare this solution (red points) with NDSolve[] solution (solid lines)

Clear[x, y, z];

sol = NDSolve[{x'[t] == f[t, x[t], y[t - tau], z[t]], 
    y'[t] == g[t, x[t - tau], y[t], z[t]], 
    z'[t] == d[t, x[t - tau], y[t], z[t]], x[t /; t <= 0] == 0.1, 
    y[t /; t <= 0] == 4, z[t /; t <= 0] == 0.5}, {x, y, z}, {t, 0, 
    40}];


{Show[ParametricPlot3D[{x[t], y[t], z[t]} /. sol[[1]], {t, .0, 6}, 
   AxesLabel -> {x, y, z}], ListPointPlot3D[lst1, PlotStyle -> Red]], 
 Show[ParametricPlot[{x[t], y[t]} /. sol[[1]], {t, .0, 6}, 
   Frame -> True, Axes -> False, FrameLabel -> {x, y}], 
  ListPlot[lst2, PlotStyle -> Red]]}

Figure 1

We can improve this code and make it running faster but it is not so precise algorithm with of $h^2$ error. When we compare with NDSolve we need to decries h down to .0025 for t=15, but solutions are differ since it unstable for this set of parameters

a1 = 3; b1 = 0.1; 
c1 = 1; 
f[t_, x_, y_, z_] := x*(y - a1) + z; 
g[t_, x_, y_, z_] := -(b1*y) - x^2 + 1; 
d[t_, x_, y_, z_] := -(c1*z) - x; 
\[Alpha] =1; 
h = 0.0025; 
x0 = .1; 
y0 =4; 
z0 =1/2; tau = 0.15; nn = Round[tau/h]; smax = Round[15/tau];n=nn*smax; 
Do[x[i] = x0; y[i] = y0; z[i] = z0, {i, -nn, 0}]; 
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); ]; 
Do[Do[x1[i] = x[i - nn]; y1[i] = y[i - nn]; , {i, 0, s*nn}]; For[j = (s-1)*nn+1, j <= s*nn, j++, p[j] = (h^\[Alpha]*Sum[b[j - A]*f[A*h, x[A], y[A], z[A]], {A, 0, j - 1}])/Gamma[\[Alpha] + 1] + x[0]; 
     l[j] = (h^\[Alpha]*Sum[b[j - B]*g[B*h, x1[B], y[B], z[B]], {B, 0, j - 1}])/Gamma[\[Alpha] + 1] + y[0]; r[j] = (h^\[Alpha]*Sum[b[j - C]*d[C*h, x1[C], y[C], z[C]], {C, 0, j - 1}])/Gamma[\[Alpha] + 1] + z[0]; 
     x[j] = (h^\[Alpha]*(Sum[a[j - K]*f[h*K, x[K], y1[K], z[K]], {K, 1, j - 1}] + f[h*j, p[j], l[j], r[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*f[0, x[0], y[0], z[0]]))/Gamma[\[Alpha] + 2] + x[0]; 
     y[j] = (h^\[Alpha]*(Sum[a[j - F]*g[F*h, x1[F], y[F], z[F]], {F, 1, j - 1}] + g[h*j, p[j], l[j], r[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*g[0, x[0], y[0], z[0]]))/Gamma[\[Alpha] + 2] + y[0]; 
     z[j] = (h^\[Alpha]*(Sum[a[j - G]*d[G*h, x1[G], y[G], z[G]], {G, 1, j - 1}] + d[h*j, p[j], l[j], r[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*d[0, x[0], y[0], z[0]]))/Gamma[\[Alpha] + 2] + z[0]; ]; , 
  {s, 1, smax}]

lst1 = Table[{x[j], y[j], z[j]}, {j, n}]; lst = 
 Table[{x[j], y[j], z[j]}, {j, 10, n, 10}]; lst2 = 
 Table[{x[j], y[j]}, {j, n}];


Clear[x, y, z];


sol = NDSolve[{x'[t] == f[t, x[t], y[t - tau], z[t]], 
    y'[t] == g[t, x[t - tau], y[t], z[t]], 
    z'[t] == d[t, x[t - tau], y[t], z[t]], x[t /; t <= 0] == 0.1, 
    y[t /; t <= 0] == 4, z[t /; t <= 0] == 0.5}, {x, y, z}, {t, 0, 15}];


{Show[ParametricPlot3D[{x[t], y[t], z[t]} /. sol[[1]], {t, .0, 15}, 
   AxesLabel -> {x, y, z}], ListPointPlot3D[lst1, PlotStyle -> Red]], 
 Show[ListPlot[lst2, PlotStyle -> Red, Frame -> True, Axes -> False, 
   FrameLabel -> {x, y}, AspectRatio -> 1], 
  ParametricPlot[{x[t], y[t]} /. sol[[1]], {t, .0, 15}]]}

Figure 2

Alex Trounev
  • 44,369
  • 3
  • 48
  • 106
  • The above code works good for a three dimensional system but in the case of four dimensional system. I get the following error message General::ovfl: Overflow occurred in computation. Could anyone please give suggestions to solve this issue – vicky Feb 20 '21 at 12:14
  • @vicky We can give you suggestion if you show your system. – Alex Trounev Feb 20 '21 at 16:55
  • Here is my system $D^\alpha x[t]=x[t](y[t-\tau]-a1)+(m_1 w[t])+z[t];\ D^\alpha y[t]= - b1y[t]+(m_2 w[t])-x[t-\tau]^2+1;\ D^\alpha z[t] = -x[t - \tau] - c1z[t]+(m_3 w[t]);\ D^\alpha w[t] = -x[t]y[t]z[t].\$ $\alpha=0.93; a1 = 2.1; b1 = 0.01; c1 = 2.6; m_1 = 8.4; m_2 = 6.4; m_3= 2.2$ and for $t\leq 0,x[t] = 2;y[t] = 0;z[t] = -0.2;w[t] = 0.4.$ – vicky Feb 20 '21 at 17:00
  • @vicky Can you start a new topic with explanation of this system and with code you run? – Alex Trounev Feb 20 '21 at 17:10
  • In the above four dimensional system when $\tau=0.10$ and $n=4500$(smax=300 in the above code) I get error message – vicky Feb 20 '21 at 17:13
  • Sure will do it – vicky Feb 20 '21 at 17:14
  • Here I share my new question https://mathematica.stackexchange.com/questions/240347/overflow-during-computation – vicky Feb 20 '21 at 19:14