哪位同志有下图的MAPLE或MATLAB程序,先谢了!
哪位同志有下图的MAPLE或MATLAB程序,先谢了!3维ODE方程zhou1920@126.com
[ 本帖最后由 无水1324 于 2007-8-27 09:58 编辑 ] 你这个是什么图啊,光给一个图其他什么也没有,怎么编写程序给你
回复 #2 咕噜噜 的帖子
是分支图呀回复 #1 xueyongzhou 的帖子
写出方程,自己在maple和matlab中编写一个,与你作同样模型的人才可能有程序的 方程为dS/dt=r*S*(1-(S+I)/K)-b*S*IdI/dt=b*S*I-c*I-c1*I*y/(I+K1)
dy/dt=(a2-c2*I*y/(I+K2))*y
请大家帮忙!谢谢先.
[ 本帖最后由 xueyongzhou 于 2007-8-27 09:11 编辑 ]
回复 #4 无水1324 的帖子
程序没有通用性吗?回复 #6 xueyongzhou 的帖子
你这个可以参考一下论坛里面的一些分岔程序,分岔针对自己的系统可以采用不同的方法如频闪、最大值法等等,所以通用性的还是比较少的 你的系统虽说是3个自由度,可是方程并不算复杂,因此算好了方法编写程序并不困难
论坛有很多分叉程序,你可以参考,有什么具体问题可以上来提问 自己模仿liliangbiao的程序画的图.不好看.请帮忙修改!
function dx=abc(t,X)
x=X(1);
y=X(2);
z=X(3);
r=5;a2=1.2;c=1.5;c1=2;c2=1;K=3;k1=0.6;k2=0.5;
global beta;
dx=zeros(3,1);
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;
clear;
global beta;
range=;
k=0;
YY=[];
for beta=range
beta
y0=rand(1,3);
k=k+1;
tspan=;
=ode45('abc',tspan,y0);
count=find(t>100);
Y=Y(count,:);
% 画x的分岔图。
j=1;
n=length(Y(:,1));
for i=2:n-1
if Y(i-1,1)+eps<Y(i,1) && Y(i,1)>Y(i+1,1)+eps% 简单的取出局部最大值。
YY(k,j)=Y(i,1);
j=j+1;
end
end
if j>1
plot(beta,YY(k,),'k.','markersize',1);
end
hold on;
index(k)=j-1;
end
xlabel('beta');
ylabel('x');
%title('bifurcation diagram'); 你这好像没有发生分叉啊,可能是你的参数选取有问题,或者参数范围不合适 请问程序有问题吗?
但从相图上看在这范围有分支的.
[ 本帖最后由 xueyongzhou 于 2007-8-28 18:01 编辑 ]
回复 #11 xueyongzhou 的帖子
可是你的分叉图上可真是没有发生分叉啊,贴出你的相图 clc;clear allglobal r a2 beta c c1 c2 k k1 k2
r=5;a2=1.2;beta=6;c=1.5;c1=2;c2=1;k=3;k1=0.6;k2=0.5;
options=odeset('RelTol',1e-4,'AbsTol',);
=ode45('epidemic',,,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');
function dy = rigid(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);
取beta=5时稳定;
取beta=8时产生HOPF分支周期解!
回复 #13 xueyongzhou 的帖子
:@L 我说的是相图,不过贴出来程序也好,看了一下,发现你的相图和分叉图的系统微分方程不一样回复 #14 咕噜噜 的帖子
微分方程是一样的呀!
页:
[1]
2