2
$\begingroup$

The following code display the Predictor-Corrector approach to solve the four dimensional Fractional-Order Delay Differential system,

The below code works fine for some values of $\tau$ but it gives an overflow error message for other values of $\tau$ for example, when $\tau=0.10$ and n=4500(smax=300 in the code) I get error message which says that an overflow occurred during computation.

a1 = 2.1; b1 = 0.01;
c1 = 2.6;
Subscript[m, 1] = 8.4;
Subscript[m, 2] = 6.4;
Subscript[m, 3] = 2.2;
f[t_, x_, y_, z_, w_] := x*(y - a1) + (Subscript[m, 1]* w) + z;
g[t_, x_, y_, z_, w_] := -(b1*y) + (Subscript[m, 2]* w) - x^2 + 1;
d[t_, x_, y_, z_, w_] := -(c1*z) + (Subscript[m, 3]* w) - x;
q[t_, x_, y_, z_, w_] := -x*y*z;
\[Alpha] = 0.93;
h = 0.01;
t[0] = 0;
x0 = 2;
y0 = 0;
z0 = -0.2;
w0 = 0.4;
tau = 0.10; nn = Round[tau/h]; smax = 300; n = nn*smax;
Do[x[i] = x0; y[i] = y0; z[i] = z0; w[i] = w0;, {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++, t[j] = j*h; 
  p[j] = (h^\[Alpha]*
       Sum[b[j - A]*f[A*h, x[A], y1[A], z[A], w[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], w[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], w[C]], {C, 0, j - 1}])/
     Gamma[\[Alpha] + 1] + z[0];
  o[j] = (h^\[Alpha]*
       Sum[b[j - H]*q[H*h, x[H], y[H], z[H], w[H]], {H, 0, j - 1}])/
     Gamma[\[Alpha] + 1] + w[0];
  x[j] = (h^\[Alpha]*(Sum[
          a[j - K]*f[h*K, x[K], y1[K], z[K], w[K]], {K, 1, j - 1}] + 
         f[h*j, p[j], l[j], r[j], 
          o[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
             j^\[Alpha])*f[0, x[0], y[0], z[0], w[0]]))/
     Gamma[\[Alpha] + 2] + x[0];
  y[j] = (h^\[Alpha]*(Sum[
          a[j - F]*g[F*h, x1[F], y[F], z[F], w[F]], {F, 1, j - 1}] + 
         g[h*j, p[j], l[j], r[j], 
          o[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
             j^\[Alpha])*g[0, x[0], y[0], z[0], w[0]]))/
     Gamma[\[Alpha] + 2] + y[0];
  z[j] = (h^\[Alpha]*(Sum[
          a[j - G]*d[G*h, x1[G], y[G], z[G], w[G]], {G, 1, j - 1}] + 
         d[h*j, p[j], l[j], r[j], 
          
          o[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
             j^\[Alpha])*d[0, x[0], y[0], z[0], w[0]]))/
     Gamma[\[Alpha] + 2] + z[0];
  w[j] = (h^\[Alpha]*(Sum[
          a[j - R]*q[R*h, x1[R], y[R], z[R], w[R]], {R, 1, j - 1}] + 
         q[h*j, p[j], l[j], r[j], 
          o[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
             j^\[Alpha])*q[0, x[0], y[0], z[0], w[0]]))/
     Gamma[\[Alpha] + 2] + w[0];];, {s, 1, smax}]
lst = Table[{t[i], x[i]}, {i, n}];
lst1 = Table[{x[i], y[i], z[i]}, {i, n}];
ListLinePlot[lst, PlotRange -> Full, 
  Frame -> True, PlotLabel -> "\[Alpha]=0.93,\[Tau]=0.10", 
 FrameStyle -> Black, 
  FrameLabel -> {t, x}, LabelStyle -> Directive[Blue, Bold], 
  PlotStyle -> Purple]
 ListPointPlot3D[lst1, PlotRange -> Full, 
  AxesLabel -> {x, y, z}, PlotLabel -> "\[Alpha]=0.93,\[Tau]=0.10", 
  LabelStyle -> Directive[Black, Bold, FontFamily -> "Helvetica"], 
  TicksStyle -> Black, PlotStyle -> Purple]

Could anyone please solve this issue.

$\endgroup$
2
  • $\begingroup$ The four dimensional system which I considered $ D^\alpha x[t]=x[t]*(y[t-\tau]-a1)+(m_1 * w[t])+z[t];\\ D^\alpha y[t]= - b1*y[t]+(m_2 * w[t])-x[t-\tau]^2+1;\\ D^\alpha z[t] = -x[t - \tau] - c1*z[t]+(m_3 * w[t]);\\ D^\alpha w[t] = -x[t]*y[t]*z[t].\\$ $\alpha\in(0,1);\tau>0; 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.$ $\endgroup$ – vicky Feb 20 at 19:08
  • $\begingroup$ Did you test this code for $\alpha =1$? $\endgroup$ – Alex Trounev Feb 20 at 20:14
2
$\begingroup$

We can use last code from my answer here extended with taken account new variable and functions. First, we test code for $alpha=1$ to compare with NDSolve[]

a1 = 2.1; b1 = 0.01; 
c1 = 2.6; m1 = 8.4; m2 = 6.4; m3 = 2.2; 
f1[t_, x_, y_, z_, w_] := x*(y - a1) + z + m1*w; 
f2[t_, x_, y_, z_, w_] := -(b1*y) - x^2 + 1 + m2*w; 
f3[t_, x_, y_, z_, w_] := -(c1*z) - x + m3*w; 
f4[t_, x_, y_, z_, w_] := (-x)*y*z
\[Alpha] = 1; 
h = 0.01; 
x0 = 2; 
y0 = 0; 
z0 = -0.2; w0 = 0.4; tau = 0.1; tmax = 5; nn = Round[tau/h]; smax = Round[tmax/tau]; 
  n = nn*smax; Do[x[i] = x0; y[i] = y0; z[i] = z0; w[i] = w0; , {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++, 
    p1[j] = (h^\[Alpha]*Sum[b[j - A]*f1[A*h, x[A], y1[A], z[A], w[A]], {A, 0, j - 1}])/
        Gamma[\[Alpha] + 1] + x[0]; 
     p2[j] = (h^\[Alpha]*Sum[b[j - B]*f2[B*h, x1[B], y[B], z[B], w[B]], {B, 0, j - 1}])/
        Gamma[\[Alpha] + 1] + y[0]; 
     p3[j] = (h^\[Alpha]*Sum[b[j - C1]*f3[C1*h, x1[C1], y[C1], z[C1], w[C1]], {C1, 0, j - 1}])/
        Gamma[\[Alpha] + 1] + z[0]; 
     p4[j] = (h^\[Alpha]*Sum[b[j - C2]*f4[C2*h, x[C2], y[C2], z[C2], w[C2]], {C2, 0, j - 1}])/
        Gamma[\[Alpha] + 1] + w[0]; 
     x[j] = (1/Gamma[\[Alpha] + 2])*(h^\[Alpha]*(Sum[a[j - K]*f1[h*K, x[K], y1[K], z[K], w[K]], 
           {K, 1, j - 1}] + f1[h*j, p1[j], p2[j], p3[j], p4[j]] + 
          ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*f1[0, x[0], y[0], z[0], w[0]])) + 
       x[0]; y[j] = (1/Gamma[\[Alpha] + 2])*h^\[Alpha]*
        (Sum[a[j - K]*f2[h*K, x1[K], y[K], z[K], w[K]], {K, 1, j - 1}] + 
         f2[h*j, p1[j], p2[j], p3[j], p4[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*
          f2[0, x[0], y[0], z[0], w[0]]) + y[0]; 
     z[j] = (1/Gamma[\[Alpha] + 2])*h^\[Alpha]*(Sum[a[j - K]*f3[h*K, x1[K], y[K], z[K], w[K]], 
          {K, 1, j - 1}] + f3[h*j, p1[j], p2[j], p3[j], p4[j]] + 
         ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*f3[0, x[0], y[0], z[0], w[0]]) + z[0]; 
     w[j] = (1/Gamma[\[Alpha] + 2])*h^\[Alpha]*(Sum[a[j - K]*f4[h*K, x[K], y[K], z[K], w[K]], 
          {K, 1, j - 1}] + f4[h*j, p1[j], p2[j], p3[j], p4[j]] + 
         ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*f4[0, x[0], y[0], z[0], w[0]]) + 
       w[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}];

Numerical solution with NDSolve

Clear[x, y, z, w];

sol = NDSolve[{x'[t] == f1[t, x[t], y[t - tau], z[t], w[t]], 
    y'[t] == f2[t, x[t - tau], y[t], z[t], w[t]], 
    z'[t] == f3[t, x[t - tau], y[t], z[t], w[t]], 
    w'[t] == f4[t, x[t], y[t], z[t], w[t]], x[t /; t <= 0] == x0, 
    y[t /; t <= 0] == y0, z[t /; t <= 0] == z0, 
    w[t /; t <= 0] == w0}, {x, y, z, w}, {t, 0, tmax}];

Compare two solution

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

Figure 1

Second, we run code for $\alpha =0.93, smax=300$ and plot numerical solution

a1 = 2.1; b1 = 0.01; 
c1 = 2.6; m1 = 8.4; m2 = 6.4; m3 = 2.2; 
f1[t_, x_, y_, z_, w_] := x*(y - a1) + z + m1*w; 
f2[t_, x_, y_, z_, w_] := -(b1*y) - x^2 + 1 + m2*w; 
f3[t_, x_, y_, z_, w_] := -(c1*z) - x + m3*w; 
f4[t_, x_, y_, z_, w_] := (-x)*y*z
\[Alpha] = 0.93; 
h = 0.01; 
x0 = 2; 
y0 = 0; 
z0 = -0.2; w0 = 0.4; tau = 0.1; tmax = 30; nn = Round[tau/h]; smax = Round[tmax/tau]; 
  n = nn*smax; Do[x[i] = x0; y[i] = y0; z[i] = z0; w[i] = w0; , {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++, 
    p1[j] = (h^\[Alpha]*Sum[b[j - A]*f1[A*h, x[A], y1[A], z[A], w[A]], {A, 0, j - 1}])/
        Gamma[\[Alpha] + 1] + x[0]; 
     p2[j] = (h^\[Alpha]*Sum[b[j - B]*f2[B*h, x1[B], y[B], z[B], w[B]], {B, 0, j - 1}])/
        Gamma[\[Alpha] + 1] + y[0]; 
     p3[j] = (h^\[Alpha]*Sum[b[j - C1]*f3[C1*h, x1[C1], y[C1], z[C1], w[C1]], {C1, 0, j - 1}])/
        Gamma[\[Alpha] + 1] + z[0]; 
     p4[j] = (h^\[Alpha]*Sum[b[j - C2]*f4[C2*h, x[C2], y[C2], z[C2], w[C2]], {C2, 0, j - 1}])/
        Gamma[\[Alpha] + 1] + w[0]; 
     x[j] = (1/Gamma[\[Alpha] + 2])*(h^\[Alpha]*(Sum[a[j - K]*f1[h*K, x[K], y1[K], z[K], w[K]], 
           {K, 1, j - 1}] + f1[h*j, p1[j], p2[j], p3[j], p4[j]] + 
          ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*f1[0, x[0], y[0], z[0], w[0]])) + 
       x[0]; y[j] = (1/Gamma[\[Alpha] + 2])*h^\[Alpha]*
        (Sum[a[j - K]*f2[h*K, x1[K], y[K], z[K], w[K]], {K, 1, j - 1}] + 
         f2[h*j, p1[j], p2[j], p3[j], p4[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*
          f2[0, x[0], y[0], z[0], w[0]]) + y[0]; 
     z[j] = (1/Gamma[\[Alpha] + 2])*h^\[Alpha]*(Sum[a[j - K]*f3[h*K, x1[K], y[K], z[K], w[K]], 
          {K, 1, j - 1}] + f3[h*j, p1[j], p2[j], p3[j], p4[j]] + 
         ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*f3[0, x[0], y[0], z[0], w[0]]) + z[0]; 
     w[j] = (1/Gamma[\[Alpha] + 2])*h^\[Alpha]*(Sum[a[j - K]*f4[h*K, x[K], y[K], z[K], w[K]], 
          {K, 1, j - 1}] + f4[h*j, p1[j], p2[j], p3[j], p4[j]] + 
         ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])*f4[0, x[0], y[0], z[0], w[0]]) + 
       w[0]; ]; , {s, 1, smax}]

{ListPointPlot3D[Table[{x[j], y[j], z[j]}, {j, n}], PlotRange -> Full, 
  AxesLabel -> {"x", "y", "z"}],
 ListLinePlot[Table[{x[j], y[j]}, {j, n}], PlotRange -> Full, Frame -> True, 
  Axes -> False, FrameLabel -> {"x", "y"}, AspectRatio -> 1, 
  PlotLabel -> Row[{"\[Tau] = ", tau}]]}

Figure 2

$\endgroup$
2
  • $\begingroup$ yes the above code works for tau=0.1 but still overflow in computation occurs when tau takes values like 0.8 $\endgroup$ – vicky Feb 22 at 16:36
  • 1
    $\begingroup$ @vicky Please, run first test with $\alpha =1$, tmax=10 and compare with NDSolve. Pay attention for the big discrepancies and reduce h to 0.0025. Run test for tmax=15, pay attention that NDSolve can't solve this range and also there are discrepancies of two solutions for t<=10. Reduce $h$ to 0.001 and run test again. If we can't solve this problem with NDSolve at $\alpha =1$, then we probably can't solve it with this algorithm for $\alpha < 1$. $\endgroup$ – Alex Trounev Feb 22 at 17:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.