请教：如何画出位移随参数变化的分岔图（附有原始程序）
%%Part 1clc
clear all
fno=('shuchu.txt');
i=1;
global Ij;
fid=fopen(fno,'w');
for j=1030:1040
Ij=j;
=ode45(@yundongguiji,,);
=size(u);
uuu(:,i)=real(u(:,1));
i=i+1;
end
=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变化的分岔图，请问程序如何实现。
麻烦赐教！ 这么久也没人说，那我就自己说一下吧。
这些天看了论坛里面不少关于分岔与poincare映射的帖子，对相关概念的理解又上升了一层，觉得对自己的帮助很大。尤其是octopussheng 非常无私地将源码与经验介绍给大家；无水1324的回答简短但精炼；中原的解释往往一语中的......
我是搞水轮发电机组轴系耦联振动的，在大连理工大学读博。所研究方向涉及到了有限元、非线性转子动力学、不平衡电磁力、水力等相关领域的问题。
这段时间通过看书、看帖子发现我贴出来的程序是有问题的，第一部分可以不要。因为和分岔没有关系。我的系统应该称之为非自治系统，里面存在显含时间t项cos(omega*t)以及sin(omega*t)，激励频率应该为2*pi/omega。
按照octopussheng的帖子“非自治系统分岔图绘制实例”，应对Ij定义范围
function xxxx
clear;
global Ij;
omega=10;
Ij=;%Ij从1030变化到1120的范围，步长没取太大，可以调试
period=2*pi/omega;%%% 周期
%% 其余的部分应该同之前的帖子差不多了
k=0;
YY1=[];
step=2*pi/100;%步长
for F=Ij %%%%考虑的就是参数F对YY1的影响
y0=; %%%%初始值
Ij
k=k+1;
% discard the first 60 periodic data;
%除去前面60个周期的数据，并将最后的结果作为下一次积分的初值
tspan=;%%%%%取前60个周期
=ode45(@xxxx,tspan,y0);
y0=Y(end,:);%%%%取结果最后一行的值，也就是3个数值作为下一次积分的初值
j=1;
for i=60:200 %%%从第60个周期开始算到第200个周期
tspan=;
=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映射的同窗非常少，所以希望在这里同大家互相学习交流。 t,Y]=ode45(@xxxx,tspan,y0);
在所有的计算中参数Ij是全局变化的，然后你又赋值给了F，是不是F才是真正的需要全局化的变量呢，总感觉你这个参数传递有点问题 恩，当时帖子发的有点仓促。其实后来想一下不需要这样，把求解程序改为
%%%%%求解运动轨迹的分岔图
clc
clear
% global Ij
Ij=1050:1070;
for i=1:length(Ij)
disp('Ij(i)=');
disp(Ij(i));
period=2*pi/10;
=ode45(@yundongguiji,,);
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之后整百的点，那是不是浪费计算时间呢。如果只计算整周期的点，而画图时选择后面几十个周期的点是不是更好呢。
回复 地板 chunshui2003 的帖子
1、 =ode45(@yundongguiji,,);，这里只是初始参数的变化。2、你上面的程序就是“只计算整周期的点，而画图时选择后面几十个周期的点”
[ 本帖最后由 无水1324 于 2009-12-18 08:07 编辑 ] 恩，谢谢无水指点，又把困惑弄明白了一点。看来仔细研究还是有好处的。
页:
[1]