应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 计算模拟 » matlab处理常微分方程作图问题
81/1返回列表
查看: 1266  |  回复: 7
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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语言。
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2013-02-25 13:02:03
已阅   回复此楼   关注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 的主题更新
81/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭