中心差分法Matlab源程序

clear n=100 %根据实际可任意取n值 M=rand(n) %假设M为n阶随机矩阵 C=rand(n) %假设M为n阶随机矩阵 K=rand(n) %假设M为n阶随机矩阵 Q1=rand(n,1) %假设Q1为nx1阶随机矩阵 d1=zeros(n,1) %假设d1为nx1阶零矩阵 v1=zeros(n,1) %假设v1为nx1阶零矩阵 a1=inv(M)*(Q1-C*v1-K*d1) %求得初始加速度 delta=0.28 %选取时间步长（可根据实际选取） c0=1./(delta^2) c1=1./(2.*delta) c2=2.*c0 c3=1./c2 d0=d1-delta*v1+c3*a1 %起步位移计算 M1=c0*M+c1*C i=1 for ti=0:delta:5 %选定时间范围（可任意） d(i,1:n)=d1 %存储各个时刻ti位移d(t) a(i,1:n)=a1 %存储各个时刻ti速度v(t) v(i,1:n)=v1 %存储各个时刻ti加速度a(t) t(i)=ti %存储各个时刻ti Q2=Q1-(K-c2*M)*d1-(c0*M-c1*C)*d0 %计算t时刻有效载荷 d2=inv(M1)*Q2 %求t+delta时刻位移 v1=c1*(-d0+d2) %计算t时刻速度 a1=c0*(d0-2*d1+d2) %计算t时刻加速度 d0=d1 d1=d2 i=i+1 end %这里将第一条位移、速度、加速度曲线画在同一面内 subplot(2,3,1),plot(t,d(:,1),'r','linewidth',2),xlabel('t'),ylabel('d1（t)') subplot(2,3,2),plot(t,v(:,1),'g','linewidth',2),xlabel('t'),ylabel('v1(t)') subplot(2,3,3),plot(t,a(:,1),'b','linewidth',2),xlabel('t'),ylabel('a1(t)') %因为这里质量矩阵M，刚度矩阵K，阻尼矩阵C随机选取，故出的图也是随机的