1
$\begingroup$

I have adopted the method in the link to solve a much simpler system but it takes ages to execute.

a1 = 2/3; b1 = 4/3; c1 = 1; d1 = 1; x0 = 0.9; y0 = 1.8;
\[Alpha] = 1; h = .0025;

tau = 0.15; nn = Round[tau/h]; smax = Round[15/tau]; n = nn*smax;

f[t_, x_, y_] := a1*x - b1*x*y;
g[t_, x_, y_] := c1*x*y - d1*y;

Do[x[i] = x0; y[i] = y0, {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]], {A, 0, j - 1}])/
Gamma[\[Alpha] + 1] + x[0];
l[j] = (h^\[Alpha]*
Sum[b[j - B]*g[B*h, x1[B], y[B]], {B, 0, j - 1}])/
Gamma[\[Alpha] + 1] + y[0];
x[j] = (h^\[Alpha]*(Sum[
a[j - K]*f[h*K, x[K], y1[K]], {K, 1, j - 1}] +
f[h*j, p[j],
l[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
j^\[Alpha])*f[0, x[0], y[0]]))/Gamma[\[Alpha] + 2] + x[0];
y[j] = (h^\[Alpha]*(Sum[
a[j - F]*g[F*h, x1[F], y[F]], {F, 1, j - 1}] +
g[h*j, p[j],
l[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
j^\[Alpha])*g[0, x[0], y[0]]))/Gamma[\[Alpha] + 2] +
y[0];];, {s, 1, smax}]

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

Clear[x, y];
sol = NDSolve[{x'[t] == f[t, x[t], y[t - tau]],
y'[t] == g[t, x[t - tau], y[t]], x[t /; t <= 0] == 0.9,
y[t /; t <= 0] == 1.8}, {x, y}, {t, 0, 40}];

{Show[ParametricPlot[{x[t], y[t]} /. sol[[1]], {t, .0, 20},
PlotRange -> All, Frame -> True, FrameLabel -> {x, y}],
ListPlot[lst1, PlotStyle -> Red, PlotRange -> All, Frame -> True]]}

Moreover for smaller steps and for alpha other than 1, the solution doesnt make anysense.

Improve this question
$\endgroup$
1
  • 3
    $\begingroup$ Could you add more written explanation of your problem? And what does "FDEs" stand for - fractional differential equations, functional differential equations, fuzzy differential equations, or something else? $\endgroup$
    – Chris K
    Jan 31 at 2:07

1 Answer 1

Reset to default
6
+50
$\begingroup$

There is a typo in the code. We should use y1[A] in p[j] definition. Also for $\alpha<1$ we should take very small step h, since predictor-corrector algorithm depends on $h^{\alpha}$. For example, if for $\alpha =1$ solution converges at $h=0.01$, then for $\alpha =1/2$ solution converges at $h=10^{-4}$. To test code we put $h=1/100, \alpha =1, 4/5$. As benchmark solution we use NDSolve at $\alpha=1$, so we have

Clear[x, y];
a1 = 2/3; b1 = 4/3; c1 = 1; d1 = 1; x0 = 0.9; y0 = 1.8;
\[Alpha] = 1; h = .01; par1 = h^\[Alpha]/Gamma[\[Alpha] + 1]; par2 = 
 h^\[Alpha]/Gamma[\[Alpha] + 2];

tau = 0.15; nn = Round[tau/h]; n = 40 nn;

f[t_, x_, y_] := a1*x - b1*x*y;
g[t_, x_, y_] := c1*x*y - d1*y;

Do[x[i] = x0; y[i] = y0, {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] = (par1 Sum[b[j - A]*f[A*h, x[A], y1[A]], {A, 0, j - 1}]) + 
     x[0];
   l[j] = (par1 Sum[b[j - B]*g[B*h, x1[B], y[B]], {B, 0, j - 1}]) + 
     y[0];
   x[j] = (par2 (Sum[a[j - K]*f[h*K, x[K], y1[K]], {K, 1, j - 1}] + 
         f[h*j, p[j], 
          l[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
             j^\[Alpha])*f[0, x[0], y[0]])) + x[0];
   y[j] = (par2 (Sum[a[j - F]*g[F*h, x1[F], y[F]], {F, 1, j - 1}] + 
         g[h*j, p[j], 
          l[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
             j^\[Alpha])*g[0, x[0], y[0]])) + y[0];];, {s, 1, 
   40}] // AbsoluteTiming


lst1 = Table[{x[j], y[j]}, {j, n}];

Clear[x, y];

\[Alpha] = 4/5; h = .01; par1 = h^\[Alpha]/Gamma[\[Alpha] + 1]; par2 =
  h^\[Alpha]/Gamma[\[Alpha] + 2];

tau = 0.15; nn = Round[tau/h]; n = 40 nn;

f[t_, x_, y_] := a1*x - b1*x*y;
g[t_, x_, y_] := c1*x*y - d1*y;

Do[x[i] = x0; y[i] = y0, {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] = (par1 Sum[b[j - A]*f[A*h, x[A], y1[A]], {A, 0, j - 1}]) + 
     x[0];
   l[j] = (par1 Sum[b[j - B]*g[B*h, x1[B], y[B]], {B, 0, j - 1}]) + 
     y[0];
   x[j] = (par2 (Sum[a[j - K]*f[h*K, x[K], y1[K]], {K, 1, j - 1}] + 
         f[h*j, p[j], 
          l[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
             j^\[Alpha])*f[0, x[0], y[0]])) + x[0];
   y[j] = (par2 (Sum[a[j - F]*g[F*h, x1[F], y[F]], {F, 1, j - 1}] + 
         g[h*j, p[j], 
          l[j]] + ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*
             j^\[Alpha])*g[0, x[0], y[0]])) + y[0];];, {s, 1, 
   40}] // AbsoluteTiming
lst2 = Table[{x[j], y[j]}, {j, n}];

Clear[x, y];

a1 = 2/3; b1 = 4/3; c1 = 1; d1 = 1; x0 = 0.9; y0 = 1.8; 
f[t_, x_, y_] := a1*x - b1*x*y;
g[t_, x_, y_] := c1*x*y - d1*y; tau = 0.15;
sol = NDSolve[{x'[t] == f[t, x[t], y[t - tau]], 
    y'[t] == g[t, x[t - tau], y[t]], x[t /; t <= 0] == 0.9, 
    y[t /; t <= 0] == 1.8}, {x, y}, {t, 0, 40}];


pl1 = ParametricPlot[{x[t], y[t]} /. sol[[1]], {t, .0, n h}, 
   PlotRange -> All, Frame -> True, FrameLabel -> {x, y}, 
   AspectRatio -> 1/2, PlotStyle -> Green];

Visualization

Show[pl1, ListLinePlot[lst1, PlotStyle -> {Red, Dashed}], 
 ListLinePlot[lst2, PlotStyle -> {Blue, Dashed}]]

Figure 1

With small modification of the code above we can compute up to t=15 as follows

a1 = 2/3; b1 = 
 4/3; c1 = 1; d1 = 1; x0 = 0.9; y0 = 1.8; \[Alpha] = 1; h = .01; par1 \
= h^\[Alpha]/Gamma[\[Alpha] + 1]; par2 = 
 h^\[Alpha]/Gamma[\[Alpha] + 2];

tau = 0.15; nn = Round[tau/h]; smax = 100; n = smax nn;

f[t_, x_, y_] := a1*x - b1*x*y;
g[t_, x_, y_] := c1*x*y - d1*y;

Do[x[i] = x0; y[i] = y0, {i, -nn, 0}]; f0 = f[0, x[0], y[0]]; g0 = 
 g[0, x[0], y[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); 
  c[k] = ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])];

Do[Do[x1[i] = x[i - nn]; y1[i] = y[i - nn];, {i, 0, s*nn}];
  For[j = 1, j <= s*nn, j++, bj = Table[b[j - A], {A, 0, j - 1}]; 
   f1j = Table[f[A*h, x[A], y1[A]], {A, 0, j - 1}]; 
   p[j] = par1 bj . f1j + x[0]; 
   g1j = Table[g[B*h, x1[B], y[B]], {B, 0, j - 1}];
   l[j] = par1 bj . g1j + y[0]; aj = Table[a[j - i], {i, 1, j - 1}]; 
   f2j = Table[f[h*Kk, x[Kk], y1[Kk]], {Kk, 1, j - 1}];
   x[j] = par2 (aj . f2j + f[h*j, p[j], l[j]] + c[j] f0) + x[0]; 
   g2j = Table[g[F*h, x1[F], y[F]], {F, 1, j - 1}];
   y[j] = 
    par2 (aj . g2j + g[h*j, p[j], l[j]] + c[j] g0) + y[0];];, {s, 1, 
   smax}] // AbsoluteTiming

This part with $\alpha=1$ takes 444s on my laptop.

lst1 = Table[{x[j], y[j]}, {j, n}];


Clear[x, y];

a1 = 2/3; b1 = 4/3; c1 = 1; d1 = 1; x0 = 0.9; y0 = 1.8; \[Alpha] = 
 4/5; h = .01; par1 = h^\[Alpha]/Gamma[\[Alpha] + 1]; par2 = 
 h^\[Alpha]/Gamma[\[Alpha] + 2];

tau = 0.15; nn = Round[tau/h]; n = smax nn;

f[t_, x_, y_] := a1*x - b1*x*y;
g[t_, x_, y_] := c1*x*y - d1*y;

Do[x[i] = x0; y[i] = y0, {i, -nn, 0}];
f0 = f[0, x[0], y[0]]; g0 = g[0, x[0], y[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); 
  c[k] = ((j - 1)^(\[Alpha] + 1) - (-\[Alpha] + j - 1)*j^\[Alpha])];

Do[Do[x1[i] = x[i - nn]; y1[i] = y[i - nn];, {i, 0, s*nn}];
  For[j = 1, j <= s*nn, j++, bj = Table[b[j - A], {A, 0, j - 1}]; 
   f1j = Table[f[A*h, x[A], y1[A]], {A, 0, j - 1}]; 
   p[j] = par1 bj . f1j + x[0]; 
   g1j = Table[g[B*h, x1[B], y[B]], {B, 0, j - 1}];
   l[j] = par1 bj . g1j + y[0]; aj = Table[a[j - i], {i, 1, j - 1}]; 
   f2j = Table[f[h*Kk, x[Kk], y1[Kk]], {Kk, 1, j - 1}];
   x[j] = par2 (aj . f2j + f[h*j, p[j], l[j]] + c[j] f0) + x[0]; 
   g2j = Table[g[F*h, x1[F], y[F]], {F, 1, j - 1}];
   y[j] = 
    par2 (aj . g2j + g[h*j, p[j], l[j]] + c[j] g0) + y[0];];, {s, 1, 
   smax}] // AbsoluteTiming

Part with $\alpha=4/5$ takes 739s. Why this code so slow compare to NDSolve? Actually we solve the system of delay integrodifferential equations, while with NDSolve we solve DDEs. Visualization

Clear[x, y];

a1 = 2/3; b1 = 4/3; c1 = 1; d1 = 1; x0 = 0.9; y0 = 1.8; 
f[t_, x_, y_] := a1*x - b1*x*y;
g[t_, x_, y_] := c1*x*y - d1*y; tau = 0.15;
sol = NDSolve[{x'[t] == f[t, x[t], y[t - tau]], 
    y'[t] == g[t, x[t - tau], y[t]], x[t /; t <= 0] == 0.9, 
    y[t /; t <= 0] == 1.8}, {x, y}, {t, 0, 15}];


pl1 = ParametricPlot[{x[t], y[t]} /. sol[[1]], {t, .0, n h}, 
   PlotRange -> All, Frame -> True, FrameLabel -> {x, y}, 
   AspectRatio -> 1/2, PlotStyle -> Green];

Show[pl1, 
 ListLinePlot[lst1, PlotStyle -> {Red, Dashed}, PlotRange -> Full], 
 ListLinePlot[lst2, PlotStyle -> {Blue, Dashed}, PlotRange -> Full]]

Figure 2

Update 1. This problem also can be solved with wavelets collocation method. In the code below we use Haar wavelets. The main limitation of this algorithm is FindRoot usage. In this example we use 32 collocation points to compute system of fractional DDE with $\alpha=q=4/5$ up to t=tmax=10. It takes about 15s on my laptop.

tmax = 10; n = 4; q = 8/10; a1 = 2/3; b1 = 4/3; c1 = 1;
  d1 = 1; x0 = 0.9; y0 = 1.8; tau = 0.15;
 M = Sum[1, {j, 0, n, 1}, {i, 0, 2^j - 1, 1}] + 1; dx = 1/M; 
 xl = Table[l dx, {l, 0, M}]; 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, M + 1}]; 
 xcolN = xcol - 1.*10^-10;
 psi1[x_] := WaveletPsi[HaarWavelet[], x];
 psi2[x_] := WaveletPhi[HaarWavelet[], x];
 psi1jk[x_, j_, k_] := psi1[j x - k];
 psi2jk[x_, j_, k_] := psi2[j x - k];
 psijk[x_, j_, k_] := Sqrt[j] (psi1jk[x, j, k] + psi2jk[x, j, k]); 
 p1i = Table[
   Flatten[Table[
     NIntegrate[psijk[t, 2^j, k], {t, 0, xcol[[i]]}], {j, 0, n, 
      1}, {k, 0, 2^j - 1, 1}]], {i, Length[xcol]}]; 
 pci = Table[
   Flatten[Table[
     NIntegrate[
      psijk[t, 2^j, k]/(xcol[[i]] - t)^q, {t, 0, xcolN[[i]]}], {j, 0, 
      n, 1}, {k, 0, 2^j - 1, 1}]], {i, Length[xcol]}];
 var = Table[
   Flatten[Table[b[i, j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]], {i, 
    2}]; ic = {.9, 1.8};
 fp = Table[
   Table[var[[j]] . p1i[[i]] + b0[j] xcol[[i]] + ic[[j]], {i, 
     Length[xcol]}], {j, 2}]; 
 fpc = Table[
   Table[var[[j]] . pci[[i]] + b0[j] xcol[[i]]^(1 - q)/(1 - q), {i, 
     Length[xcol]}], {j, 2}]; tn = (1/tmax)^q/Gamma[1 - q];
 p1tau = 
  Table[Flatten[
    Table[NIntegrate[
      psijk[t - tau/tmax, 2^j, k], {t, 0, xcol[[i]]}], {j, 0, n, 
      1}, {k, 0, 2^j - 1, 1}]], {i, Length[xcol]}]; 
 fptau = Table[
   Table[var[[j]] . p1tau[[i]] + b0[j] xcol[[i]] + ic[[j]], {i, 
     Length[xcol]}], {j, 2}];
 eq1 = Table[-fpc[[1, i]] tn + a1*fp[[1, i]] - 
     b1*fp[[1, i]]*fptau[[2, i]] == 0, {i, Length[xcol]}]; 
 eq2 = Table[-fpc[[2, i]] tn - d1*fp[[2, i]] + 
     c1*fp[[2, i]]*fptau[[1, i]] == 0, {i, Length[xcol]}];
 varM = Join[var[[1]], var[[2]], {b0[1], b0[2]}]; 
 sol1 = FindRoot[{eq1, eq2}, 
   Table[{varM[[i]], 1/10}, {i, Length[varM]}]];

Visualization along and together with Figure 2

lst = Join[{ic}, 
   Table[{fp[[1, i]], fp[[2, i]]} /. sol1, {i, Length[xcol]}]]; 
ListPlot[lst, PlotRange -> Full, PlotStyle -> {Red, PointSize[.015]}, 
 AxesLabel -> {"x", "y"}]

Figure3

Update 2. We can improve code above using exact solution for Caputo derivative from here as follows

AbsoluteTiming[
 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)))];
 tmax = 10; q = 8/10; tn = 1/tmax^q/Gamma[1 - q]; a1 = 2/3; b1 = 4/3; 
 c1 = 1; d1 = 1; x0 = 0.9; y0 = 1.8; tau = 0.15; J = 5; 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}];
 x1[t_, q_] := 
  Sum[v[i, j] pc[t, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   v1 pc1[t, q];
 x[t_] := 
  Sum[v[i, j] p[t, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   v1 p1[t] + v0;
 y1[t_, q_] := 
  Sum[u[i, j] pc[t, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   u1 pc1[t, q];
 y[t_] := 
  Sum[u[i, j] p[t, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   u1 p1[t] + 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[{-x1[t, q] tn + a1 x[t] - b1 x[t] y[t - tau/tmax] == 
      0, -y1[t, q] tn - d1 y[t] + c1 y[t] x[t - tau/tmax] == 0}, {t, 
     xcol}]]; ic = {x[0] == .9, y[0] == 1.8};
 sol = FindRoot[Join[eq[q], ic], 
   Table[{varM[[i]], 1/10}, {i, Length[varM]}]];]

It takes 4s only for 64 collocation points. As a bonus we have interpolation function for x[t], y[t]. Visualization

ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 1}, PlotRange -> All, 
 AspectRatio -> 1/2, PlotStyle -> {Red, Dashed}, Frame -> True, 
 FrameLabel -> {"x", "y"}]

Figure 4

Improve this answer
$\endgroup$
11
  • $\begingroup$ Thansk for your response. Can you please try to check the result for {t,0,15} for both alpha=1 and alpha=4/5? $\endgroup$
    – zhk
    Feb 1 at 7:46
  • $\begingroup$ @zhk Please, see update to my answer. $\endgroup$
    – Alex Trounev
    Feb 1 at 8:59
  • $\begingroup$ Would it be faster with wavelets method? $\endgroup$
    – zhk
    Feb 2 at 1:54
  • $\begingroup$ @zhk Actually we don't know, but we can try to solve this problem using the Euler wavelets collocation method and compare time with predictor-corrector algorithm. $\endgroup$
    – Alex Trounev
    Feb 2 at 3:28
  • $\begingroup$ If possible for you then please? $\endgroup$
    – zhk
    Feb 2 at 3:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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