声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2113|回复: 0

[编程技巧] 求助:如何解带时滞的微分方程组

[复制链接]
发表于 2007-4-16 11:08 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
请教高手指点,如何解如下带时滞的微分方程组:

ds(t)/dt=δ*N2(t)-β(t)*s(t)*i(t)+α*i(t)-μ*s(t-T)*j(t-T)-d1*s(t)
di(t)/dt=(t)*s(t)*i(t)-λ*i(t-T)-α*i(t)-d2*i(t)
dr(t-T)/dt=λ*i(t-T)-d3*r(t-T)
dq(t-T)/dt=μ*s(t-T)*j(t-T)-d4*q(t-T)

其中,β(t)=β0*[1-i(t)/N1(t)]^η
     j(t)=i(t)+r(t)
     N1(t)=q(t)+s(t)+i(t)+r(t)
     N=N1(t)+N2(t)

在仿真中取参数N1(0)==300000, =0.1/N1(0),μ=0.04/N1(0),η=2,i(0)=1,s(0)=N1(0)-i(0),λ=0.006,δ=0.003,α=0.001,d1=d2=d3=d4=0.0005,t<T时,r(t)=q(t)=0;

我是这样解得:
ffdu.m
function dy=ffdu(t,y,Z)
N=6*10^5;
N1=y(1)+y(2)+y(3)+y(4);
N2=N-N1;
beta0=0.1/(3*10^5);
mu=0.04/(3*10^5);
lambda=0.006;
delta=0.003;
alpha=0.001;
d1=d2=d3=d4=0.0005;
eta=2;
beta=beta0*(1-y(2)/N1)^eta;


dy=[delta*N2-beta*y(1)*y(2)+alpha*y(2)-mu*Z(1)*(Z(2)+Z(3))-d1*y(1);...
    beta*y(1)*y(2)-lambda*Z(2)-alpha*y(2)-d2*y(2);...
    lambda*Z(2)-d3*Z(3);...
    mu*Z(1)*(Z(2)+Z(3))-d4*Z(4)];


main.m
clear
tt=0:1:800;
options=ddeset('RelTol',1e-4,'AbsTol',1e-7);
sol=dde23('ffdu',[50;50;50;50],[3*10^5-1;1;0;0],tt,options);
plot(sol.x,sol.y(1,:),':k',sol.x,sol.y(2,:),'-r',sol.x,sol.y(3,:),'-.m');
axis([0 800 0 4*10^5]);


提示错误:
??? Error: File: D:\matlab\work\ffdu.m Line: 10 Column: 6
The expression to the left of the equals sign is not a valid target for an assignment.


Error in ==> dde23 at 224
f0 = feval(ddefun,t0,y0,Z0,varargin{:});


Error in ==> du at 5
sol=dde23('ffdu',[50;50;50;50],[3*10^5-1;1;0;0],tt,options);


请高手指点一下,我的解法存在什么问题,多谢!

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-13 02:25 , Processed in 0.052769 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表