请教：如何画出位移随参数变化的分岔图（附有原始程序）

chunshui2003

%%Part 1

clc
clear all
fno=('shuchu.txt');
i=1;
global Ij;
    fid=fopen(fno,'w');              
for j=1030:1040
    Ij=j;                           
    [t,u]=ode45(@yundongguiji,[0:0.01:20],[1 1 1 1 1 1 1 1]);                  
    [n m]=size(u);                                 
   
    uuu(:,i)=real(u(:,1));  

    i=i+1;                                       
end         
   
[n m]=size(uuu);

    for j=1:m
      
        fprintf(fid,'%d(R)     ',1029+j);
    end
for k=1:n

   for j=1:m
        %fprintf(fid,'%f     %f   ',uuu(k,2*j-1), uuu(k,2*j));  
        %fprintf(fid,'%f ',uuu(k,2*j-1));
        fprintf(fid,'%f ',uuu(k,j));
    end
    fprintf(fid,'\n');                          
end
    status=fclose(fid);                       
   
   
%% Part 2
function uu=yundongguiji(t,u)

%求解微分中涉及到的一些计算参数

m1=1.5*10^4;               
m2=1.1*10^4;               
c1=6.0*10^4;            
c2=7.0*10^4;              
k1=8.5*10^7;              
k2=6.5*10^7;               
k3=3.5*10^7;               
delta0=8*10^-3;            
e2=0.3*10^-3;            
mu0=4*pi*10^-7;           
kj=5.2;                  
R=1.2;                     
L=0.5;                 
omega=10;        

K1=25*k1+9*k2+k3;
K2=-5*k1+3*k2+3*k3;
K3=k1+k2+9*k3;
e1=0.5*10-3;           
B1=e1*omega^2*cos(omega*t)/delta0;
B2=e1*omega^2*sin(omega*t)/delta0;
B3=e2*omega^2*cos(omega*t)/delta0;
B4=e2*omega^2*sin(omega*t)/delta0;

global Ij;              

F0=R*L*pi*kj^2*Ij^2*mu0/(m1*delta0^3);
Z1=(1-sqrt(1-u(1).^2-u(3).^2))/sqrt(u(1).^2+u(3).^2);
Z2=sqrt(u(5).^2+u(7).^2)/sqrt(u(1).^2+u(3).^2);

uu=zeros(8,1);
uu(1)=u(2);
uu(2)=B1+F0./(1-u(1).^2-u(3).^2)*(Z1+Z1.^3+Z1.^5).*u(1)./sqrt(u(1).^2+u(3).^2)-c1./m1.*u(2)-u(1).*(K1+K2.*Z2)./(16.*m1);
uu(3)=u(4);
uu(4)=B2+F0./(1-u(1).^2-u(3).^2)*(Z1+Z1.^3+Z1.^5).*u(3)./sqrt(u(1).^2+u(3).^2)-c1./m1.*u(4)-u(3).*(K1+K2.*Z2)./(16.*m1);
uu(5)=u(6);
uu(6)=B3-c2./m2.*u(6)-u(5).*(K3+K2.*Z2)./(16.*m2);
uu(7)=u(8);
uu(8)=B4-c2./m2.*u(8)-u(7)*(K3+K2.*Z2)./(16.*m2);

   
看着很长，其实就两部分。
将第二部分写成.m文件，运行第一部分。即得到每一个Ij（从1030到1040）对应的一组u结果。   
  现在希望得到u(1)随Ij变化的分岔图，请问程序如何实现。
麻烦赐教！

chunshui2003

这么久也没人说，那我就自己说一下吧。
这些天看了论坛里面不少关于分岔与poincare映射的帖子，对相关概念的理解又上升了一层，觉得对自己的帮助很大。尤其是octopussheng 非常无私地将源码与经验介绍给大家；无水1324的回答简短但精炼；中原的解释往往一语中的......

我是搞水轮发电机组轴系耦联振动的，在大连理工大学读博。所研究方向涉及到了有限元、非线性转子动力学、不平衡电磁力、水力等相关领域的问题。

这段时间通过看书、看帖子发现我贴出来的程序是有问题的，第一部分可以不要。因为和分岔没有关系。我的系统应该称之为非自治系统，里面存在显含时间t项cos(omega*t)以及sin(omega*t)，激励频率应该为2*pi/omega。

按照octopussheng的帖子“非自治系统分岔图绘制实例”，应对Ij定义范围

function xxxx
clear;
global Ij;
omega=10;

Ij=[1030:0.1:1120];  %Ij从1030变化到1120的范围，步长没取太大，可以调试
period=2*pi/omega;  %%% 周期
%% 其余的部分应该同之前的帖子差不多了
k=0;
YY1=[];
step=2*pi/100;  %步长
for F=Ij     %%%%考虑的就是参数F对YY1的影响
    y0=[1 1 1 1 1 1 1 1 ];    %%%%初始值
    Ij
    k=k+1;
    % discard the first 60 periodic data;
    %除去前面60个周期的数据，并将最后的结果作为下一次积分的初值
    tspan=[0:step:60*period];  %%%%%取前60个周期
    [t,Y]=ode45(@xxxx,tspan,y0);
    y0=Y(end,:);  %%%%取结果最后一行的值，也就是3个数值作为下一次积分的初值
    j=1;
    for i=60:200 %%%从第60个周期开始算到第200个周期
        tspan=[i*period:step:(i+1)*period];
        [t,Y]=ode45(@xxxx,tspan,y0);
        YY1(k,j)=Y(end,1);   % get the omega data from every period end
        j=j+1;               %取出每一个周期内的第一个解的最后一个值。
        y0=Y(end,:);
    end
end
bifdata=YY1(:,end-51:end);
plot(Ij,bifdata,'k.','markersize',1);

这里面我把微分方程的定义m.文件省略了。

不知道这样的改进是否正确，请教octopussheng、无水及中原。
另外，实际上运行octopussheng的程序并不能得到帖子中的图，请问是哪里出现了问题。
由于懂分岔与poincare映射的同窗非常少，所以希望在这里同大家互相学习交流。

无水1324

t,Y]=ode45(@xxxx,tspan,y0);

在所有的计算中参数Ij是全局变化的，然后你又赋值给了F，是不是F才是真正的需要全局化的变量呢，总感觉你这个参数传递有点问题

chunshui2003

恩，当时帖子发的有点仓促。其实后来想一下不需要这样，把求解程序改为
%%%%%求解运动轨迹的分岔图
clc
clear
% global Ij

Ij=1050:1070;
for i=1:length(Ij)
    disp('Ij(i)=');
    disp(Ij(i));
   
    period=2*pi/10;
    [t,u]=ode45(@yundongguiji,[0:period/100:100*period],[1 1 1 1 1 1 1 1 Ij(i)]);
    plot(Ij(i),real(u(6000:100:end,1)),'k','markersize',2);
    hold on;
    xlabel('Ij/A');
    ylabel('d\x\dt');
end

而Ij在整体运动微分方程中作为u(9)，以初值的形式传递进去即可。
请问无水，如果这样做的话是不是应该就没有问题。
另外，又想到每次计算10000个点，而画图的时候只取6000之后整百的点，那是不是浪费计算时间呢。如果只计算整周期的点，而画图时选择后面几十个周期的点是不是更好呢。

无水1324

1、 [t,u]=ode45(@yundongguiji,[0:period/100:100*period],[1 1 1 1 1 1 1 1 Ij(i)]);，这里只是初始参数的变化。
2、你上面的程序就是“只计算整周期的点，而画图时选择后面几十个周期的点”

[ 本帖最后由 无水1324 于 2009-12-18 08:07 编辑 ]

chunshui2003

恩，谢谢无水指点，又把困惑弄明白了一点。看来仔细研究还是有好处的。