matlab实现数值解中心差分法 程序源码.zip

%两自由度中心差分法 dzb=xlsread('D:\Elcentrobo.xls','sheet1'); dzb1=dzb(1:800,2); p=0.01*2*0.287*dzb1; m1=2.7*10^5; m2=2.6*10^5; k1=8.56*10^7; k2=9.6*10^7; a1=0.05; a2=0.05; m=[m1 0;0 m2]; k=[k1+k2 -k2;-k2 k2]; [x,d]=eig(k,m); w1=sqrt(d(1,1)); w2=sqrt(d(2,2)); t=0.02; A1=(1/t^2+2*w1*a1/(2*t)); B1=(w1^2-2/(t^2)); C1=(1/t^2-2*w1*a1/(2*t)); A2=(1/t^2+2*w2*a2/(2*t)); B2=(w2^2-2/(t^2)); C2=(1/t^2-2*w2*a2/(2*t)); b1=x(:,1); b2=x(:,2); r1=(b1'*m*[1;1])/(b1'*m*b1); r2=(b2'*m*[1;1])/(b2'*m*b2); f1=-r1*p; f2=-r2*p; u1(1,1)=0; u1(2,1)=0; for i=2:800 p1(i-1,1)=f1(i-1,1)-B1*u1(i,1)-C1*u1(i-1,1); u1(i+1,1)=p1(i-1,1)/A1; end u2(1,1)=0; u2(2,1)=0; for i=2:800 p2(i-1,1)=f2(i-1,1)-B2*u2(i,1)-C2*u2(i-1,1); u2(i+1,1)=p2(i-1,1)/A2; end for i=1:801 y1(i,1)=x(1,:)*[u1(i,:);u2(i,:)]; y2(i,1)=x(2,:)*[u1(i,:);u2(i,:)]; end o=0:0.02:16; plot(o,y1) o=0:0.02:16; plot(o,y2) y11=abs(y1); y111=max(max(y11)) [rows,cols]=find(y11==max(max(y11))) y22=abs(y2); y222=max(max(y22)) [rows,cols]=find(y22==max(max(y22)))