数学建模 - Matlab ode45-for 循环

Question

我需要编写代码来获取具有“历史”的数学模型的相图。我会在代码之后解释。

close all;
clear all;

times = 1990:1:2015;
hold on

b=zeros(1,26); %75-2000 per 5 years
b(1:5)=0.0358;
b(6:10)=0.0339;
b(11:15)=0.0311;
b(16:20)=0.0275;
b(21:26)=0.0249;

m=zeros(1,26); %90-2015 per 5 years
m(1:5)=0.008;
m(6:10)=0.0031;
m(11:15)=0.0137;
m(16:20)=0.0147;
m(21:26)=0.0125;

l=zeros(1,26); %90-2015 per 5 years
l(1:5)=0.015;
l(6:10)=0.031;
l(11:15)=0.026;
l(16:20)=0.015;
l(21:26)=0.014;

u=zeros(1,26); %90-2015 per 5 years
u(1:5)=0.04;
u(6:10)=0.02;
u(11:15)=0.038;
u(16:20)=0.05;
u(21:26)=0.035;

S=zeros(1,26);
I=zeros(1,26);
N=zeros(1,26);
S(1)=18442000;
I(1)=186000; %1990
N(1)=18628000;
P=zeros(1,26); %15 years before S
P(1:5)=12788000; 
P(6:10)=14731000;
P(11:15)=16968000;
P(16:20)=19696000;
P(21:26)=22893000;

for i=1:26
    [time, xy] = ode45('test_func',times,[S(i) I(i) N(i) P(i) b(i) m(i) l(i) u(i)]);
    plot(time,xy(:,1),'-g',time,xy(:,2),'-r',time,xy(:,3),'-b');
end

function rhs = test_func(t,xx)

S = xx(1);
I = xx(2);
N = xx(3);
P = xx(4);
b = xx(5);
m = xx(6);
l = xx(7);
u = xx(8);

Sdot=b*P-m*S-l*S*I; 
Idot=l*S*I-(m+u)*I;
Ndot=Sdot+Idot;

rhs = [Sdot; Idot; Ndot; P; b; m; l; u;];

end

详细列表：

• S=健康人群
• I=感染人群
• N=总人口
• b=出生率（S 之前 15 年）
• m=死亡率
• l=接触后感染几率
• u=疾病死亡率

P 和 S 代表相同的事物，只是在不同的时间段（P=S 之前的 15 年），还给出了所有 P 值。

代码需要返回S,I和N的相图。我绝对不能 100% 确定我的代码是否适合我的目标，但这就是我想出的。当前代码运行但永不结束。欢迎对我的代码提出任何建议或帮助修复错误。

如果需要，我还考虑在 ode45 和 plot 之间的 for 循环中添加以下内容：

if i<26
    xy(i+1)=S(i+1);
    xy(i+27)=I(i+1);
    xy(i+53)=N(i+1);
end

Answer 1

ode45 用于求解常微分方程。一个常微分方程问题的设置将包含一些基本组成部分。

xdot = f(t, x)

是微分方程。

x(t=t0) = x0

是初始条件。

t 是自变量，对应于您的实现中的time。

x 是因变量，对应于您的实现中的S、I 和N。

初始条件下的t0和x0对应times(1)=1990和S(1)、I(1)和N(1)。

剩下的任务是以 MATLAB 理解的方式定义 f。拥有所有这些组件后，ode45 就可以使用了。

定义f 可能是最困难的部分。在您的实现中，它对应于test_func 和test_func 所需的附加参数（P、b、m、l、u）。

请务必注意，在您的情况下，这些参数还取决于时间。也许将它们写成P(t)、b(t)、m(t)、l(t) 和u(t) 会更清楚。

在这种情况下，看看interp1 函数可能会有所帮助，它是一个内置的线性插值函数。鉴于您的数据点，MATLAB 可以估计 P(t)、b(t)、m(t)、l(t) 和 u(t) 的值，而 t 不是每 5 年一次。

function q41895153_ode45()

time_range = 1990:0.1:2015;
time_hist = 1990:5:2015;

b=zeros(1,6); %75-2000 per 5 years
b(1)=0.0358;
b(2)=0.0339;
b(3)=0.0311;
b(4)=0.0275;
b(5)=0.0249;
b(6)=0.0249;

m=zeros(1,6); %90-2015 per 5 years
m(1)=0.008;
m(2)=0.0031;
m(3)=0.0137;
m(4)=0.0147;
m(5)=0.0125;
m(6)=0.0125;

l=zeros(1,6); %90-2015 per 5 years
%{
l(1)=0.015;
l(2)=0.031;
l(3)=0.026;
l(4)=0.015;
l(5)=0.014;
l(6)=0.014;
%}

l(1)=0.001;
l(2)=0.001;
l(3)=0.001;
l(4)=0.001;
l(5)=0.001;
l(6)=0.001;

u=zeros(1,6); %90-2015 per 5 years
u(1)=0.04;
u(2)=0.02;
u(3)=0.038;
u(4)=0.05;
u(5)=0.035;
u(6)=0.035;

S0=18442;
I0=186; %1990
P=zeros(1,6); %15 years before S
P(1)=12788; 
P(2)=14731;
P(3)=16968;
P(4)=19696;
P(5)=22893;
P(6)=22893;

[time, xy] = ode45(@test_func,time_range,[S0 I0],odeset(),time_hist,P,b,m,l,u);

S = xy(:,1)
I = xy(:,2)
N = S + I

plot(time,xy);

end

function rhs = test_func(t,xx,time_hist,P,b,m,l,u)

S = xx(1);
I = xx(2);

% Interpolate to find b(t), m(t), l(t), u(t), P(t)
bt = interp1(time_hist,b,t);
mt = interp1(time_hist,m,t);
lt = interp1(time_hist,l,t);
ut = interp1(time_hist,u,t);
Pt = interp1(time_hist,P,t);

Sdot=bt*Pt-mt*S-lt*S*I; 
Idot=lt*S*I-(mt+ut)*I;

rhs = [Sdot; Idot];

end