请教octopussheng

蓝云天

你好：
     我有一个4维系统，用你空间里给出的最大值方法做分叉图，程序改正如下：
      function dy = Chen(t,y)
% Chen系统
% 系统微分方程：
%        dx/dt = a(y-x)
%        dy/dt = dx-xz+cy-u
%        dz/dt = xy-bz
%        du/dt = x+k
%    a=y(5)
%    d=y(6)
%    b=y(7)
%    c=y(8)
%    k=y(9)
dy=zeros(9,1);
dy(1)=y(5)*(y(2)-y(1));
dy(2)=y(6)*y(1)-y(1)*y(3)+y(8)*y(2)-y(4);
dy(3)=y(1)*y(2)-y(7)*y(3);
dy(4)=y(1)+y(9);
dy(5)=0;
dy(6)=0;
dy(7)=0;
dy(8)=0;
dy(9)=0;


function Lorenz_bifur_r_getmax
% 最大值法求解分岔图
clear all
t0=[0 100];%积分时间
%bifurcation
for k=linspace(-4,4,1000);  
    [t,y]=ode45('Chen',t0,[1;1;1;1;36;-16;28;3;k]);
    [Xmax]=getmax(y(:,1));
    plot(k,Xmax,'b','markersize',1)
    hold on
    clear Xmax
end

红色字体是我做的修改，运行后提示：??? One or more output arguments not assigned during call to ‘getmax.m (getmax)’

麻烦你帮我看一下是怎么回事？
其他的版主大侠，或研友有时间也欢迎帮忙看看。
谢谢。

[ 本帖最后由 蓝云天 于 2009-10-27 19:27 编辑 ]

蓝云天

getmax函数为：
function [Xmax] = getmax(y)
a=length(y);
j=1;
for i=(a-1)/2:a
    b=(y(i,1)-y(i-2,1))/2;
    c=(y(i,1)+y(i-2,1))/2-y(i-1,1);
    if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0
        Xmax(j)=y(i-1,1);
        j=j+1;
    elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)
        Xmax(j)=y(i-1,1)-b^2/(4*c);
        j=j+1;
    end
end

octopussheng

程序没有问题，我复制过来可以顺利运行。

蓝云天

是可以画出图来，但是，
第一：与我看得文章中的图像不符。
第二：运行后提示：??? One or more output arguments not assigned during call to ‘getmax.m (getmax)’

系统在文中参数下有一部分是混沌和超混沌的。但是显示的图像没有这一特征。是怎么回事呢

谢谢

octopussheng

对于这个现象，我只能说，getmax这个函数不是很好用，建议用庞加莱截面方法来计算你这个系统的分岔图。

蓝云天

你说我用下面这个能行么？

clear;
global k;
range=[-4:0.01:4];


for k=range
    k
    y0=[10,10,10,10];
   
    tspan=[0:0.05:200];
    [t,Y]=ode45('Chen1',tspan,y0);
    [Xmax]=getmax(Y(:,1));
    plot(k,Xmax,'b','markersize',1)
    hold on
    clear Xmax
   
end

octopussheng

结果如何？有没有提示错误？

但是我理解，本质上是没有区别的。

[ 本帖最后由 octopussheng 于 2009-10-27 21:38 编辑 ]

蓝云天

呵呵，还行，刚才试了试，可以跑过去，不过换成了Chen1函数。谢谢  octopussheng   版主的指点。
其实，我很想用poincare section方法做分叉图，理论上是最好的。但是，你空间里提供的程序中没弄明白那个截面选取是怎么回事。
唉！还是需要努力的学习啊，

再次谢谢  octopussheng   版主的指点和帮助

octopussheng

可以自己先尝试用Lorenz庞加莱截面的做法，做做你这个系统的庞加莱截面，然后再试试计算分岔图。