Numerical Fractional Optimal Control of Respiratory Syncytial Virus Infection in Octave/MATLAB
Abstract
:1. Introduction
2. A Fractional-Order RSV Model
3. Numerical Resolution of the Fractional RSV Model
3.1. The fde12 Solver
>> N = 400; alpha = 0.995; >> [t,y] = model_SEIRS_fde12(N, alpha)
3.2. Fractional Forward Euler’s Method
- figure
- hold on
- plot( t,yf(1,:),t,ye(1,:),'--')
- xlabel('time')
- ylabel('S(t)')
- legend( '\it fde12','Euler');
- legend('boxoff')
- set(gca,'XTick',[0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5])
- hold off
3.3. PECE Algorithm
4. Fractional Optimal Control of RSV Transmission
Numerical Resolution of the RSV Fractional Optimal Control Problem
5. Conclusions
Author Contributions
Funding
Institutional Review Board Statement
Informed Consent Statement
Data Availability Statement
Acknowledgments
Conflicts of Interest
Appendix A. Resolution of the IVP with the fde12 Function
function [t, y] = model_SEIRS_fde12(N,alpha) % initial conditions y0=[0.426282; 0.0109566; 0.0275076; 0.535254]; % Values of parameters miu = 0.0113; niu = 36; epsilon = 91; b0 = 85; b1 = 0.167; c1 = 0.167; gama = 1.8; tfinal = 5; phi = pi/2; ft = linspace(0,tfinal,N); h = tfinal/(N-1); % Correction of values of parameters miu_ = miu^alpha; niu_ = niu^alpha; epsilon_ = epsilon^alpha; gama_ = gama^alpha; % time-dependent parameters flambda = @(t) miu.^alpha.*(1 + c1.* cos( 2.* pi.* t + phi) ); fbeta = @(t) b0.^alpha.* (1 + b1.* cos( 2.* pi.* t + phi ) ); % Differential system of equations of the model fdefun = @(t,y,ft)[flambda(t)-miu_*y(1)-fbeta(t)*y(1)*y(3)+gama_*y(4); ... fbeta(t)*y(1)*y(3)-(miu_+epsilon_)*y(2); epsilon_*y(2)-(miu_+niu_)*y(3); ... niu_*y(3)-miu_*y(4)-gama_*y(4)]; % resolution of system with solver fde12 [t,y] = fde12(alpha,fdefun,0,tfinal,y0,h,ft); end
Appendix B. Resolution of the IVP with the Forward Euler Method
function [t,y] = model_SEIRS_EULER(N,alpha) % Values of parameters miu = 0.0113; niu = 36; epsilon = 91; gama = 1.8; tfinal = 5; b0 = 85; b1 = 0.167; c1 = 0.167; phi = pi/2; % initial conditions S0 = 0.426282; E0 = 0.0109566; I0 = 0.0275076; R0 = 0.535254; % Correction of values of parameters miu_ = miu^alpha; niu_ = niu^alpha; epsilon_ = epsilon^alpha; gama_ = gama^alpha; % time-dependent parameters flambda = @(t) miu_*(1 + c1 * cos( 2 * pi * t + phi) ); fbeta = @(t) b0^alpha.* (1 + b1 * cos( 2 * pi * t + phi ) ); % Initialization of variables t = linspace(0,tfinal,N); h = tfinal/N; init = zeros(1,N); S = init; E = init; I = init; R = init; S(1) = S0; E(1) = E0; I(1) = I0; R(1) = R0; beta = fbeta(t); lambda = flambda(t); for j = 2:N aux_s = 0; aux_e = 0; aux_i = 0; aux_r = 0; for k = 1:j-1 bk = (j-k+1)^alpha-(j-k)^alpha; % Differential system of equations of the model aux_s = aux_s+bk*(lambda(k)-miu_*S(k)-beta(k)*S(k)*I(k) ... +gama_*R(k)); aux_e = aux_e+bk*(beta(k)*S(k)*I(k)-(miu_+epsilon_)*E(k)); aux_i = aux_i+bk*(epsilon_*E(k)-(miu_+niu_)*I(k)); aux_r = aux_r+bk*(niu_*I(k)-miu_*R(k)-gama_*R(k)); end S(j) = S0+h^alpha/gamma(1+alpha)*aux_s; E(j) = E0+h^alpha/gamma(1+alpha)*aux_e; I(j) = I0+h^alpha/gamma(1+alpha)*aux_i; R(j) = R0+h^alpha/gamma(1+alpha)*aux_r; end y(1,:) = S; y(2,:) = E; y(3,:) = I; y(4,:) = R; end
Appendix C. Resolution of the IVP with the PECE Method
function [t,y] = model_SEIRS_PECE(N,alpha) % Values of parameters miu = 0.0113; niu = 36; epsilon = 91; gama = 1.8; tfinal = 5; b0 = 85; b1 = 0.167; c1=0.167; phi = pi/2; % initial conditions S0 = 0.426282; E0 = 0.0109566; I0 = 0.0275076; R0 = 0.535254; % Correction of values of parameters miu_ = miu^alpha; niu_ = niu^alpha; epsilon_ = epsilon^alpha; gama_ = gama^alpha; % time-dependent parameters flambda = @(t) miu_*(1 + c1 * cos( 2 * pi * t + phi) ); fbeta = @(t) b0^alpha.* (1 + b1 * cos( 2 * pi * t + phi ) ); % Initialization of variables t = linspace(0,tfinal,N); h = tfinal/N; init = zeros(1,N); beta = fbeta(t); lambda = flambda(t); S = init; E = init; I = init; R = init; b = init; a = init; S(1) = S0; E(1) = E0; I(1) = I0; R(1) = R0; Sp = S; Ep = E; Ip = I; Rp = R; % computation of coefficients a_k and b_k for k = 1:N b(k) = k^alpha-(k-1)^alpha; a(k) = (k+1)^(alpha+1)-2*k^(alpha+1)+(k-1)^(alpha+1); end for j = 2:N % First part: prediction aux_s = 0; aux_e = 0; aux_i = 0; aux_r = 0; for k = 1:j % Differential system of equations of the model aux_s = aux_s+b(j-k+1)*(lambda(k)-miu_*S(k)... -beta(k)*S(k)*I(k)+gama_*R(k)); aux_e = aux_e+b(j-k+1)*(beta(k)*S(k)*I(k)... -(miu_+epsilon_)*E(k)); aux_i = aux_i+b(j-k+1)*(epsilon_*E(k)-(miu_+niu_)*I(k)); aux_r = aux_r+b(j-k+1)*(niu_*I(k)-miu_*R(k)-gama_*R(k)); end Sp(j) = S0+h^alpha/gamma(1+alpha)*aux_s; Ep(j) = E0+h^alpha/gamma(1+alpha)*aux_e; Ip(j) = I0+h^alpha/gamma(1+alpha)*aux_i; Rp(j) = R0+h^alpha/gamma(1+alpha)*aux_r; % Second part: correction aux_ss = lambda(j)-miu_*Sp(j)-beta(j)*Sp(j)*Ip(j)+gama_*Rp(j); aux_ee = beta(j)*Sp(j)*Ip(j)-(miu_+epsilon_)*Ep(j); aux_ii = epsilon_*Ep(j)-(miu_+niu_)*Ip(j); aux_rr = niu_*Ip(j)-miu_*Rp(j)-gama_*Rp(j); auxx = ((j-1)^(alpha+1)-(j-1-alpha)*j^alpha); aux_s0 = auxx*(lambda(1)-miu_*S(1)-beta(1)*S(1)*I(1)+gama_*R(1)); aux_e0 = auxx* (beta(1)*S(1)*I(1)-(miu_+epsilon_)*E(1)); aux_i0 = auxx*(epsilon_*E(1)-(miu_+niu_)*I(1)); aux_r0 = auxx*(niu_*I(1)-miu_*R(1)-gama_*R(1)); aux_s = 0; aux_e = 0; aux_i = 0; aux_r = 0; for k = 1:j-1 % Differential system of equations of the model aux_s = aux_s+a(j-k)*(lambda(k)-miu_*S(k)-beta(k)*S(k)*I(k)... +gama_*R(k)); aux_e = aux_e+a(j-k)*(beta(k)*S(k)*I(k)-(miu_+epsilon_)*E(k)); aux_i = aux_i+a(j-k)*(epsilon_*E(k)-(miu_+niu_)*I(k)); aux_r = aux_r+a(j-k)*(niu_*I(k)-miu_*R(k)-gama_*R(k)); end S(j) = S0+h^alpha/gamma(2+alpha)*(aux_ss+aux_s0+aux_s); E(j) = E0+h^alpha/gamma(2+alpha)*(aux_ee+aux_e0+aux_e); I(j) = I0+h^alpha/gamma(2+alpha)*(aux_ii+aux_i0+aux_i); R(j) = R0+h^alpha/gamma(2+alpha)*(aux_rr+aux_r0+aux_r); end y(1,:) = S; y(2,:) = E; y(3,:) = I; y(4,:) = R; end
Appendix D. Numerical Resolution of the Fractional Optimal Control Problem
function [t,y] = FOCP_PECE(N,alpha); % values assumed as global global tfinal miu niu epsilon gama b0 b1 c1 phi k1 k2 S0 E0 I0 R0; % Values of parameters miu = 0.0113; niu = 36; epsilon = 91; gama = 1.8; tfinal = 5; b0 = 85; b1 = 0.167; phi = pi/2; c1 = .167; % parameters of the algorithm k1 = 1; k2 = 0.001; trmax = 1.0; tol = 0.001; test = 1; % initial conditions S0 = 0.426282; E0 = 0.0109566; I0 = 0.0275076; R0 = 0.535254; % initialization of variables t = linspace(0,tfinal,N); init = zeros(1,N); S = init; E = init; I = init; R = init; p1 = init; p2 = init; p3 = init; p4 = init; Ta = init; % iterations of the numerical method while test>tol, oldS = S; oldE = E; oldI = I; oldR = R; oldp1 = p1; oldp2 = p2; oldp3 = p3; oldp4 = p4; oldTa = Ta; % forward PECE iterations [y1] = system1_control(Ta,t,N,alpha); S = y1(1,:); E = y1(2,:); I = y1(3,:); R = y1(4,:); % backward PECE iterations [y2] = system2_adjoint(S,I,Ta,t,N,alpha); p1 = y2(1,:); p2 = y2(2,:); p3 = y2(3,:); p4 = y2(4,:); % new control Ta = projection((p3-p4).*I/(2*k2),trmax); Ta = ( Ta + oldTa ) / 2; % Relative error values for convergence vector = [max(abs(S-oldS))/(max(abs(S))),... max(abs(oldE-E))/(max(abs(E))),... max(abs(oldI-I))/(max(abs(I))),... max(abs(oldR-R))/(max(abs(R))),... max(abs(oldp1-p1))/(max(abs(p1))),... max(abs(oldp2-p2))/(max(abs(p2))),... max(abs(oldp3-p3))/(max(abs(p3))),... max(abs(oldp4-p4))/(max(abs(p4))), ... max(abs(oldTa-Ta))/(max(abs(Ta)))]*100; test = max(vector); end y(1,:) = S; y(2,:) = E; y(3,:) = I; y(4,:) = R; y(5,:) = Ta; y(6,:) = p1; y(7,:) = p2; y(8,:) = p3; y(9,:) = p4; end % function II: resolution of the fractional control system function [y]= system1_control(Ta,t,N,alpha) global b0 b1 c1 phi miu gama epsilon niu tfinal S0 E0 I0 R0; % time-dependent parameters flambda = @(t) miu^alpha*(1 + c1 * cos( 2 * pi * t + phi) ); fbeta = @(t) b0^alpha.* (1 + b1 * cos( 2 * pi * t + phi ) ); % Correction of values of parameters miu_ = miu^alpha; niu_ = niu^alpha; epsilon_ = epsilon^alpha; gama_ = gama^alpha; % initialization of variables beta = fbeta(t); lambda = flambda(t); h = tfinal/N; init = zeros(1,N); S = init; E = init; I = init; R = init; a = init; b = init; S(1) = S0; E(1) = E0; I(1) = I0; R(1) = R0; Sp = init; Ep = init; Ip = init; Rp = init; % computation of coefficients a_k and b_k for k = 1:N b(k) = k^alpha-(k-1)^alpha; a(k) = (k+1)^(alpha+1)-2*k^(alpha+1)+(k-1)^(alpha+1); end for j = 2:N % First part: predict % differential equations of control system aux_s = 0; aux_e = 0; aux_i = 0; aux_r = 0; for k = 1:j aux_s = aux_s+b(j-k+1)*(lambda(k)-miu_*S(k)... -beta(k)*S(k)*I(k)+gama_*R(k)); aux_e = aux_e+b(j-k+1)*(beta(k)*S(k)*I(k)... -(miu_+epsilon_)*E(k)); aux_i = aux_i+b(j-k+1)*(epsilon_*E(k)-(miu_+niu_+Ta(k))*I(k)); aux_r = aux_r+b(j-k+1)*(niu_*I(k)-miu_*R(k)-gama_*R(k)... +Ta(k)*I(k)); end Sp(j) = S0+h^alpha/gamma(1+alpha)*aux_s; Ep(j) = E0+h^alpha/gamma(1+alpha)*aux_e; Ip(j) = I0+h^alpha/gamma(1+alpha)*aux_i; Rp(j) = R0+h^alpha/gamma(1+alpha)*aux_r; % Second part: correct aux_ss = lambda(j)-miu_*Sp(j)-beta(j)*Sp(j)*Ip(j)+gama_*Rp(j); aux_ee = beta(j)*Sp(j)*Ip(j)-(miu_+epsilon_)*Ep(j); aux_ii = epsilon_*Ep(j)-(miu_+niu_+Ta(j))*Ip(j); aux_rr = niu_*Ip(j)-miu_*Rp(j)-gama_*Rp(j)+Ta(j)*Ip(j); auxx = ((j-1)^(alpha+1)-(j-1-alpha)*j^alpha); aux_s0 = auxx*(lambda(1)-miu_*S(1)-beta(1)*S(1)*I(1)+gama_*R(1)); aux_e0 = auxx* (beta(1)*S(1)*I(1)-(miu_+epsilon_)*E(1)); aux_i0 = auxx*(epsilon_*E(1)-(miu_+niu_+Ta(1))*I(1)); aux_r0 = auxx*(niu_*I(1)-miu_*R(1)-gama_*R(1)+Ta(1)*I(1)); aux_s = 0; aux_e = 0; aux_i = 0; aux_r = 0; for k = 1:j-1 aux_s = aux_s+a(j-k)*(lambda(k)-miu_*S(k)-beta(k)*S(k)*I(k)... +gama_*R(k)); aux_e = aux_e+a(j-k)*(beta(k)*S(k)*I(k)-(miu_+epsilon_)*E(k)); aux_i = aux_i+a(j-k)*(epsilon_*E(k)-(miu_+niu_+Ta(k))*I(k)); aux_r = aux_r+a(j-k)*(niu_*I(k)-miu_*R(k)... -gama_*R(k)+Ta(k)*I(k)); end S(j) = S0+h^alpha/gamma(2+alpha)*(aux_ss+aux_s0+aux_s); E(j) = E0+h^alpha/gamma(2+alpha)*(aux_ee+aux_e0+aux_e); I(j) = I0+h^alpha/gamma(2+alpha)*(aux_ii+aux_i0+aux_i); R(j) = R0+h^alpha/gamma(2+alpha)*(aux_rr+aux_r0+aux_r); end y(1,:) = S; y(2,:) = E; y(3,:) = I; y(4,:) = R; end % function III: resolution of the fractional adjoint system function [y] = system2_adjoint(S,I,Ta,t,N,alpha) global miu gama epsilon niu tfinal k1 b0 b1 phi; % time-dependent parameter fbeta = @(t) b0^alpha.* (1 + b1 * cos( 2 * pi * t + phi ) ); % Correction of values of parameters miu_=miu^alpha; niu_=niu^alpha; epsilon_=epsilon^alpha; gama_=gama^alpha; % initialization of variables beta = fbeta(t); h = tfinal/N; init = zeros(1,N); a = init; b = init; p1 = init; p2 = init; p3 = init; p4 = init; p1p = init; p2p = init; p3p = init; p4p = init; % First part: predict S = S(end:-1:1); I = I(end:-1:1); Ta = Ta(end:-1:1); beta = beta(end:-1:1); % computation of coefficients a_k and b_k for k = 1:N b(k) = k^alpha-(k-1)^alpha; a(k) = (k+1)^(alpha+1)-2*k^(alpha+1)+(k-1)^(alpha+1); end for j = 2:N % differential equations of adjoint system aux_p1 = 0; aux_p2 = 0; aux_p3 = 0; aux_p4 = 0; for k = 1:j aux_p1 = aux_p1+b(j-k+1)*(-1)*(p1(k)*(miu_+beta(k)*I(k))- ... beta(k)*I(k)*p2(k)); aux_p2 = aux_p2+b(j-k+1)*(-1)*(p2(k)*(miu_+epsilon_)... -epsilon_*p3(k)); aux_p3 = aux_p3+b(j-k+1)*(-1)*(-k1+beta(k)*p1(k)*S(k)... -p2(k)*beta(k)*S(k)+p3(k)*(miu_+niu_+Ta(k))... -p4(k)*(niu_+Ta(k))); aux_p4 = aux_p4+b(j-k+1)*(-1)*(-gama_*p1(k)... +p4(k)*(miu_+gama_)); end p1p(j) = h^alpha/gamma(1+alpha)*aux_p1; p2p(j) = h^alpha/gamma(1+alpha)*aux_p2; p3p(j) = h^alpha/gamma(1+alpha)*aux_p3; p4p(j) = h^alpha/gamma(1+alpha)*aux_p4; % Second part: correct aux_pp1 = (-1)*(p1p(j)*(miu_+beta(j)*I(j))-beta(j)*I(j)*p2p(j)); aux_pp2 = (-1)*(p2p(j)*(miu_+epsilon_)-epsilon_*p3p(j)); aux_pp3 = (-1)*(-k1+beta(j)*p1p(j)*S(j)-p2p(j)*beta(j)*S(j)... +p3p(j)*(miu_+niu_+Ta(j))-p4p(j)*(niu_+Ta(j))); aux_pp4 = (-1)*(-gama_*p1p(j)+p4p(j)*(miu_+gama_)); auxx = (-1)*((j-1)^(alpha+1)-(j-1-alpha)*j^alpha); aux_p10 = auxx*(p1(1)*(miu_+beta(1)*I(1))-beta(1)*I(1)*p2(1)); aux_p20 = auxx*(p2(1)*(miu_+epsilon_)-epsilon_*p3(1)); aux_p30 = auxx*(-k1+beta(1)*p1(1)*S(1)-p2(1)*beta(1)*S(1)... +p3(1)*(miu_+niu_+Ta(1))-p4(1)*(niu_+Ta(1))); aux_p40 = auxx*(-gama_*p1(1)+p4(1)*(miu_+gama_)); aux_p1 = 0; aux_p2 = 0; aux_p3 = 0; aux_p4 = 0; for k = 1:j-1 aux_p1 = aux_p1+a(j-k)*(-1)*(p1(k)*(miu_+beta(k)*I(k))- ... beta(k)*I(k)*p2(k)); aux_p2 = aux_p2+a(j-k)*(-1)*( p2(k)*(miu_+epsilon_)... -epsilon_*p3(k)); aux_p3 = aux_p3+a(j-k)*(-1)*( -k1+beta(k)*p1(k)*S(k)... -p2(k)*beta(k)*S(k)+p3(k)*(miu_+niu_+Ta(k))... -p4(k)*(niu_+Ta(k))); aux_p4 = aux_p4+a(j-k)*(-1)*(-gama_*p1(k)+p4(k)*(miu_+gama_)); end p1(j) = h^alpha/gamma(2+alpha)*(aux_pp1+aux_p10+aux_p1); p2(j) = h^alpha/gamma(2+alpha)*(aux_pp2+aux_p20+aux_p2); p3(j) = h^alpha/gamma(2+alpha)*(aux_pp3+aux_p30+aux_p3); p4(j) = h^alpha/gamma(2+alpha)*(aux_pp4+aux_p40+aux_p4); end y(1,:) = p1(end:-1:1); y(2,:) = p2(end:-1:1); y(3,:) = p3(end:-1:1); y(4,:) = p4(end:-1:1); end % function IV: control projection over the set of admissible controls function [v] = projection(vect,trmax) isNeg = vect<0; vect(isNeg) = 0; isHuge = vect>trmax; vect(isHuge) = trmax; v = vect; end
References
- Leibniz, G.W. Letter from Hanover, Germany to GFA L’Hospital, September 30, 1695. Math. Schriften 1849, 2, 301–302. [Google Scholar]
- Podlubny, I. An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. In Mathematics in Science and Engineering; Academic Press, Inc.: San Diego, CA, USA, 1999; Volume 198, p. xxiv+340. [Google Scholar]
- Rezapour, S.; Etemad, S.; Agarwal, R.P.; Nonlaopon, K. On a Lyapunov-Type Inequality for Control of a Ψ-Model Thermostat and the Existence of Its Solutions. Mathematics 2022, 10, 4023. [Google Scholar] [CrossRef]
- Shah, K.; Ali, A.; Zeb, S.; Khan, A.; Alqudah, M.A.; Abdeljawad, T. Study of fractional order dynamics of nonlinear mathematical model. Alex. Eng. J. 2022, 61, 11211–11224. [Google Scholar] [CrossRef]
- Sintunavarat, W.; Turab, A. Mathematical analysis of an extended SEIR model of COVID-19 using the ABC-fractional operator. Math. Comput. Simul. 2022, 198, 65–84. [Google Scholar] [CrossRef] [PubMed]
- Area, I.; Losada, J.; Ndaïrou, F.; Nieto, J.J.; Tcheutia, D.D. Mathematical modeling of 2014 Ebola outbreak. Math. Methods Appl. Sci. 2017, 40, 6114–6122. [Google Scholar] [CrossRef]
- Carvalho, A.R.M.; Pinto, C.M.A.; Baleanu, D. HIV/HCV coinfection model: A fractional-order perspective for the effect of the HIV viral load. Adv. Differ. Equ. 2018, 1, 1–22. [Google Scholar] [CrossRef]
- Zafar, Z.U.A.; Rehan, K.; Mushtaq, M. HIV/AIDS epidemic fractional-order model. J. Differ. Equ. Appl. 2017, 23, 1298–1315. [Google Scholar] [CrossRef]
- Bandeira, T.; Carmo, M.; Lopes, H.; Gomes, C.; Martins, M.; Guzman, C.; Bangert, M.; Rodrigues, F.; Januário, G.; Tomé, T.; et al. Burden and severity of children’s hospitalizations by respiratory syncytial virus in Portugal, 2015–2018. Influenza Other Respir. Viruses 2022, 17, e13066. [Google Scholar] [CrossRef]
- Rosa, S.; Torres, D.F.M. Parameter estimation, sensitivity analysis and optimal control of a periodic epidemic model with application to HRSV in Florida. Stat. Optim. Inf. Comput. 2018, 6, 139–149. [Google Scholar] [CrossRef] [Green Version]
- Rosa, S.; Torres, D.F.M. Optimal control of a fractional order epidemic model with application to human respiratory syncytial virus infection. Chaos Solitons Fractals 2018, 117, 142–149. [Google Scholar] [CrossRef] [Green Version]
- Campos, C.; Silva, C.J.; Torres, D.F.M. Numerical optimal control of HIV transmission in Octave/MATLAB. Math. Comput. Appl. 2020, 25, 1. [Google Scholar] [CrossRef] [Green Version]
- An, G. The Crisis of Reproducibility, the Denominator Problem and the Scientific Role of Multi-scale Modeling. Bull. Math. Biol. 2018, 80, 3071–3080. [Google Scholar] [CrossRef] [PubMed] [Green Version]
- Weber, A.; Weber, M.; Milligan, P. Modeling epidemics caused by respiratory syncytial virus (RSV). Math. Biosci. 2001, 172, 95–113. [Google Scholar] [CrossRef]
- Zhang, T.; Liu, J.; Teng, Z. Existence of positive periodic solutions of an SEIR model with periodic coefficients. Appl. Math. 2012, 57, 601–616. [Google Scholar] [CrossRef] [Green Version]
- Mateus, J.P.; Rebelo, P.; Rosa, S.; Silva, C.M.; Torres, D.F.M. Optimal control of non-autonomous SEIRS models with vaccination and treatment. Discret. Contin. Dyn. Syst. Ser. S 2018, 11, 1179–1199. [Google Scholar] [CrossRef] [Green Version]
- Almeida, R. Analysis of a fractional SEIR model with treatment. Appl. Math. Lett. 2018, 84, 56–62. [Google Scholar] [CrossRef]
- Carvalho, A.R.M.; Pinto, C.M.A. Immune response in HIV epidemics for distinct transmission rates and for saturated CTL response. Math. Model. Nat. Phenom. 2019, 14, 307. [Google Scholar] [CrossRef] [Green Version]
- FLHealthCHARTS. 2015. Available online: http://www.flhealthcharts.com/charts/default.aspx (accessed on 31 January 2023).
- Florida Department of Health. Available online: http://www.floridahealth.gov/diseases-and-conditions/respiratory-syncytial-virus/ (accessed on 3 May 2017).
- Eaton, J.W.; Bateman, D.; Hauberg, S.; Wehbring, R. GNU Octave Version 7.3.0 Manual: A High-Level Interactive Language for Numerical Computations; Network Theory Limited: Boston, MA, USA, 2022. [Google Scholar]
- Garrapa, R. Predictor-Corrector PECE Method for Fractional Differential Equations. In Matlab Central File Exchange, File ID: 32918; MathWorks: Natick, MA, USA, 2023. [Google Scholar]
- Diethelm, K.; Ford, N.J.; Freed, A.D. Detailed error analysis for a fractional Adams method. Numer. Algorithms 2004, 36, 31–52. [Google Scholar] [CrossRef] [Green Version]
- Garrappa, R. On linear stability of predictor-corrector algorithms for fractional differential equations. Int. J. Comput. Math. 2010, 87, 2281–2290. [Google Scholar] [CrossRef]
- Li, C.; Zeng, F. Numerical Methods for Fractional Calculus; CRC Press: Boca Raton, FL, USA, 2015; Volume 24. [Google Scholar]
- Diethelm, K.; Ford, N.J.; Freed, A.D.; Luchko, Y. Algorithms for the fractional calculus: A selection of numerical methods. Comput. Methods Appl. Mech. Eng. 2005, 194, 743–773. [Google Scholar] [CrossRef] [Green Version]
- Diethelm, K.; Ford, N.J.; Freed, A.D. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. 2002, 29, 3–22. [Google Scholar] [CrossRef]
- Salati, A.B.; Shamsi, M.; Torres, D.F.M. Direct transcription methods based on fractional integral approximation formulas for solving nonlinear fractional optimal control problems. Commun. Nonlinear Sci. Numer. Simul. 2019, 67, 334–350. [Google Scholar] [CrossRef] [Green Version]
- Almeida, R.; Pooseh, S.; Torres, D.F.M. Computational Methods in the Fractional Calculus of Variations; Imperial College Press: London, UK, 2015; p. xii+266. [Google Scholar] [CrossRef]
- Bergounioux, M.; Bourdin, L. Pontryagin maximum principle for general Caputo fractional optimal control problems with Bolza cost and terminal constraints. ESAIM Control Optim. Calc. Var. 2020, 26, 35. [Google Scholar] [CrossRef] [Green Version]
- Lenhart, S.; Workman, J.T. Optimal Control Applied to Biological Models; Chapman & Hall/CRC Mathematical and Computational Biology Series; Chapman & Hall/CRC: Boca Raton, FL, USA, 2007; p. xii+261. [Google Scholar]
0.0113 | 36 | 1.8 | 91 | 85 | 0.167 | 0.167 |
Variable | ||||
---|---|---|---|---|
3.89849 | 0.297874 | 0.776241 | 3.8815 | |
0.221838 | 0.0196 | 0.0512966 | 0.22177 | |
0.0197738 | 0.00285878 | 0.00745357 | 0.019802 |
Variable | ||||
---|---|---|---|---|
0.191041 | 0.0162465 | 0.0415721 | 0.18758 | |
0.0115686 | 0.00106802 | 0.00269032 | 0.0113171 | |
0.00133593 | 0.000135793 | 0.000322856 | 0.00128548 |
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
Share and Cite
Rosa, S.; Torres, D.F.M. Numerical Fractional Optimal Control of Respiratory Syncytial Virus Infection in Octave/MATLAB. Mathematics 2023, 11, 1511. https://doi.org/10.3390/math11061511
Rosa S, Torres DFM. Numerical Fractional Optimal Control of Respiratory Syncytial Virus Infection in Octave/MATLAB. Mathematics. 2023; 11(6):1511. https://doi.org/10.3390/math11061511
Chicago/Turabian StyleRosa, Silvério, and Delfim F. M. Torres. 2023. "Numerical Fractional Optimal Control of Respiratory Syncytial Virus Infection in Octave/MATLAB" Mathematics 11, no. 6: 1511. https://doi.org/10.3390/math11061511
APA StyleRosa, S., & Torres, D. F. M. (2023). Numerical Fractional Optimal Control of Respiratory Syncytial Virus Infection in Octave/MATLAB. Mathematics, 11(6), 1511. https://doi.org/10.3390/math11061511