17
$\begingroup$

As it is known in biological system with memory it would be suitable to use fractional derivatives to describe evolution of the system. In a current version of Mathematica 12.1 there is no special solver for integrodifferential equations. Here we show solver with using Haar wavelets for dynamic system (13) presented in a paper M.A. Khan, A. Atangana, Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative, Alexandria Eng. J. (2020)
Figure 1

with differential operator replaced with the Caputo definition for fractional derivative as follows $$\frac {d f}{dt}\rightarrow \frac {1}{\Gamma (1-\rho)}\int_0^t{\frac{f'(x)dx}{(t-x)^{\rho}}}$$ The code below allows us to reproduce Figure 7 from the paper linked above. Let define functions

h[x_, k_, m_] := WaveletPsi[HaarWavelet[], m x - k];
h1[x_] := WaveletPhi[HaarWavelet[], x]

Let take $\rho =9/10$, and then we can calculate integrals

Integrate[h[t, k, m], {t, 0, x}, Assumptions -> {k >= 0, m > 0, x > 0}]

Integrate[h1[t], {t, 0, x}, Assumptions -> {x > 0}]

Integrate[h[x, k, m]/(t - x)^(9/10), {x, 0, t}, 
 Assumptions -> {t > 0, k >= 0, m > 0}]

Integrate[h1[x]/(t - x)^(9/10), {x, 0, t}, Assumptions -> {t > 0}]

With these integrals let define functions

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]

p1[x_] := Piecewise[{{1, x > 1}}, x]

pc[t_, k_, m_] := Piecewise[{{10*t^(1/10), k == 0 && 1/m - 2*t >= 0 && m > 0 && t > 0 && 
      1/m + (2*k)/m - 2*t >= 0 && 1/m + k/m - t >= 0}, {(10*(-k + m*t)^(1/10))/m^(1/10), 
     k > 0 && 1/m + (2*k)/m - 2*t >= 0 && k/m - t < 0 && m > 0 && 1/m + k/m - t >= 0}, 
    {(10*((-k + m*t)^(1/10) - 2^(9/10)*(-1 - 2*k + 2*m*t)^(1/10)))/m^(1/10), 
     k > 0 && 1/m + (2*k)/m - 2*t < 0 && 1/m + k/m - t >= 0 && m > 0}, 
    {(10*((-1 - k + m*t)^(1/10) + (-k + m*t)^(1/10) - 2^(9/10)*(-1 - 2*k + 2*m*t)^(1/10)))/
      m^(1/10), k > 0 && 1/m + (2*k)/m - 2*t < 0 && 1/m + k/m - t < 0 && m > 0}, 
    {(5*(2*(m*t)^(1/10) - 2^(9/10)*(-1 + 2*m*t)^(1/10) - 2^(9/10)*(-1 - 2*k + 2*m*t)^(1/10)))/
      m^(1/10), k == 0 && 1/m - 2*t < 0 && 1/m + (2*k)/m - 2*t < 0 && 1/m + k/m - t >= 0 && m > 0}, 
    {(5*(2*(m*t)^(1/10) + 2*(-1 - k + m*t)^(1/10) - 2^(9/10)*(-1 + 2*m*t)^(1/10) - 
        2^(9/10)*(-1 - 2*k + 2*m*t)^(1/10)))/m^(1/10), k == 0 && 1/m - 2*t < 0 && 
      1/m + k/m - t < 0 && m > 0}}, 0]

pc1[t_] := Piecewise[{{-10*((-1 + t)^(1/10) - t^(1/10)), t >= 1}}, 10*t^(1/10)]

Now we have all functions to solve a problem

AbsoluteTiming[ J = 4; M = 2^J; dx = 1/(2*M); 
 Np0 = 8266000; 
  μp (*Natural mortality rate*)= 
  1/(76.79 365); Πp (*Birth rate*)= μp Np0 ; ηp \
(*Contact rate*)= 0.05; ψ (*Transmissibility multiple*) = 
  0.02; ηw (*Disease transmission coefficient*)= 
  0.000001231; θp (*The proportion of asymptomatic \
infection*)= 0.1243; ωp (*Incubation period*)= 
  0.00047876;  ρp (*Incubation period*)= 
  0.005;  τp (*Removal or recovery rate of Ip*)= 
  0.09871;  τap (*Removal or recovery rate of Ap *)= 
  0.854302; ϱp (*Contribution of the virus to M by Ip*)= 
  0.000398; ϖp (*Contribution of the virus to M by Ap*) = 
  0.001; πp(*Removing rate of virus from M*) = 0.01;

 var1 = {Sp1, Ep1, Ip1, Ap1, Rp1, Mp1}; 
 var = {Sp, Ep, Ip, Ap, Rp, Mp}; aco = {aS, aE, aI, aA, aR, aM}; 
 aco1 = {aS1, aE1, aI1, aA1, aR1, aM1}; 
 aco0 = {aS0, aE0, aI0, aA0, aR0, aM0};
 A = 0; xl = Table[A + l dx, {l, 0, 2 M}]; 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
 Sp1[x_] := 
  Sum[aS[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 pc1[x]; 
 Sp[x_] := 
  Sum[aS[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 p1[x] + aS0; 
 Ep1[x_] := 
  Sum[aE[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 pc1[x]; 
 Ep[x_] := 
  Sum[aE[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 p1[x] + aE0; 
 Ip1[x_] := 
  Sum[aI[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 pc1[x]; 
 Ip[x_] := 
  Sum[aI[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 p1[x] + aI0; 
 Ap1[x_] := 
  Sum[aA[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 pc1[x]; 
 Ap[x_] := 
  Sum[aA[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 p1[x] + aA0; 
 Rp1[x_] := 
  Sum[aR[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 pc1[x]; 
 Rp[x_] := 
  Sum[aR[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 p1[x] + aR0; 
 Mp1[x_] := 
  Sum[aM[i, j] pc[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 pc1[x]; 
 Mp[x_] := 
  Sum[aM[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 p1[x] + aM0;

 varM = Join[aco0, aco1, 
   Flatten[Table[{aS[i, j], aE[i, j], aI[i, j], aA[i, j], aR[i, j], 
      aM[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]];
 ρ = 9/10; tn = (1/120);
 eq1[t_] := -tn/Gamma[1 - ρ] Sp1[t] + Πp/
    Np0 - μp Sp[t] - ηp Sp[
     t] (Ip[t] + ψ Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) - Np0 ηw Sp[t] Mp[t]; 
 eq2[t_] := -tn/Gamma[1 - ρ] Ep1[t] + ηp  Sp[
     t] (Ip[t] + ψ Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) + 
   Np0 ηw Sp[t] Mp[t] - (1 - θp) ωp Ep[
     t] - θp ρp Ep[t] - μp Ep[t];
 eq3[t_] := -tn/Gamma[1 - ρ] Ip1[
     t] + (1 - θp) ωp Ep[t] - (τp + μp) Ip[t]; 
 eq4[t_] := -tn/Gamma[1 - ρ] Ap1[t] + θp ρp Ep[
     t] - (τap + μp) Ap[t]; 
 eq5[t_] := -tn/Gamma[1 - ρ] Rp1[t] + τp Ip[
     t] + τap Ap[t] - μp Rp[t]; 
 eq6[t_] := -tn/Gamma[1 - ρ] Mp1[t] + ϱp Ip[
     t] + ϖp Ap[t] - πp Mp[t];


 eq = Flatten[
   ParallelTable[{eq1[t] == 0, eq2[t] == 0, eq3[t] == 0, eq4[t] == 0, 
     eq5[t] == 0, eq6[t] == 0}, {t, xcol}]];
 Do[icv[i] = {Sp[0] == 8065518/Np0/8 i, Ep[0] == 200000/Np0, 
    Ip[0] == 282/Np0, Ap[0] == 200/Np0, Rp[0] == 0, 
    Mp[0] == 50000/Np0};
  eqM = Join[eq, icv[i]];
  solv[i] = 
   FindRoot[eqM, Table[{varM[[j]], .1}, {j, Length[varM]}], 
    MaxIterations -> 1000];
  lstSv[i] = 
   Table[{x 120 , Np0 Evaluate[Sp[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstEv[i] = 
   Table[{x 120, Np0 Evaluate[Ep[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstIv[i] = 
   Table[{x 120, Np0 Evaluate[Ip[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstAv[i] = 
   Table[{x 120, Np0 Evaluate[Ap[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstRv[i] = 
   Table[{x 120, Np0 Evaluate[Rp[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstMv[i] = 
   Table[{x 120, Np0 Evaluate[Mp[x] /. solv[i]]}, {x, 0, 
     1, .01}];, {i, 1, 8}]]

Finally we visualize solution

{ListLinePlot[Table[lstSv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(S\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstEv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(E\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstIv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(I\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstAv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(A\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstRv[i], {i, 1, 8}], Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(R\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[Table[lstMv[i], {i, 1, 8}], Frame -> True, 
  FrameLabel -> {"t, days", "M"}, 
     PlotRange -> All, PlotLegends -> Automatic]}

Figure 2

The question is about how to add $\rho$ in this code as parameter?

Update 1. The straight forward solution of this problem is simply to include $\rho$ in the pc, pc1functions definitions as follows (here $\rho$ is replaced by q) :

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

With these functions we can calculate Figure 6 from the paper above with the next piece of code

AbsoluteTiming[J = 4; M = 2^J; dx = 1/(2*M); 
 Np0 = 8266000; 
  \[Mu]p (*Natural mortality rate*)= 
  1/(76.79 365); \[CapitalPi]p (*Birth rate*)= \[Mu]p Np0 ; \[Eta]p \
(*Contact rate*)= 0.05; \[Psi] (*Transmissibility multiple*) = 
  0.02; \[Eta]w (*Disease transmission coefficient*)= 
  0.000001231; \[Theta]p (*The proportion of asymptomatic \
infection*)= 0.1243; \[Omega]p (*Incubation period*)= 
  0.00047876;  \[Rho]p (*Incubation period*)= 
  0.005;  \[Tau]p (*Removal or recovery rate of Ip*)= 
  0.09871;  \[Tau]ap (*Removal or recovery rate of Ap *)= 
  0.854302; \[CurlyRho]p (*Contribution of the virus to M by Ip*)= 
  0.000398; \[CurlyPi]p (*Contribution of the virus to M by Ap*) = 
  0.001; \[Pi]p(*Removing rate of virus from M*) = 0.01;

 var1 = {Sp1, Ep1, Ip1, Ap1, Rp1, Mp1}; 
 var = {Sp, Ep, Ip, Ap, Rp, Mp}; aco = {aS, aE, aI, aA, aR, aM}; 
 aco1 = {aS1, aE1, aI1, aA1, aR1, aM1}; 
 aco0 = {aS0, aE0, aI0, aA0, aR0, aM0};
 A = 0; xl = Table[A + l dx, {l, 0, 2 M}]; 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
 Sp1[x_, q_] := 
  Sum[aS[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 pc1[x, q]; 
 Sp[x_] := 
  Sum[aS[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 p1[x] + aS0; 
 Ep1[x_, q_] := 
  Sum[aE[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 pc1[x, q]; 
 Ep[x_] := 
  Sum[aE[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 p1[x] + aE0; 
 Ip1[x_, q_] := 
  Sum[aI[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 pc1[x, q]; 
 Ip[x_] := 
  Sum[aI[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 p1[x] + aI0; 
 Ap1[x_, q_] := 
  Sum[aA[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 pc1[x, q]; 
 Ap[x_] := 
  Sum[aA[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 p1[x] + aA0; 
 Rp1[x_, q_] := 
  Sum[aR[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 pc1[x, q]; 
 Rp[x_] := 
  Sum[aR[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 p1[x] + aR0; 
 Mp1[x_, q_] := 
  Sum[aM[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 pc1[x, q]; 
 Mp[x_] := 
  Sum[aM[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 p1[x] + aM0;

 varM = Join[aco0, aco1, 
   Flatten[Table[{aS[i, j], aE[i, j], aI[i, j], aA[i, j], aR[i, j], 
      aM[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]];
 tn[q_] := (1/120)^q;
 eq1[t_, q_] := -tn[q]/Gamma[1 - q] Sp1[t, q] + \[CapitalPi]p/
    Np0 - \[Mu]p Sp[t] - \[Eta]p Sp[
     t] (Ip[t] + \[Psi] Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) - Np0 \[Eta]w Sp[t] Mp[t]; 
 eq2[t_, q_] := -tn[q]/Gamma[1 - q] Ep1[t, q] + \[Eta]p  Sp[
     t] (Ip[t] + \[Psi] Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) + 
   Np0 \[Eta]w Sp[t] Mp[t] - (1 - \[Theta]p) \[Omega]p Ep[
     t] - \[Theta]p \[Rho]p Ep[t] - \[Mu]p Ep[t];
 eq3[t_, q_] := -tn[q]/Gamma[1 - q] Ip1[t, 
     q] + (1 - \[Theta]p) \[Omega]p Ep[t] - (\[Tau]p + \[Mu]p) Ip[t]; 
 eq4[t_, q_] := -tn[q]/Gamma[1 - q] Ap1[t, q] + \[Theta]p \[Rho]p Ep[
     t] - (\[Tau]ap + \[Mu]p) Ap[t]; 
 eq5[t_, q_] := -tn[q]/Gamma[1 - q] Rp1[t, q] + \[Tau]p Ip[
     t] + \[Tau]ap Ap[t] - \[Mu]p Rp[t]; 
 eq6[t_, q_] := -tn[q]/Gamma[1 - q] Mp1[t, q] + \[CurlyRho]p Ip[
     t] + \[CurlyPi]p Ap[t] - \[Pi]p Mp[t];


 eq[q_] := 
  Flatten[ParallelTable[{eq1[t, q] == 0, eq2[t, q] == 0, 
     eq3[t, q] == 0, eq4[t, q] == 0, eq5[t, q] == 0, 
     eq6[t, q] == 0}, {t, xcol}]];
 Do[icv[i] = {Sp[0] == 8065518/Np0, Ep[0] == 200000/Np0, 
    Ip[0] == 282/Np0, Ap[0] == 200/Np0, Rp[0] == 0, 
    Mp[0] == 50000/Np0};
  eqM[i] = Join[eq[i], icv[i]];
  solv[i] = 
   FindRoot[eqM[i], Table[{varM[[j]], .1}, {j, Length[varM]}], 
    MaxIterations -> 1000];
  lstSv[i] = 
   Table[{x 120 , Np0 Evaluate[Sp[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstEv[i] = 
   Table[{x 120, Np0 Evaluate[Ep[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstIv[i] = 
   Table[{x 120, Np0 Evaluate[Ip[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstAv[i] = 
   Table[{x 120, Np0 Evaluate[Ap[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstRv[i] = 
   Table[{x 120, Np0 Evaluate[Rp[x] /. solv[i]]}, {x, 0, 1, .01}]; 
  lstMv[i] = 
   Table[{x 120, Np0 Evaluate[Mp[x] /. solv[i]]}, {x, 0, 
     1, .01}];, {i, {99/100, 9/10, 8/10, 7/10, 6/10}}];] 

We can check that it is run 4-5 times longer than code with a fixed $\rho$. Visualization:

{ListLinePlot[Table[lstSv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(S\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstEv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(E\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstIv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(I\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstAv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(A\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstRv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, 
     FrameLabel -> {"t, days", "\!\(\*SubscriptBox[\(R\), \(p\)]\)"}, 
  PlotRange -> All], 
   ListLinePlot[
  Table[lstMv[i], {i, {99/100, 9/10, 8/10, 7/10, 6/10}}], 
  Frame -> True, FrameLabel -> {"t, days", "M"}, 
     PlotRange -> All, PlotLegends -> Automatic]}  

Figure 2

Update 2. We can reduce time by 3-4 times simply replace where it is possible function definition f[x_,...]:=... with f=Compile[{{x,_Real},{...}},...]. So in the last code we have to replace first part as follows

h = Compile[{{x, _Real}, {k, _Integer}, {m, _Integer}}, 
   WaveletPsi[HaarWavelet[], m x - k]];

p = Compile[{{x, _Real}, {k, _Integer}, {m, _Integer}}, 
   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 = Compile[{{x, _Real}}, WaveletPhi[HaarWavelet[], x]];

p1 = Compile[{{x, _Real}}, Piecewise[{{1, x > 1}}, x]];

pc = Compile[{{t, _Real}, {k, _Integer}, {m, _Integer}, {q, _Real}}, 
   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 = Compile[{{t, _Real}, {q, _Real}}, 
  Piecewise[{{-(t^(1 - q)/(-1 + q)), t <= 1}}, 
     -(((-1 + t)^q*t + t^q - t^(1 + q))/((-1 + t)^q*
        t^q*(-1 + q)))]]; tn = Compile[{{q, _Real}}, (1/120)^q];
$\endgroup$
  • 2
    $\begingroup$ "it would be suitable to use fractional derivatives to describe evolution of the system" - just out of curiosity, do you have a good reference that demonstrates this? $\endgroup$ – Chris K May 10 at 19:12
  • 2
    $\begingroup$ @Szabolcs This code is quite good for particular $\rho$ (it takes few second to calculate solution). To make it faster I have introduced functions like pc[t,k,m], pc1[t,k,m] calculated with a given $\rho$. But when I try to use it in general form as a functions of $\rho$, it takes hundred times longer. $\endgroup$ – Alex Trounev May 10 at 20:09
  • 3
    $\begingroup$ @ChrisK I have few papers about application of fractional derivatives to biological systems including this one arxiv.org/abs/2005.01820 $\endgroup$ – Alex Trounev May 10 at 20:28
  • 4
    $\begingroup$ Note that there is now ResourceFunction["FractionalIntegrate"] in the Wolfram Function Repository. $\endgroup$ – Daniel Lichtblau May 11 at 15:57
  • 1
    $\begingroup$ @Daniel, I think Alex would use ResourceFunction["FractionalD"] instead for this, since he wants the Caputo derivative. $\endgroup$ – J. M.'s technical difficulties May 11 at 16:54
7
$\begingroup$

WaveletPsi is not compilable, so I changed some of your compile definitions. Changing the parallel evaluation a little bit your code runs now in about a second on a simple 4-core machine. Using the experimental FunctionCompile does not help at all, since both the Piecewise and the HaarWavelet function generate compile errors. Looks like those functions are not among the "approximately 2000 functions which cover 31 functionality areas" [Compiler paper].

h = Function[{x, k, m}, WaveletPsi[HaarWavelet[], m x - k]]; 

p = Compile[{{x, _Real}, {k, _Integer}, {m, _Integer}}, 
   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 = Function[{x}, WaveletPhi[HaarWavelet[], x]];

p1 = Function[x, Piecewise[{{1, x > 1}}, x]];

pc = Compile[{{t, _Real}, {k, _Integer}, {m, _Integer}, {q, _Real}}, 
   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 = Compile[{{t, _Real}, {q, _Real}}, 
   Piecewise[{{-(t^(1 - q)/(-1 + q)), 
      t <= 1}}, -(((-1 + t)^q*t + t^q - t^(1 + q))/((-1 + t)^q*
         t^q*(-1 + q)))]];

tn = Function[{q}, (1/120)^q];


(*Now we have all functions to solve a problem with the given \
parameres*)

LaunchKernels[] (* Launch parallel kernels before doing \
AbsoluteTiming *)

AbsoluteTiming[
 J = 4; M = 2^J; dx = 1/(2*M);
 Np0 = 8266000;
 \[Mu]p (*Natural mortality rate*)= 
  1/(76.79 365); \[CapitalPi]p (*Birth rate*)= \[Mu]p Np0; \[Eta]p \
(*Contact rate*)= 0.05; \[Psi] (*Transmissibility multiple*)= 
  0.02; \[Eta]w (*Disease transmission coefficient*)= 
  0.000001231; \[Theta]p (*The proportion of asymptomatic infection*)=
   0.1243; \[Omega]p (*Incubation period*)= 
  0.00047876; \[Rho]p (*Incubation period*)= 
  0.005; \[Tau]p (*Removal or recovery rate of Ip*)= 
  0.09871; \[Tau]ap (*Removal or recovery rate of Ap*)= 
  0.854302; \[CurlyRho]p (*Contribution of the virus to M by Ip*)= 
  0.000398; \[CurlyPi]p (*Contribution of the virus to M by Ap*)= 
  0.001; \[Pi]p(*Removing rate of virus from M*)= 0.01;
 var1 = {Sp1, Ep1, Ip1, Ap1, Rp1, Mp1};
 var = {Sp, Ep, Ip, Ap, Rp, Mp}; aco = {aS, aE, aI, aA, aR, aM};
 aco1 = {aS1, aE1, aI1, aA1, aR1, aM1};
 aco0 = {aS0, aE0, aI0, aA0, aR0, aM0};
 A = 0; xl = Table[A + l dx, {l, 0, 2 M}];
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
 Sp1[x_, q_] := 
  Sum[aS[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 pc1[x, q];
 Sp[x_] := 
  Sum[aS[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aS1 p1[x] + aS0;
 Ep1[x_, q_] := 
  Sum[aE[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 pc1[x, q];
 Ep[x_] := 
  Sum[aE[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aE1 p1[x] + aE0;
 Ip1[x_, q_] := 
  Sum[aI[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 pc1[x, q];
 Ip[x_] := 
  Sum[aI[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aI1 p1[x] + aI0;
 Ap1[x_, q_] := 
  Sum[aA[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 pc1[x, q];
 Ap[x_] := 
  Sum[aA[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aA1 p1[x] + aA0;
 Rp1[x_, q_] := 
  Sum[aR[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 pc1[x, q];
 Rp[x_] := 
  Sum[aR[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aR1 p1[x] + aR0;
 Mp1[x_, q_] := 
  Sum[aM[i, j] pc[x, i, 2^j, q], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 pc1[x, q];
 Mp[x_] := 
  Sum[aM[i, j] p[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
   aM1 p1[x] + aM0;
 varM = Join[aco0, aco1, 
   Flatten[Table[{aS[i, j], aE[i, j], aI[i, j], aA[i, j], aR[i, j], 
      aM[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]];

 eq1[t_, q_] := -tn[q]/Gamma[1 - q] Sp1[t, q] + \[CapitalPi]p/
    Np0 - \[Mu]p Sp[t] - \[Eta]p Sp[
     t] (Ip[t] + \[Psi] Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) - Np0 \[Eta]w Sp[t] Mp[t];
 eq2[t_, q_] := -tn[q]/Gamma[1 - q] Ep1[t, q] + \[Eta]p Sp[
     t] (Ip[t] + \[Psi] Ap[t])/(Sp[t] + Ep[t] + Ip[t] + Ap[t] + 
       Rp[t]) + 
   Np0 \[Eta]w Sp[t] Mp[t] - (1 - \[Theta]p) \[Omega]p Ep[
     t] - \[Theta]p \[Rho]p Ep[t] - \[Mu]p Ep[t];
 eq3[t_, q_] := -tn[q]/Gamma[1 - q] Ip1[t, 
     q] + (1 - \[Theta]p) \[Omega]p Ep[t] - (\[Tau]p + \[Mu]p) Ip[t];
 eq4[t_, q_] := -tn[q]/Gamma[1 - q] Ap1[t, q] + \[Theta]p \[Rho]p Ep[
     t] - (\[Tau]ap + \[Mu]p) Ap[t];
 eq5[t_, q_] := -tn[q]/Gamma[1 - q] Rp1[t, q] + \[Tau]p Ip[
     t] + \[Tau]ap Ap[t] - \[Mu]p Rp[t];
 eq6[t_, q_] := -tn[q]/Gamma[1 - q] Mp1[t, q] + \[CurlyRho]p Ip[
     t] + \[CurlyPi]p Ap[t] - \[Pi]p Mp[t];
 eq[q_] := 
  Flatten[Table[{eq1[t, q] == 0, eq2[t, q] == 0, eq3[t, q] == 0, 
     eq4[t, q] == 0, eq5[t, q] == 0, eq6[t, q] == 0}, {t, xcol}]];
 ParallelDo[ (* It is more efficient to parallelize this Do loop *)
  icv[i] = {Sp[0] == 8065518/Np0, Ep[0] == 200000/Np0, 
    Ip[0] == 282/Np0, Ap[0] == 200/Np0, Rp[0] == 0, 
    Mp[0] == 50000/Np0};
  eqM[i] = Join[eq[i], icv[i]];
  solv[i] = 
   FindRoot[eqM[i], Table[{varM[[j]], .1}, {j, Length[varM]}], 
    MaxIterations -> 1000];
  lstSv[i] = 
   Table[{x 120, Np0 Evaluate[Sp[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstEv[i] = 
   Table[{x 120, Np0 Evaluate[Ep[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstIv[i] = 
   Table[{x 120, Np0 Evaluate[Ip[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstAv[i] = 
   Table[{x 120, Np0 Evaluate[Ap[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstRv[i] = 
   Table[{x 120, Np0 Evaluate[Rp[x] /. solv[i]]}, {x, 0, 1, .01}];
  lstMv[i] = 
   Table[{x 120, Np0 Evaluate[Mp[x] /. solv[i]]}, {x, 0, 1, .01}];
  , {i, {99/100, 9/10, 8/10, 7/10, 6/10}}];
 (* Collect the definitions from the parallel kernels to the main \
one, e.g. by: (in principle this should work also by SharedFunction, \
but tat does not work well )*)
 Table[With[{lst = lst}, 
   DownValues[lst] = 
    Flatten@ParallelEvaluate[DownValues[lst]]], {lst, {lstSv, lstEv, 
    lstIv, lstAv, lstRv, lstMv}}];
 ]
| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you! Actually I didn't see difference with my code if a piece of code in AbsoluteTiming[] using many times for parameters evaluation. In this case ParallelTable[] running even faster. $\endgroup$ – Alex Trounev May 18 at 10:27
  • $\begingroup$ It can be very tricky to optimize parallel evaluations in Mathematica. It also depends on the nature of the problem and also on hardware. $\endgroup$ – Rolf Mertig May 18 at 11:07
  • $\begingroup$ If memory serves, WaveletPsi[] and WaveletPhi[] are supposed to produce Piecewise[] functions, so if need be, someone could convert the function (by hand?) into something with Which[] (or re-express in terms of UnitStep[]) so that it could be compileable. $\endgroup$ – J. M.'s technical difficulties May 18 at 11:58
  • 1
    $\begingroup$ @Alex, that must have been a recent thing (I'm still on version 11); if it is working for you, then my concerns are moot. But just to check: you can use CompilePrint[] on your functions and check for MainEvaluate calls in the output. $\endgroup$ – J. M.'s technical difficulties May 18 at 17:01
  • 1
    $\begingroup$ Ah, so the compilation did not succeed in that case, @Alex. $\endgroup$ – J. M.'s technical difficulties May 18 at 17:22

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.