请教：ode45迭代一直buzy，算不了

chunshui2003

%%%%% 微分方程
function y=fangcheng(t,x)    % 求解定值问题

m=20;             %转子质量
e=0.0009;         %圆盘偏心距
xi=0.13;
xi_z=0.15;
omega_0=16*2*pi;      
omega_e=22*2*pi;      
delta=0.0018;     
kc=25*10^5;       %定子的径向刚度
f=0.1;            %转子与定子之间的摩擦系数
R=0.1;            %碰摩点到转子形心的距离
Fe=9;             %轴向推力的幅值
omega_z=20*2*pi;
g=9.80;           %重力加速度
global omega   

Vr=R*omega+(x(2).*x(5)-x(3).*x(4))/sqrt(x(2).^2+x(4).^2);
x(1)=t;

y=zeros(7,1);
y(1)=t;
y(2)=x(3);
y(3)=e*omega.^2*cos(omega.*x(1))-2.*xi*omega_0.*x(3)-omega_0^2.*x(2)-kc/m*(1-delta/sqrt(x(2).^2+x(4).^2)).*(x(2)-f*Vr/sqrt(Vr^2+x(7).^2).*x(4));
y(4)=x(5);
y(5)=e*omega.^2*sin(omega.*x(1))-2.*xi*omega_0.*x(5)-omega_0^2.*x(4)-kc/m*(1-delta/sqrt(x(2).^2+x(4).^2)).*(x(4)+f*Vr/sqrt(Vr^2+x(7).^2).*x(2));
y(6)=x(7);
y(7)=Fe/m*cos(omega_e*x(1))-g-2*xi_z*omega_z*x(7)-omega_z^2*x(6)-kc/m*f*(sqrt(x(2)^2+x(4)^2)-delta)*x(7)/sqrt(Vr^2+x(7)^2);


%%%%%%%

clc
clear
global omega   
tic

omega=2.51*16*2*pi;
period=2*pi/omega;
tspan=[0:period/100:10*period];
y0=[1 1 1 1 1 1 1 ];
[t,x]=ode45(@fangcheng,tspan,y0);

toc

问题出现在y(7)行的红色部分，之前将方程写错了。“*x(7)”写成了“-x(7)”，但是可以算。结果改过来“*x(7)”反倒是不能算了，一直提示buzy，开始以为是初值的问题，改为了y0=[10 10 10 10 10 10 10 ]也还是运行不了，不知道是哪里错了，请大家给指点一下。

chunshui2003

请大家帮忙看一下啊，我又试了ode113、23、23s和15s，只有15s能够运行。但结果不让人满意，并且出现提示：Warning: Failure at t=1.615971e-001.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.741082e-016) at time t.
我又把时间t的步长给进一步缩小至10^17，但提示“Maximum variable size allowed by the program is exceeded.”
真的头疼了，当看到希望时又无法前进感觉有点痛苦，希望各位友人能够给解答一下，谢谢了！

无水1324

方程应该是刚性的了，所以只有15s可以同，因为它可以求解刚性方程

chunshui2003

请问无水用其他的语言吗，比如fortran或者C。因为方程如果是刚性的，并且matlab的ode15s又不理想，那么在fortran里面计算是不是效果更好一点呢？另外，当涉及到循环画分岔图的时候，matlab效率很低，上次算了一个分岔图，循环了3000次，每次里面又有140个小循环，居然算了50分钟。
如果无水遇到这种问题，会选择用fortran替代matlab进行前期计算而用matlab进行后处理吗？

无水1324

呵呵，如果是刚性，就不是语言的问题了，而是算法的问题。你可以搜索一下google，关于刚性方程的求解 ，现在都有很多文章在研究。市一个比较麻烦的问题。

另外画分岔图慢，你可以考虑用Fortran来计算，matlab来处理，这是一个很好的策略

chunshui2003

恩，谢谢无水。本来打算在这学期结束前写出一篇论文，但越学越发现需要学的东西太多。也好，遇到问题解决问题是一个很难受但又快乐的过程。很高兴能向你讨教，另外无水的研究方向是什么呢？

无水1324

做齿轮的

chunshui2003

除了用能量法之外，你用有限元吗？另外，非线性油膜力和密封力有涉及吗？如果有的话，以后希望可以交流一下。

无水1324

没有做有限元的

yangxiaokang

我最近也在做转轴非线性油膜力和密封力方面的，不过是刚刚起步，希望今后能够多多指点，谢谢:@)

chunshui2003

恩，可以互相学习交流。

yangxiaokang

原帖由 chunshui2003 于 2010-1-9 10:50 发表

                               
登录/注册后可看大图

恩，可以互相学习交流。

请多多指教，我也在做转子碰摩，发了一个帖子请帮我看看。你的我也看过几遍，不过小弟才疏学浅，没看明白
http://www.chinavib.com/forum/thread-89514-1-1.html