24小时热门版块排行榜    

查看: 1294  |  回复: 7

xwndf250

银虫 (小有名气)

[求助] matlab处理常微分方程作图问题

流行病模型 sir
dI/dt=a*S*I-b*I
dS/dt=-a*S*I
dR/dt=b*I
S+I+R=1且0 想要做一个横坐标t=[0,50], 纵坐标dI/dt、dS/dt、dR/dt的3条曲线,求具体程序,本人写的程序
m文件:function y=SIR(t,x)
a=0.2;b=0.1;
y=[a*x(1)*x(2)-b*x(1);
-a*x(1)*x(2);
b*x(3)];
end
命令:应该怎样写,才是t与dI/dt、dS/dt、dR/dt的图像(就是S、I、R的变化率随时间的变化),注意,不是t与S、I、R的图像!最好用matlab语言。
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

csgt0

荣誉版主 (著名写手)

彩色挂图

【答案】应助回帖

★ ★
感谢参与,应助指数 +1
dbb627: 金币+2, 感谢应助 2013-02-26 14:07:58
CODE:
function xwn
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
intvalue=[0.1 0.4 0.5];       %t=0时的初值
[T,Y] = ode45(@rigid,[0 50],intvalue,options);
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
title('T-Y图')
for i=1:length(T)     
vdy(:,i)=rigid(T(i),Y(i,:));
end
figure
plot(T,vdy(1,:),'-',T,vdy(2,:),'-.',T,vdy(3,:),'.')
title('T-dY图')
end

function dy = rigid(t,y)
dy = zeros(3,1);    % a column vector
a=0.2;
b=0.1;
dy(1) = a*y(2)*y(1)-b * y(1);
dy(2) = -a*y(2) * y(1);
dy(3) = b * y(1);
end

showmethemoney
2楼2013-02-26 11:06:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xwndf250

银虫 (小有名气)

引用回帖:
2楼: Originally posted by csgt0 at 2013-02-26 11:06:11
function xwn
options = odeset('RelTol',1e-4,'AbsTol',);
intvalue=;       %t=0时的初值
= ode45(@rigid,,intvalue,options);
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
title('T-Y图')
for i=1: ...

function dy = rigid(t,y)后面的我看懂了,但是前面function xwn没有看懂,我记得求导貌似用diff函数的,能不能详细讲解下odeset()函数,intvalue作用,“for i=1:length(T) 和vdy(:,i)=rigid(T(i),Y(i,);”的意思。我是小白,求指教。好的话我多给金币。
3楼2013-02-26 21:06:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

csgt0

荣誉版主 (著名写手)

彩色挂图

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
xwndf250: 金币+15, ★★★很有帮助, 非常感谢您的帮助,给您金币,顺便问下,t0初值怎么确定的?也就是0.1,0.4,0.5是怎么确定的? 2013-02-27 10:00:25
fegg7502: 金币+1, 应助指数+1, 鼓励交流 2013-04-02 09:30:14
odeset()函数设置数值解的计算精度,其实一般不用也可以
intvalue是初值啊,t=0时的初值
“for i=1:length(T) 和vdy(:,i)=rigid(T(i),Y(i,);”的为了循环计算每个t下的3个导数值
showmethemoney
4楼2013-02-27 09:54:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

csgt0

荣誉版主 (著名写手)

彩色挂图

【答案】应助回帖


fegg7502: 金币+1, 应助指数+1, 鼓励交流 2013-04-02 09:30:24
引用回帖:
3楼: Originally posted by xwndf250 at 2013-02-26 21:06:43
function dy = rigid(t,y)后面的我看懂了,但是前面function xwn没有看懂,我记得求导貌似用diff函数的,能不能详细讲解下odeset()函数,intvalue作用,“for i=1:length(T) 和vdy(:,i)=rigid(T(i),Y(i,);”的 ...

就是t=0是的ISR啊,具体多少得看你的实际情况,我只是随意写的个数。
showmethemoney
5楼2013-02-27 10:37:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

xwndf250

银虫 (小有名气)

引用回帖:
5楼: Originally posted by csgt0 at 2013-02-27 10:37:52
就是t=0是的ISR啊,具体多少得看你的实际情况,我只是随意写的个数。...

好的 谢谢你
6楼2013-02-27 11:09:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

657801288

铁虫 (初入文坛)

引用回帖:
2楼: Originally posted by csgt0 at 2013-02-26 11:06:11
function xwn
options = odeset('RelTol',1e-4,'AbsTol',);
intvalue=;       %t=0时的初值
= ode45(@rigid,,intvalue,options);
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
title('T-Y图')
for i=1: ...

你好,两种方法做出的图怎么不一样呢?求指教。
7楼2013-03-30 08:49:39
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

csgt0

荣誉版主 (著名写手)

彩色挂图

【答案】应助回帖


fegg7502: 金币+1, 应助指数+1, 鼓励交流 2013-04-02 09:30:38
不是两种方法啊,前面画的是T-Y图,后面是你要的T-dY图。
showmethemoney
8楼2013-04-01 10:48:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 xwndf250 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 336求调剂 +4 收到VS 2026-03-20 4/200 2026-03-23 19:02 by macy2011
[考研] 上海电力大学材料防护与新材料重点实验室招收调剂研究生(材料、化学、电化学,环境) +3 我爱学电池 2026-03-23 3/150 2026-03-23 17:16 by AZMK
[考研] 材料与化工考研调剂 +4 孅華 2026-03-22 4/200 2026-03-23 16:13 by 一休哥FU
[考研] 307求调剂 +3 余意卿 2026-03-21 3/150 2026-03-23 10:32 by Iveryant
[考研] 08工科 320总分 求调剂 +11 梨花珞晚风 2026-03-17 11/550 2026-03-22 17:42 by luoyongfeng
[考研] 求调剂 +6 十三加油 2026-03-21 6/300 2026-03-22 17:00 by i_cooler
[考研] 生物学调剂 +5 Surekei 2026-03-21 5/250 2026-03-22 14:39 by tcx007
[考研] 0805 316求调剂 +3 大雪深藏 2026-03-18 3/150 2026-03-21 18:55 by 学员8dgXkO
[考研] 296求调剂 +4 www_q 2026-03-20 4/200 2026-03-21 17:26 by 学员8dgXkO
[考研] 一志愿重庆大学085700资源与环境总分308求调剂 +7 墨墨漠 2026-03-20 7/350 2026-03-21 16:36 by barlinike
[考研] 0805材料320求调剂 +3 深海物语 2026-03-20 3/150 2026-03-21 15:46 by 无际的草原
[考研] 268求调剂 +9 简单点0 2026-03-17 9/450 2026-03-21 15:37 by lature00
[考研] 304求调剂 +6 曼殊2266 2026-03-18 6/300 2026-03-21 00:32 by JourneyLucky
[考研] 304求调剂 +7 司空. 2026-03-18 7/350 2026-03-20 23:08 by JourneyLucky
[考研] 考研调剂求学校推荐 +3 伯乐29 2026-03-18 5/250 2026-03-20 22:59 by JourneyLucky
[考研] 材料学硕297已过四六级求调剂推荐 +11 adaie 2026-03-19 11/550 2026-03-20 21:30 by laoshidan
[考研] 一志愿西南交通 专硕 材料355 本科双非 求调剂 +5 西南交通专材355 2026-03-19 5/250 2026-03-20 21:10 by JourneyLucky
[考研] 320求调剂0856 +3 不想起名字112 2026-03-19 3/150 2026-03-19 22:53 by 学员8dgXkO
[考研] 材料工程专硕调剂 +5 204818@lcx 2026-03-17 6/300 2026-03-18 22:55 by 204818@lcx
[考研] 考研求调剂 +3 橘颂. 2026-03-17 4/200 2026-03-17 21:43 by 有只狸奴
信息提示
请填处理意见