数学的モデリング - 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 =出生率（15年前Sへ）
• 012死亡率による疾患への接触
• u上は=死亡率
• l =感染の確率=

PとSは（SにP = 15年前）だけ異なる時間帯に同じことを表しまた、すべてP値が与えられます。

S,IおよびNのフェイズポートレートを返す必要があります。私は100％確実に私のコードが私の目的であると確信しているわけではありませんが、これは私が思いついたものです。現在、コードは実行されますが、終了しません。私のコード上の提案やエラーを修正するための助けを歓迎します。

Iはまた、右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に対応します。初期状態で

t0とx0はtimes(1)=1990とS(1)、I(1)、及びN(1)に対応します。

残りのタスクは、fをMATLABが理解できるように定義することです。これらのすべてのコンポーネントが用意できたら、ode45をすぐに使用できます。

fの定義はおそらく最も難しい部分です。実装では、test_funcと、test_func（P,b,,,u）に必要な追加パラメータに対応します。

状況によっては、これらのパラメータも時間に依存することに注意することが重要です。おそらく、それらをP(t),b(t),m(t),l(t)、およびu(t)と書くことは明らかです。

この場合、組み込みの線形補間関数であるinterp1関数を見てみると役に立ちます。データポイントが与えられた場合、tは5年ごとではない場合、P(t),b(t),m(t),l(t)、およびu(t)の値を見積もることができます。

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