使用 matlab 的 ODE45 和 Runge-Kutta 方法对电晕进行优化控制

Question

如果有人可以帮助解决以下问题，我将不胜感激。我有以下 ODE：

这是我为它创建 Matlab 代码的项目

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7376536/#Equ1

人口 N(t)=S(t)+E(t)+IA(t)+IS(t)+R(t)。

模型参数名称符号值

Birth rate of the human population  b   0.00018 days−1

Natural human death rate    μ   4.563×10−5 days−1

Human life expectancy   1μ  21915 days or 60 years

Natural death rate of pathogens in the environment  μP  0.1724 days−1

Life expectancy of pathogens in the environment 1μP 5.8 days

Proportion of interaction with an infectious environment    α1  0.10

Proportion of interaction with an infectious individual α1  0.10

Rate of transmission from S to E due to contact with P  β1  0.00414

Rate of transmission from S to E due to contact with IA and/or IS   β2  0.0115

Proportion of symptomatic infectious people δ   0.7

Progression rate from E back to S due to robust immune system   ψ   0.0051

Progression rate from E to either IA or IS  ω   0.09

Death rate due to the coronavirus   σ   0.0018

Rate of recovery of the symptomatic population  γS  0.05 days−1 or 120days

Rate of recovery of the asymptomatic human population   γA  0.0714 days−1

Rate of virus spread to environment by symptomatic infectious individuals   ηS  0.1 days−1 or 110days

Rate of virus spread to environment by asymptomatic infectious individuals  ηA  0.05 days−1 or 120days

我正在尝试使用 4 阶 Runge-Kutta 方法求解耦合 ODE 系统，用于我在 Matlab 中的项目工作。

我被困在第一个循环中，无法理解在方法中使用 k1,k2,k3,k4 的后续步骤

function y = code(S0,E0,IA0,IS0,R0,P0)

% Birth rate of the human population     
b= 0.00018;

% Natural human death rate 
u = 0.00004563

% Human Life expectancy
inverse_u = 1/u;

% Natural deapth rate of pathogens in the environment
up= 0.1724;

%Life expectancy of pathogens in environment
inverse_up = 1/up;

%Proportion of interaction with an infectious environment
a1= 0.1;

% Proportion of interaction with an infectious individual
a2= 0.1;

% Rate of transmission from S to E  due to contact with P
B1= 0.00414;

% Rate of transmission from S to E  due to contact with Ι_Α and/or Ι_S
B2= 0.0115;

% δ Proportion of symptomatic infectious people
Delta = 0.7;

%ψ Progression rate from E back to S due to robust immune system
psi =  0.0051;

% ω Progression rate from E to either Ι_Αor Ι_S
Omega =  0.09;

% σ  Death rate due to the Coronavirus
Sigma = 0.0018

% νs Rate of recovery of the symptomatic population

Vs = 0.05;

% νa Rate of recovery of the symptomatic human population

Va = 0.0714;

% ηs Rate of virus spread to environment by symptomatic infectious individuals
Eta_s = 0.1;

% ηa  Rate of virus spread to environment by asymptomatic infectious individuals
Eta_a = 0.1;

test = -1;

delta = 0.001;
M = 90;
t=linspace(0,10,M+1);
h=1/M;
h2 = h/2;

S=zeros(1,M+1);
E=zeros(1,M+1);
IA=zeros(1,M+1);
IS=zeros(1,M+1);
R=zeros(1,M+1);
P=zeros(1,M+1);


S(1)=S0;
E(1)=E0;
IA(1)=IA0;
IS(1)= IS0;
R(1)=R0;
P(1)=P0;

 %????(????) = ????(????) + ????(????) +????????(????) + ????????(????) + ????(????).
N(1)=S0+E0+IA0+IS0+R0

lambda1=zeros(1,M+1);
lambda2=zeros(1,M+1);
lambda3=zeros(1,M+1);
lambda4=zeros(1,M+1);
lambda5=zeros(1,M+1);

while(test < 0)
    
    oldP = P;
    oldS = S;
    oldE = E;
    oldIA = IA;
    oldIS = IS;
    oldR = R;
    oldN = N;
    oldlambda1 = lambda1;
    oldlambda2 = lambda2;
    oldlambda3 = lambda3;
    oldlambda4 = lambda4;
     oldlambda5 =lambda5;   

    for i = 1:M
        m11 = b -( ( B1 *S(i)*P(i) ) / (1+a1*P(i) ) )  - ( ( B2*S(i) *(IA(i) + IS(i) )  ) / ( 1+a2*(IA(i) + IS(i))  ) ) + psi*E(i) - u*S(i) ;
        m12 =   ( ( B1 *S(i)*P(i) ) / (1+a1*P(i) ) ) +  ( ( B2*S(i) *(IA(i) + IS(i) )  ) / ( 1+a2*(IA(i) + IS(i))  ) )  -psi*E(i)   - u*E(i) - Omega *E(i);
        m13 = (1 - Sigma ) *Omega*E(i) - (u + Sigma  ) * IA(i) - Va *IA(i) ;
        m14 =  Sigma * Omega*E(i) - ( u + Sigma ) * IS(i) - Vs*IS(i) ;
        m15 = Vs*IS(i) +  Va *IA(i) - u*R(i) ;
        m16 =  Eta_a * IA(i) +  Eta_s * IS(i) - up*P(i);
 

Answer 1

正如我所见，方程是耦合的，因此最好以矩阵形式处理它们或使用本视频中提到的方法
https://www.youtube.com/watch?v=0LzDiScAcJI