kpgkpg 发表于 2008-11-28 16:18

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

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

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


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

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

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

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 :
%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=;
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=;      
%      for t=0:0.1:500;%500个周期(omega)
            =ode45('right',tao,initial,[],coef);
            cf=;
%      ccf=;
%   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 发表于 2008-11-28 18:46

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


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

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

kpgkpg 发表于 2008-11-28 19:04

无水兄,请问怎么样把图片贴出来呢?

无水1324 发表于 2008-11-28 19:27

回复 板凳 kpgkpg 的帖子

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

kpgkpg 发表于 2008-12-2 18:36

轴心轨迹图

轴心轨迹图

jgwang 发表于 2008-12-26 20:42

学习学习

多谢二位的奉献

无水1324 发表于 2008-12-29 11:44

回复 6楼 kpgkpg 的帖子

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

shipo507 发表于 2009-1-4 10:01

其他图形如何作出

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

shipo507 发表于 2009-1-4 10:01

回复 楼主 kpgkpg 的帖子

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

woshiaq 发表于 2010-12-26 17:35

学习了~~~~~

octopussheng 发表于 2010-12-27 13:02

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

fanhuasijin 发表于 2011-3-23 09:53

无水1324 发表于 2008-11-28 16:52 static/image/common/back.gif
main
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%系统方程(圆盘不平衡)


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

hsfy919 发表于 2011-3-23 15:30

回复 13 # fanhuasijin 的帖子

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

fanhuasijin 发表于 2011-3-23 17:33

回复 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这两个自变量了,如果有两个自变量还能算吗,如果可以的话可不可以指定些参考资料,
问题比较大,希望先给我指个方向,我再琢磨琢磨,谢谢
页: [1] 2
查看完整版本: 圆盘非线性微分方程组使用RK4计算的结果