Message Boards

WOLFRAM COMMUNITY

15751 Views

19 Replies

13 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Staff Picks Biological Sciences Control Systems Mathematics Medical Sciences Calculus Wolfram Language Events & Media

Solver for COVID-19 epidemic model with the Caputo fractional derivatives

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 4 years ago

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus First version of this code been published on https://mathematica.stackexchange.com/questions/221609/solver-for-covid-19-epidemic-model-with-the-caputo-fractional-derivatives 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 (3) presented in a paper M.A. Khan, A. Atangana, Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative, Alexandria Eng. J. (2020) with differential operator replaced with the Caputo definition for fractional derivative as follows $$\frac {d f}{dt}\rightarrow \frac {1}{\Gamma (1-q)}\int_0^t{\frac{f'(x)dx}{(t-x)^{q}}}$$ 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]; 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)^q, {x, 0, t}, Assumptions -> {t > 0, k >= 0, m > 0, q < 1}] Integrate[h1[x]/(t - x)^q, {x, 0, t}, Assumptions -> {t > 0, q < 1}] With these integrals let define functions p[x_, k_, m_] := Piecewise[{{(1 + k - mx)/m, k >= 0 && 1/m + (2k)/m - 2x < 0 && 1/m + k/m - x >= 0 && m > 0}, {(-k + mx)/m, k >= 0 && 1/m + (2k)/m - 2x >= 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_, q_] := Piecewise[{{-(t^(1 - q)/(-1 + q)), k == 0 && 1/m - 2t >= 0 && m > 0 && t > 0 && 1/m - t >= 0}, {-((m^(-1 + q)(1/(-k + mt))^(-1 + q))/(-1 + q)), k > 0 && 1/m + (2k)/m - 2t > 0 && k/m - t < 0 && m > 0 && 1/m + k/m - t > 0}, {(-t^q + 2mt^(1 + q) - mt(-(1/(2m)) + t)^q)/ (t^q(-(1/(2m)) + t)^q(m(-1 + q))), k == 0 && m > 0 && 1/m - 2t < 0 && 1/m - t >= 0}, {(1/(-1 + q))((2^(-1 + q)m^(-1 + 2q)(-(-(k/m) + t)^q - 2k(-(k/m) + t)^q + 2mt(-(k/m) + t)^q + 2k(-((1/2 + k)/m) + t)^q - 2mt(-((1/2 + k)/m) + t)^ q))/((1 + 2k - 2mt)(k - mt))^q), k > 0 && 1/m + (2k)/m - 2t == 0 && m > 0 && 1/m + k/m - t > 0}, {-((1/(-1 + q))((2^(-1 + q)m^(-1 + 2q) (-2(-((1/2 + k)/m) + t)^ q((1 + 2k - 2mt)(k - mt))^ q - 2k(-((1/2 + k)/m) + t)^q ((1 + 2k - 2mt)(k - mt))^q + 2mt(-((1/2 + k)/m) + t)^q((1 + 2k - 2mt)* (k - mt))^q + (-((1 + k)/m) + t)^q ((1 + 2k - 2mt)(k - mt))^q + 2k(-((1 + k)/m) + t)^q((1 + 2k - 2mt)(k - mt))^ q - 2mt(-((1 + k)/m) + t)^q* ((1 + 2k - 2mt)(k - mt))^ q + (-(k/m) + t)^q ((1 + 2k - 2mt)(1 + k - mt))^q + 2k(-(k/m) + t)^q((1 + 2k - 2mt)(1 + k - mt))^q - 2mt(-(k/m) + t)^q((1 + 2k - 2mt)(1 + k - mt))^ q - 2k(-((1/2 + k)/m) + t)^q* ((1 + 2k - 2mt)(1 + k - mt))^q + 2mt(-((1/2 + k)/m) + t)^q((1 + 2k - 2mt)* (1 + k - mt))^ q))/(((1 + 2k - 2mt)(k - mt))^q* ((1 + 2k - 2mt)(1 + k - mt))^q))), k > 0 && m > 0 && 1/m + (2k)/m - 2t <= 0 && 1/m + k/m - t <= 0}, {-((1/(2m(-1 + q)))((2^qm^(2q)t^q(-(1/m) + t)^q* (-(1/(2m)) + t)^q - 2^(1 + q)m^(1 + 2q)t^(1 + q)* (-(1/m) + t)^q(-(1/(2m)) + t)^q - 2^(1 + q)m^(2q)* t^q(-(1/(2m)) + t)^(2q) + 2^(1 + q)m^(1 + 2q) t^(1 + q)(-(1/(2m)) + t)^(2q) + t^q((-1 + mt)(-1 + 2mt))^q - 2mt^(1 + q)* ((-1 + mt)(-1 + 2mt))^q + 2mt(-(1/(2m)) + t)^q* ((-1 + mt)(-1 + 2mt))^q)/(t^ q(-(1/(2m)) + t)^q* ((-1 + mt)(-1 + 2mt))^q))), k == 0 && 1/m - 2t < 0 && 1/m - t < 0 && m > 0}, {(1/(-1 + q))((2^(-1 + q)m^(-1 + q)((-m^q)(-(k/m) + t)^q - 2km^q(-(k/m) + t)^q + 2m^(1 + q)t(-(k/m) + t)^q + 2km^q(-((1/2 + k)/m) + t)^q - 2m^(1 + q)t* (-((1/2 + k)/m) + t)^ q - ((1 + 2k - 2mt)(k - mt))^q (1/(-1 - 2k + 2mt))^q - 2k((1 + 2k - 2mt)(k - mt))^q* (1/(-1 - 2k + 2mt))^q + 2mt((1 + 2k - 2mt)(k - mt))^q (1/(-1 - 2k + 2mt))^q))/((1 + 2k - 2mt)(k - mt))^ q), 1/m + (2k)/m - 2t < 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)^qt + t^q - t^(1 + q))/((-1 + t)^qt^q(-1 + q)))] Now we have all functions to solve a problem with the given parametres 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 coeﬃcient)= 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; Let define variables 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}; The problem can be solved on the unit interval, hence we calculate the collocation points as follows J = 4; M = 2^J; dx = 1/(2M); A = 0; xl = Table[A + l dx, {l, 0, 2 M}]; xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}]; We can represent our solution with functions defined above as 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; All unknown variables should be joined together in one list 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}]]]; Now we are ready to solve the problem. Since we solve system of equations on the unit interval we define scaling function (120 is the length of interval time in days) tn[q_]:= (1/120)^q; eq1[t_, q_] := -tn[q]/Gamma[1 - q] Sp1[t, q] + Π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_, q_] := -tn[q]/Gamma[1 - q] Ep1[t, q] + η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_, q_] := -tn[q]/Gamma[1 - q] Ip1[ t, q] + (1 - θp) ωp Ep[t] - (τp + μp) Ip[t]; eq4[t_, q_] := -tn[q]/Gamma[1 - q] Ap1[t, q] + θp ρp Ep[ t] - (τap + μp) Ap[t]; eq5[t_, q_] := -tn[q]/Gamma[1 - q] Rp1[t, q] + τp Ip[ t] + τap Ap[t] - μp Rp[t]; eq6[t_, q_] := -tn[q]/Gamma[1 - q] Mp1[t, q] + ϱp Ip[ t] + ϖp Ap[t] - πp Mp[t]; With these equations we can calculate Figure 6 from the paper above with the next piece of code 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}}];] Visualization of solution: {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]} Attachments:* COVID19model.nb

MODERATOR NOTE: coronavirus resources & updates: https://wolfr.am/coronavirus

First version of this code been published on https://mathematica.stackexchange.com/questions/221609/solver-for-covid-19-epidemic-model-with-the-caputo-fractional-derivatives

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 (3) presented in a paper M.A. Khan, A. Atangana, Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative, Alexandria Eng. J. (2020)

with differential operator replaced with the Caputo definition for fractional derivative as follows $$\frac {d f}{dt}\rightarrow \frac {1}{\Gamma (1-q)}\int_0^t{\frac{f'(x)dx}{(t-x)^{q}}}$$ 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];

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)^q, {x, 0, t}, 
 Assumptions -> {t > 0, k >= 0, m > 0, q < 1}]

Integrate[h1[x]/(t - x)^q, {x, 0, t}, Assumptions -> {t > 0, q < 1}]

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

Now we have all functions to solve a problem with the given parametres

 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 coeﬃcient*)= 
  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;

Let define variables

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};

The problem can be solved on the unit interval, hence we calculate the collocation points as follows

 J = 4; M = 2^J; dx = 1/(2*M);  A = 0; xl = Table[A + l dx, {l, 0, 2 M}]; 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];

We can represent our solution with functions defined above as

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;

All unknown variables should be joined together in one list

 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}]]];

Now we are ready to solve the problem. Since we solve system of equations on the unit interval we define scaling function (120 is the length of interval time in days)

  tn[q_]:= (1/120)^q;
 eq1[t_, q_] := -tn[q]/Gamma[1 - q] Sp1[t, q] + Π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_, q_] := -tn[q]/Gamma[1 - q] Ep1[t, q] + η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_, q_] := -tn[q]/Gamma[1 - q] Ip1[
     t, q] + (1 - θp) ωp Ep[t] - (τp + μp) Ip[t]; 
 eq4[t_, q_] := -tn[q]/Gamma[1 - q] Ap1[t, q] + θp ρp Ep[
     t] - (τap + μp) Ap[t]; 
 eq5[t_, q_] := -tn[q]/Gamma[1 - q] Rp1[t, q] + τp Ip[
     t] + τap Ap[t] - μp Rp[t]; 
 eq6[t_, q_] := -tn[q]/Gamma[1 - q] Mp1[t, q] + ϱp Ip[
     t] + ϖp Ap[t] - πp Mp[t];

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

 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}}];]

Visualization of solution:

{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]}

POSTED BY: Alexander Trounev

19 Replies

Sort By:

EDITORIAL BOARD

EDITORIAL BOARD, WOLFRAM

Posted 4 years ago

-- you have earned *Featured Contributor Badge* Your exceptional post has been selected for our editorial column *Staff Picks* http://wolfr.am/StaffPicks and Your Profile is now distinguished by a *Featured Contributor Badge* and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: EDITORIAL BOARD

EDITORIAL BOARD

EDITORIAL BOARD, WOLFRAM

Posted 4 years ago

This post has been listed in the main resource-hub COVID-19 thread: https://wolfr.am/coronavirus in the section Computational Publications. Please feel free to add your own comment on that discussion pointing to this post ( https://community.wolfram.com/groups/-/m/t/1976589 ) so many more interested readers will become aware of your excellent work. Thank you for your contribution!

POSTED BY: EDITORIAL BOARD

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 4 years ago

Thank you for your selection. Now I should make a last step to force this code run 10-100 times faster. And then we can incorporate real data to evaluate parameters for every region.

POSTED BY: Alexander Trounev

Vedat Erturk

Posted 3 years ago

Hello Alexander,I hope you are doing fine at Corona virus times. I guess you wrote this code: https://mathematica.stackexchange.com/questions/221609/solver-for-covid-19-epidemic-model-with-the-caputo-fractional-derivatives Could you please suggest me a reference for the Haar wavelet method? Vedat

POSTED BY: Vedat Erturk

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 3 years ago

There are many papers about Haar wavelets applications. For example, method for solving ODEs with Haar wavelets described by Ü. Lepik (2009) Haar wavelet method for solving stiff differential equations, Mathematical Modelling and Analysis, 14:4, 467-481. To link to this article: http://dx.doi.org/10.3846/1392-6292.2009.14.467-481 Method for solving integral equations proposed by Ülo Lepik and Enn Tamme, Solution of nonlinear Fredholm integral equations via the Haar wavelet method, Proc. Estonian Acad. Sci. Phys. Math., 2007, 56, 1, 17–27.

POSTED BY: Alexander Trounev

Vedat Erturk

Posted 3 years ago

Thank you for your reply. I will use your code for solving my fractional system and cite one of these references.Thank you.

POSTED BY: Vedat Erturk

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 3 years ago

Actually this code is original and it is not explained in cited papers. If you gonna use this code, then, please refer to our papers as well Mohammad, M., Trounev, A., Cattani, C.: An efficient method based on framelets for solving fractional Volterra integral equations, Entropy 2020, 22(8), 824. https://doi.org/10.3390/e22080824. Mohammad, M., Trounev, A.: On the dynamical modeling of Covid-19 involving Atangana-Baleanu fractional derivative and based on Daubechies framelet simulations, Chaos, Solitons & Fractals, 140, 2020. https://doi.org/10.1016/j.chaos.2020.110171. Mohammad, M., Trounev, A.: Fractional nonlinear Volterra–Fredholm integral equations involving Atangana–Baleanu fractional derivative: framelet applications, Advances in Difference Equations, 618, (2020). https://doi.org/10.1186/s13662-020-03042-9. Mohammad, M., Trounev, A., Cattani, C.: The dynamics of COVID-19 in the UAE based on fractional derivative modeling using Riesz wavelets simulation, Adv Differ Equ 2021, 115 (2021). https://doi.org/10.1186/s13662-021-03262-7.

POSTED BY: Alexander Trounev

Vedat Erturk

Posted 3 years ago

Hello again, Alexander, Please see the attached code for my system. I am having an error. I would be appreciate if you help me correct my mistakes.Thanks. Attachments: Haar_wavelet.nb

POSTED BY: Vedat Erturk

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 3 years ago

It looks like it can not be solved with Haar wavelets. May be we need more sophisticated method since your system of equations is not trivial.

POSTED BY: Alexander Trounev

vserturk

Posted 3 years ago

Hello Alexander, I tried to solve another system via your code. Please see the attached file for the solution. I compare the solution with other two methods. I see It's really so much different. blue curve:wavelet method black curve:adams method red curve:Nuri's method order:0.7 I need your comment. Attachments: 275150805.nb sample.pdf

POSTED BY: vserturk

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 3 years ago

Please, pay attention that in our model we use the Caputo fractional derivative. Therefore we solve system of integral equations. In you model there are ordinary derivatives, so you solve system of ordinary differential equations. Also in our model the time of integration is 120 while in your model 50. You need to scale function `tn[q_]:=(1/120)^q` as `tn[q_]:=(1/50)^q`. To be close to ordinary derivative we should put, for instance, `q=99/100`.

POSTED BY: Alexander Trounev

vserturk

Posted 3 years ago

Dear Alexander Trounev, I corrected my code with your help.Thank you. Another thing is that, If it is possible, can you please send me a paper where you have used the given haar wavelets method in the sense of caputo derivative? Because I can see that you used this method in the sense of atangana-baleanu derivative. Thank you in advance. Kind regards -Vedat

POSTED BY: vserturk

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 3 years ago

The usage of Haar wavelets and Caputo fractional derivative is explained in this post, while combination of Euler wavelets and Caputo fractional derivative is described, for instance, in our paper "A novel numerical method for solving fractional diffusion-wave and nonlinear Fredholm and Volterra integral equations with zero absolute error" published in Axioms.

POSTED BY: Alexander Trounev

vserturk

Posted 3 years ago

Dear Alexander Trounev, I hope you are doing well. Regarding your code,please see the attachment for my question. I can not get the plots of S (t) versus I (t), I (t) versus R (t), and S (t) versus R (t).Could you please help me with that?.Thank you. Attachments: f3fc2398-21c0-4c...nb

POSTED BY: vserturk

vserturk

Posted 3 years ago

I tried it in a separate code since I can not do that inside the code. .Please see the attachment for my code. What can I do to get more smooth curves? Attachments: phase_plots.nb phase_plots.pdf

POSTED BY: vserturk

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 3 years ago

What problem do you try to solve? It looks like you use my code on interval `(t,0,30000}` with number of collocation points $2^8=256$, while main dynamics takes about $t=300$ for `q=99/100`. I can recommend to reduce time interval to `{t,0,300}` and plot picture like this one

POSTED BY: Alexander Trounev

vserturk

Posted 3 years ago

I considered [0,30000] instead of [0,300]. Because I can observe the stability points of the system much better. Please see the attached file for my fractional order system. Attachments: S_I_R_plots.nb image.jpg

POSTED BY: vserturk

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 3 years ago

It is nice code, but for a large interval we need to increase number of collocation points. Therefore, for 30000 we need about $2^{10}=1024$ colocation points. Then it could be problem to get solution with `FindRoot`.

POSTED BY: Alexander Trounev

vserturk

Posted 3 years ago

Hello again, Alexander, Please see the attached code for my system. I am having an error. I would be appreciate if you help me correct my mistakes.Thanks. Attachments: 2_i.nb 2_i.pdf

POSTED BY: vserturk

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback