2016-11-27 301 views
2

我试着用dde23解决延迟微分方程但似乎我所以我写的函数有错误没有正确地理解它,我不能纠正它运行,看是否输出是正确的。延迟微分方程(在MATLAB dde23)

我想解决这个系统:

我不明白为什么最后五个公式添加到程序。我必须解决使用表1参数获得的输出,诸如图在图像1这个系统。

这是我写的代码:

clear all; 
clc; 
lags=1; 
sol=dde23(@eq24,lags,@eqh,[0 80]); 
plot(sol.x,sol.y) 


function v=eqh(t) 
v=zeros(6,1); 


function v=eq24(t,s,Ia,Is,R,N) 
Alfa=0.1; 
beta1=0.09; 
beta2=0.1; 
sigma1=0.3; 
sigma2=0.4; 
mu=0.01; 
alfa=0.2; 
rho=0.4; 
r1=0.4; 
r2=0.2; 
d1=0.2; 
d2=0.15; 
k=0.1; 
p=0.8; 
tau=1; 
T=4; 
dsdt=Alfa-beta1*((s*Ia)/(1+sigma1*s))-beta2*((s*Is)/(1+sigma2*s))-mu*s+alfa*R; 
dIadt=rho.*exp(-mu*tau).*s(((beta1.*Ia)/(1+sigma1.*s))+((beta2.*Is)/(1+sigma2.*s)))-(r1+d1+mu).*Ia; 
dIsdt=(1-rho).*exp(-mu.*tau).*s(((beta1.*Ia)/(1+sigma1.*s))+((beta2.*Is)/(1+sigma2.*s)))+(1-k).*r1.*Ia-(r2+d2+mu).*Is; 
dRdt=k.*r1.*Ia+r2.*Is-mu.*R-alfa.*R; 
dNdt=Alfa-mu.*N-d1.*Ia-d2.*Is; 
+0

是什么符号't⁺'是什么意思?你的代码中的哪个部分是条件't≠nT'?你的微分方程不符合'dde23'的期望;它应该是'dydt = ddefun(t,y,Z)'的形式,其中'Z(:,n)= y(t-τₙ)'。此外,您的代码无效; 'dIadt = rho。* exp(...)。* s(((...)))'将尝试* index *'s',而不是与s(t-τ)相乘。 –

+0

看起来这是您在MATLAB中第一次尝试DDE。我认为你最好在尝试这个问题之前实现一些*更简单的DDE *,这样你就可以进行一些练习。 –

+0

您好先生Oldenhuis是否有可能分离延迟方程从其他不是然后使用dde23解决它 –

回答

0

我设法拿出一个正在运行的代码。但我仍然对你的问题的一部分感到困惑。我的代码的结果几乎与您提供的代码相同。我希望这可以帮助你,无论你在做什么。

function sol = eq242 

clf 

global lambda beta1 beta2 sigma1 sigma2 mu alpha rho r1 r2 d1 d2 k tau 

lambda=0.1;beta1=0.09;beta2=0.1;sigma1=0.3;sigma2=0.4;mu=0.01;alpha=0.2; 
rho=0.4;r1=0.4;r2=0.2;d1=0.2;d2=0.15;k=0.1; 

opts = ddeset('RelTol',1e-5,'AbsTol',1e-8); 
sol = dde23(@eq24,tau,[1, 1, 1, 1, 1],[0, 50],opts); 

figure(1) 
plot(sol(1).x,sol(1).y(1,:),sol(1).x,sol(1).y(2,:),sol(1).x,sol(1).y(3,:),sol(1).x,sol(1).y(4,:),sol(1).x,sol(1).y(5,:)) 


function dydt=eq24(t,y,Z) 

global lambda beta1 beta2 sigma1 sigma2 mu alpha rho r1 r2 d1 d2 k tau 
s = y(1); 
Ia = y(2); 
Is = y(3); 
R = y(4); 
N = y(5); 
slag1 = Z(1,1);  
ialag2 = Z(2,1); 
islag3 = Z(3,1); 

dsdt=lambda-beta1*((s*Ia)/(1+sigma1*s))-beta2*((s*Is)/(1+sigma2*s))-mu*s+alpha*R; 
dIadt=rho*exp(-mu*tau)*slag1*(((beta1*ialag2)/(1+sigma1*slag1))+((beta2*islag3)/(1+sigma2*slag1)))-(r1+d1+mu)*Ia; 
dIsdt=(1-rho)*exp(-mu*tau)*slag1*(((beta1*ialag2)/(1+sigma1*slag1))+((beta2*islag3)/(1+sigma2*slag1)))+(1-k)*r1*Ia-(r2+d2+mu)*Is; 
dRdt=k*r1*Ia+r2*Is-mu*R-alpha*R; 
dNdt=lambda-mu*N-d1*Ia-d2*Is; 
dydt= [dsdt; dIadt; dIsdt; dRdt; dNdt]; 

enter image description here