Математическое моделирование - Matlab ode45-for loop

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 = коэффициент рождаемости (15 лет до S)
• = летальность
• l = инфекция шанс на контакт
• u = уровень смертности из-за болезни

P и S представляют то же самое только в разные периоды времени (P = 15 лет до S) , также указаны все значения P.

Код должен вернуть фазовый портрет S, I и N. Я определенно не на 100% уверен, что мой код прав для того, что я намереваюсь сделать, но это то, что я придумал. В настоящее время код работает, но никогда не заканчивается. Любые предложения по моему коду или помощь в исправлении ошибки приветствуются.

Я также думал о добавлении следующего внутри для петли прямо между ode45 и сюжетом, в случае необходимости:

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 в вашем implentation.

t0 и x0 в исходном состоянии соответствуют times(1)=1990 и S(1), I(1) и N(1).

Остальная задача - определить f так, как понимает MATLAB. После того, как у вас есть все эти компоненты, 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