Newmark法求振动方程，但是一直出现错误，这个程序是参考网上...
clc;clear;close;m=;
m=diag(m); %质量矩阵
k=[9.58*10^10 -3.15*10^10 0 0 0 0;
-3.15*10^10 9.43*10^10-6.28*10^10 0 0 0;
0 -6.28*10^10 9.11*10^10-2.83*10^10 0 0;
0 0 -2.83*10^10 9.11*10^10-6.28*10^100;
0 0 0 -6.28*10^10 1.073*10^11-4.45*10^10
0 0 0 0 -4.45*10^10 1.48*10^11];%刚度矩阵
c=[4*exp(6) -4*exp(6) 0 0 0 0 ;
-4*exp(6) 4*exp(6) 0 0 0 0;
0 0 exp(6) -exp(6) 0 0;
0 0 -exp(6) exp(6) 0 0;
0 0 0 0 0 0;
0 0 0 0 0 0]; %阻尼矩阵
f0=1.589*10^7;
dt=0.0005; %仿真步长
alfa=0.25;
beta=0.5;
tend=0.3
tt=0:dt:tend;
a0=1/alfa/dt/dt;
a1=beta/alfa/dt;
a2=1/alfa/dt;
a3=1/2/alfa-1;
a4=beta/alfa-1;
a5=dt/2*(beta/alfa-2);
a6=dt*(1-beta);
a7=dt*beta;
d=zeros(6,length(tt));
v=zeros(6,length(tt));
a=zeros(6,length(tt));
d(:,1)=';
v(:,1)=';
a(:,1)=';
for i=1:length(tt);
f(:,i)=';
end
for i=2:length(tt);
ke=k+a0*m+a1*c;
fe(:,i-1)=f(:,i)+m*(a0*d(:,i-1)+a2*v(:,i-1)+a3*a(:,i-1))+c*(a1*d(:,i-1)+a4*v(:,i-1)+a5*a(:,i-1));
d(:,i)=inv(ke)*fe(:,i);
a(:,i)=a0*(d(:,i)-d(:,i-1))-a2*v(:,i-1)-a3*a(:,i-1);
v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i);
end
plot(t,d)
页:
[1]