帮忙改个程序吧，谢谢

huerhui88

m1 ¨x 1 + c1Ûx 1 + k ( x1 - x2) =Fx ( x1 , y1 , Ûx 1 , Ûy1) ,
m1 ¨y1 + c1Ûy1 + k ( y1 - y2) =Fy ( x1 , y1 , Ûx 1 , Ûy1) - m1 g ,
m2 ¨x 2 + c2Ûx 2 + k ( x2 - x1) + k ( x2 - x 3) =m2 bω2cos (ωt) ,
m2 ¨y2 + c2Ûy2 + k ( y2 - y1) + k ( y2 - y3) =m2 bω2sin (ωt) - m2 g ,
m1 ¨x 3 + c1Ûx 3 + k ( x3 - x2) =Fx ( x3 , y3 - y4 , Ûx 3 , Ûy3 - Ûy 4) ,
m1 ¨y3 + c1Ûy3 + k ( y3 - y2) =Fy ( x3 , y3 - y4 , Ûx 3 , Ûy3 - Ûy4) - m1 g ,
m3 ¨y4 + csÛy4 + ksy4 =- Fy ( x3 , y3 - y4 , Ûx 3 , Ûy3 - Ûy 4) - m3 g ·
%运动微分方程
function d=fun(t,y,w)
N=length(y);
w=2100;
m1=4;%两端滑动轴承处等效集中质量
m2=32.1; %转子圆盘等效集中质量
m3=50.0;%轴承支座处等效集中质量
g=9.81;
e=0.00005; %偏心距
k=2.5e7;%弹性轴刚度
delta2=0.6e-3;%初始间隙
c1=1050;%转子圆盘处阻尼系数
c2=2100;%转子在轴承处阻尼系数
k1=7.5e7;
k2=2.5e9;
cb1=350;
cb2=500;
ox1=y(1);%未松动端竖直方向位移x1
oy1=y(2);%未松动端竖直方向位移y1
odx1=y(8);
ody1=y(9);
ox2=y(3);%圆盘位移x2
oy2=y(4);%圆盘位移y2
odx2=y(10);
ody2=y(11);
ox3=y(5);%松动端轴心位移x3
oy3=y(6);%松动端轴心位移y3
odx3=y(12);
ody3=y(13);
oy4=y(7);%质量m3在竖直方向位移y4
ody4=y(14);
if oy4<0
    cb=cb2;
    kb=k2;
elseif  (oy4>=0)&(oy4<=delta2)
        cb=0;
        kb=0;
    else
         cb=cb1;
         kb=k1;
end
M=[  m1   0     0    0    0     0    0;
     0    m1    0    0    0     0    0;
     0    0     m2   0    0     0    0;
     0    0     0    m2   0     0    0;
     0    0     0    0    m1    0    0;
     0    0     0    0    0     m1   0;
     0    0     0    0    0     0    m3;];
C=[  c1   0     0    0    0     0    0;
     0    c1    0    0    0     0    0;
     0    0     c2   0    0     0    0;
     0    0     0    c2   0     0    0;
     0    0     0    0    c1    0    0;
     0    0     0    0    0     c1   0;
     0    0     0    0    0     0    cb;];
  
K=[  k   0     -k    0    0     0    0;
     0    k    0    -k    0     0    0;
     -k   0    2*k    0    -k    0    0;
     0    -k   0    2*k    0    -k    0;
     0    0    -k    0     k    0    0;
     0    0    0    -k     0    k    0;
     0    0    0     0     0    0    kb;];
fx=oilx( ox1, oy1, odx1, ody1, w);
fy=oily( ox1, oy1, odx1, ody1, w);
fx1=oilx( ox3,oy3-oy4,odx3,ody3-ody4,w);
fy1=oily( ox3,oy3-oy4,odx3,ody3-ody4,w);

F=[   fx;
      fy-m1*g;
      m2*e*w^2*cos(t);
      m2*e*w^2*sin(t)-m2*g;
      fx1;
      fy1-m1*g;
      -fy1-m3*g ];
   
C=C/w;
K=K/w^2;
F=F/c/w^2;
for i=1:1:N/2
    y1(i,1)=y(i);
    y2(i,1)=y(i+N/2);
end

yy2=inv(M)*(F-C*y2-K*y1);
d=zeros(N,1);
for i=1:1:N/2
    d(i)=y2(i,1);
    d(i+N/2)=yy2(i,1);
end

%x方向油膜力
function oilforce=oilx(x,y,dx,dy,wi)
R=0.025; L=0.012; miu=0.018;  dfai=1; dert=0.00011;
e=sqrt(x*x+y*y);
delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
ppp1=(y+2.0*dx)/(x-2.0*dy);
sign1=sign(ppp1);
        
ppp2=y+2.0*dx;
sign2=sign(ppp2);
alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));
f1=sqrt(abs(abs((x-2.0*dy)*(x-2.0*dy))+abs((y+2.0*dx)*(y+2.0*dx))))/(1.0-abs(x*x)-abs(y*y));
fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);
oilforce=fx*delta;

%y方向油膜力
function oilforce=oily(x,y,dx,dy,wi)
R=0.025; L=0.012; miu=0.018;  dfai=1; dert=0.00011;
e=sqrt(x*x+y*y);
delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
ppp1=(y+2.0*dx)/(x-2.0*dy);
sign1=sign(ppp1);
      
ppp2=y+2.0*dx;
sign2=sign(ppp2);
alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));
f1=sqrt(abs(abs((x-2.0*dy)*(x-2.0*dy))+abs((y+2.0*dx)*(y+2.0*dx))))/(1.0-abs(x*x)-abs(y*y));
fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);
oilforce=fy*delta;

%主分析程序
clear all
h=pi/256;
w=2100;
tf=300000*2*pi/w;
tspan=0:h:tf;
y0=[0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5,0.05,0.5];
options=odeset('RelTol',10^-6,'AbsTol',10^-6);  
[t,y]=ode45(@fun,tspan,y0);
figure(1)
subplot(2,2,1);
plot(t,y(:,1),'r',t,y(:,2),'b')
title('未松动端位移响应');xlabel('t');ylabel('x1-red,y1-blue');
subplot(2,2,2);
plot(t,y(:,3),'r',t,y(:,4),'b')
title('圆盘位移响应');xlabel('t');ylabel('x2-red,y2-blue');
subplot(2,2,3);
plot(t,y(:,5),'r',t,y(:,6),'b')
title('松动端轴心位移响应');xlabel('t');ylabel('x3-red,y4-blue');
subplot(2,2,4);
plot(t,y(:,7),'b')
title('m3在竖直方向位移响应');xlabel('t');ylabel('y4');
figure(2)
subplot(2,2,1);
plot(y(:,1),y(:,2))
title('未松动端轴心轨迹');xlabel('x1');ylabel('y1');
subplot(2,2,2);
plot(y(:,5),y(:,6))
title('松动端轴心轨迹');xlabel('x3');ylabel('y3');
yy1=sqrt(y(:,1).^2+y(:,2).^2);%yy1=sqrt(x1^2+y1^2)
yy2=sqrt(y(:,5).^2+y(:,6).^2);%yy2=sqrt(x3^2+y3^2)
N1=256;fs=1024;stepf=fs/N1;k=pi*2/882.5056;%882.5056转子固有频率sqrt(k/m2)
n=0:stepf*k:(fs/2-stepf)*k;
YY1=abs(fft(yy1));
YY2=abs(fft(yy2));
figure(3)
plot(n,abs(YY1(1:N1/2)));grid on;
title('未松动情况下幅值谱');xlabel('频率比');ylabel('FFT幅值');
figure(4)
plot(n,abs(YY2(1:N1/2)));grid on;
title('松动情况下幅值谱');xlabel('频率比');ylabel('FFT幅值');
      
fclose('all');

dianxin05

这么长的程序怎么帮你改啊》？？
建议你自己设置断点，进行单步调试！看看每步的结果是否是自己想要的！

ChaChing

请问楼主若别人程序这麽长, LZ会帮看吗?
建议楼主看下本版规则并加强发问题方式!
待高人路过!

yufeng

什么地方错了呢？不会让我们运行吧