声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

12
返回列表 发新帖
楼主: xueyongzhou

[论坛公告] 哪位同志有下图的MAPLE或MATLAB程序,先谢了!

[复制链接]
发表于 2007-8-28 21:19 | 显示全部楼层
:@L看似一样,你在那个分叉图里面写的改成相图里面的试试看,这个写法

dx(1)=r*x*(1-x/K-y/K)-beta*x*y;
dx(2)=beta*x*y-c*y-c1*z*y/(y+k1);
dx(3)=(a2-c2*z/(y+k2))*z;
应该是错误的

dy(1)=r.*y(1).*(1-y(1)/k-y(2)/k)-beta.*y(1).*y(2);
dy(2)=beta.*y(1).*y(2)-c.*y(2)-c1.*y(3).*y(2)/(y(2)+k1);
dy(3)=(a2-c2.*y(3)/(y(2)+k2)).*y(3);
才是正确地
我不知道你的最大值法算得分叉图居然能得出结果

[ 本帖最后由 咕噜噜 于 2007-8-28 21:21 编辑 ]
回复 支持 反对
分享到:

使用道具 举报

 楼主| 发表于 2007-8-28 21:23 | 显示全部楼层
那该怎么办呀?
 楼主| 发表于 2007-8-30 11:20 | 显示全部楼层
问题仍未解决?请高手指点!
发表于 2007-8-30 11:37 | 显示全部楼层

回复 #18 xueyongzhou 的帖子

你的两个程序的主程序呢,贴出来我帮你运行一下
 楼主| 发表于 2007-8-30 15:05 | 显示全部楼层

回复 #19 咕噜噜 的帖子

程序分别在上面的#9和#13,谢谢咕噜噜
发表于 2007-8-30 15:20 | 显示全部楼层
[t,y]=ode45('epidemic',[0 100],[0.5 1 0.8],options);


function dy = rigid(t,y)
看看你10楼的程序有没有给错啊,子程序名对不上
另,用最大值法运行了,和你的结果一样,你那样的编写程序其实也可以,但是变量多时就容易出问题最好不要用

[ 本帖最后由 咕噜噜 于 2007-8-30 15:21 编辑 ]
 楼主| 发表于 2007-8-30 21:16 | 显示全部楼层

回复 #21 咕噜噜 的帖子

是贴错了.下面是对的.请问咕噜噜对我的模型有没有其他较合适的方法得到分支图?

function dy =epidemic(t,y)
global r a2 beta c c1 c2 k k1 k2
dy=zeros(3,1);    % a column vector
dy(1)=r.*y(1).*(1-y(1)/k-y(2)/k)-beta.*y(1).*y(2);
dy(2)=beta.*y(1).*y(2)-c.*y(2)-c1.*y(3).*y(2)/(y(2)+k1);
dy(3)=(a2-c2.*y(3)/(y(2)+k2)).*y(3);


clc;clear all
global r a2 beta c c1 c2 k k1 k2
r=5;a2=1.2;beta=7;c=1.5;c1=2;c2=1;k=3;k1=0.6;k2=0.5;
options=odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-4]);
[t,y]=ode45('epidemic',[0 100],[0.5 1 0.8],options);
figure(1);
plot(t,y(:,1),'b',t,y(:,2),'r');
figure(2);
plot(t,y(:,3),'g');
figure(3);
plot3(y(:,1),y(:,2),y(:,3),'r');
发表于 2007-8-31 07:49 | 显示全部楼层

回复 #22 xueyongzhou 的帖子

目前通常用的就是最大值法和相图法
 楼主| 发表于 2007-8-31 07:54 | 显示全部楼层
请问相图法你有类似的程序吗?
发表于 2007-8-31 08:20 | 显示全部楼层

回复 #24 xueyongzhou 的帖子

你搜索一下论坛,咱们版块就有
还有,你最大值法画的图,你从那里看出来beta=8时又分叉呢,你给的程序里面plot3(y(:,1),y(:,2),y(:,3),'r');这个命令画出来的图什么意思
发表于 2007-8-31 08:22 | 显示全部楼层

回复 #23 咕噜噜 的帖子

相图法是什么方法,第一次听说这个名词?
 楼主| 发表于 2007-8-31 08:29 | 显示全部楼层

回复 #25 咕噜噜 的帖子

plot3(y(:,1),y(:,2),y(:,3),'r');这个命令画出来的图是相图.
 楼主| 发表于 2007-8-31 08:36 | 显示全部楼层
帮我运行一下下面程序吧.
clear
clc
close all
global beta;
for beta=5:0.01:7;
    [T,X]=ode45('lorenz1',[0:0.5:500],[10 10 10]);
    data=X(:,1);
    n=length(data);
    N=round(n/2);
    f=data(N-2);
    g=data(N-1);
    for i=N:n
        if g>=f & g>=data(i);
            plot(beta,g);
            hold on;
            axis on;
            axis([5 7 0 10]);
        else
        end
        f=g;
        g=data(i);
    end
end





function lorenz=lorenz1(t,x)
global a; a=10;
global b; b=8/3;
global c; c=28;
lorenz=[0;0;0];
lorenz(1)=a*(x(2)-x(1));
lorenz(2)=c*x(1)-x(1)*x(3)-x(2);
lorenz(3)=x(1)*x(2)-b*x(3);
发表于 2007-8-31 08:36 | 显示全部楼层

回复 #27 xueyongzhou 的帖子

:@L :@L 这个是相图吗?那你能不能告诉我你这里面的y(:,1),y(:,2),y(:,3),分别表示什么?通常相图应该是位移x和与其相对应的速度dx的动点轨迹
发表于 2007-8-31 08:55 | 显示全部楼层

回复 #28 xueyongzhou 的帖子

这个程序中午帮你运行好吧?我刚运行了一下看上去需要很多时间
我电脑都快死机了,我正在用电脑所以有点不方便,中午一定帮你运行
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-5-19 17:29 , Processed in 0.057918 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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