哪位同志有下图的MAPLE或MATLAB程序,先谢了!

xueyongzhou

哪位同志有下图的MAPLE或MATLAB程序,先谢了!3维ODE方程
zhou1920@126.com


[ 本帖最后由 无水1324 于 2007-8-27 09:58 编辑 ]

咕噜噜

你这个是什么图啊，光给一个图其他什么也没有，怎么编写程序给你

xueyongzhou

是分支图呀

无水1324

写出方程，自己在maple和matlab中编写一个，与你作同样模型的人才可能有程序的

xueyongzhou

方程为dS/dt=r*S*(1-(S+I)/K)-b*S*I
         dI/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 编辑 ]

xueyongzhou

程序没有通用性吗?

无水1324

你这个可以参考一下论坛里面的一些分岔程序，
分岔针对自己的系统可以采用不同的方法如频闪、最大值法等等，所以通用性的还是比较少的

咕噜噜

你的系统虽说是3个自由度，可是方程并不算复杂，因此算好了方法编写程序并不困难
论坛有很多分叉程序，你可以参考，有什么具体问题可以上来提问

xueyongzhou

自己模仿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=[5:0.01:6];
k=0;
YY=[];
for beta=range
    beta
    y0=rand(1,3);
    k=k+1;
    tspan=[0:0.05:400];
    [t,Y]=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,[1:j-1]),'k.','markersize',1);
    end
    hold on;
    index(k)=j-1;
end
xlabel('beta');
ylabel('x');
%title('bifurcation diagram');

咕噜噜

你这好像没有发生分叉啊，可能是你的参数选取有问题，或者参数范围不合适

xueyongzhou

请问程序有问题吗?
但从相图上看在这范围有分支的.

[ 本帖最后由 xueyongzhou 于 2007-8-28 18:01 编辑 ]

咕噜噜

可是你的分叉图上可真是没有发生分叉啊，贴出你的相图

xueyongzhou

clc;clear all
global 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',[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');



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分支周期解!

咕噜噜

:@L 我说的是相图，不过贴出来程序也好，看了一下，发现你的相图和分叉图的系统微分方程不一样

xueyongzhou

微分方程是一样的呀!