圆盘非线性微分方程组使用RK4计算的结果

kpgkpg

这样的方程难道MATLAB也无能？

方程组见附件1。
附件2，3是我写的MATLAB程序。使用MATLAB的ode45函数求解此方程组，在我的机器上运行长达4小时，始终不能计算出结果。


检查程序，应该没有什么问题。欢迎各位高手和群里的同行指点。
我的QQ：635359794
E-mail:kpg19850521@126.com

[ 本帖最后由 kpgkpg 于 2008-11-28 16:24 编辑 ]

无水1324

main
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%系统方程（圆盘不平衡）
%    x'' + 2ξx' + x +γ( x.^2 + y.^2) x = e*ω_bar.^2cosτ
%    y'' + 2ξy' + y +γ( x.^2 + y.^2) y = e*ω_bar.^2sinτ - G
%    G =g/(ω.^2*δ) ;τ=ω*t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%初始参数
%    2ξ = 0. 24
%    ω_bar = 2. 1
%      γ = 0. 5
%偏心量e :[0. 08 ～0. 18 ]
%e = 0. 1175、0. 1384 和0. 1750
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
keci=0.12;
omega_bar=2.1;
gama=0.5;
delta=0.1;
omega=2.0;
G=9.8/(omega.^2*delta);
% e=0.08:0.01:0.18;
t=0:500;

tao=omega*t;
initial=[0,0,0,0];
cf=[];
ccf=[];

e=[ 0.1175, 0.1384 ,0.1750];
for i=1:length(e)
%     if e(i)==0.1175 || e(i)==0.1384||e(i)==0. 1750
        coef=[2*keci,1,gama,e(i)*omega_bar.^2,-G];        
%        for t=0:0.1:500;%500个周期（omega）
            [t,f]=ode45('right',tao,initial,[],coef);
            cf=[cf;f];
%        ccf=[ccf;cf];
%   end
end

figure
title('PointCare 图');
plot(cf(2,:),f(1,:))
xlabel('位移')
ylabel('速度')



function f=right(t,x,flag,coef)
f=[x(2);
   -x(2)*coef(1)-coef(2)*x(1)-coef(3)*(x(1).^2+x(3).^2)*x(1)+coef(4)*cos(t);
   x(4);
   -x(4)*coef(1)-coef(2)*x(3)-coef(3)*(x(1).^2+x(3).^2)*x(1)+coef(4)*sin(t)-coef(5)];

[ 本帖最后由 无水1324 于 2008-11-28 16:55 编辑 ]

kpgkpg

程序值得稍加修改，问题取得初步进展。修改的程序，画出的图形参见附件。


图片总是太大，裁剪了好几次，看来只能上传这一张了，其他的大家运行下MATLAB自己可以画出来的

[ 本帖最后由 kpgkpg 于 2008-11-28 19:06 编辑 ]

kpgkpg

无水兄，请问怎么样把图片贴出来呢？

无水1324

下面有一个“发表新回复”你点那个之后就知道了有一个增加附件的就可以了

kpgkpg

轴心轨迹图

轴心轨迹图

jgwang

多谢二位的奉献

无水1324

是不是你的步长太大了，这个轴心轨迹很粗糙的

shipo507

文章中的其他图形又是如何作出来呢，比如波形图，最大lyapunow指数图，相平面图，轴心轨迹图，时域波形图和幅值谱图等。

shipo507

文章中的其他图形又是如何作出来呢，比如波形图，最大lyapunow指数图，相平面图，轴心轨迹图，时域波形图和幅值谱图等。

woshiaq

学习了~~~~~

octopussheng

t的步长取小一些，试试0.001.

fanhuasijin

无水1324 发表于 2008-11-28 16:52

                               
登录/注册后可看大图

main
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%系统方程（圆盘不平衡）

请问下，就在这个模型，如果我想画无量纲振动随转速变化的分岔图，要求转速是50：1900，怎么修改下程序呢，谢谢

hsfy919

回复 13 # fanhuasijin 的帖子

那要看你的模型了，参数，初值，方程，取值等等都可能要修改

fanhuasijin

回复 14 # hsfy919 的帖子

Jeffcott刚性转子轴承系统的动力学模型。其方程如下：
mx’’=-fx+mew2sin(wt)
my’’=-fy+mew2cos(wt)+mg
上采用无量纲处理，无量纲轴颈坐标：X=x/c,Y=y/c;X’=x’/(cw), Y’=y’/(cw^2);X’’=x’’/(cw), Y’’=y’’/(cw^2);无量纲时间：tau=wt;无量纲油膜力分量：Fx=fx/…, Fy=fy/…其余无量纲量（略）.
设z1=X,Z2=Y,Z3=X’,Z4=Y’. 其中“’”表示d/dtau, 则（式-1）用状态变量表示的无量纲形式为z1’= Z3,z2’=Z4,z3’=-Fx/M+p*sin(tau),z4=-Fy/M+p*cos(tau)+G1
我要画出无量纲振动随转速变化的分岔图，那是不是就意味着要有时间t和转速w这两个自变量了，如果有两个自变量还能算吗，如果可以的话可不可以指定些参考资料，
问题比较大，希望先给我指个方向，我再琢磨琢磨，谢谢