|
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
我现在在做Wilson-θ法，程序编出来有问题！十层的结构，就是当质量由五次方改为4次方的时候，结构响应图的幅值就变得很大！好几百的次方！~~查看了程序也没查出错来:
clear;
clc;
m=25000*speye(10,10);
k=1.00E+08*sparse([1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8,9,9,9,10,10],...
[1,2,1,2,3,2,3,4,3,4,5,4,5,6,5,6,7,6,7,8,7,8,9,8,9,10,9,10],...
[2.8426,-2.0833,-2.0833,4.1667,-2.0833,-2.0833,3.75,-1.6667,-1.6667,3.3333,-1.6667,-1.6667,3.3333,-1.6667,...
-1.6667,2.52,-0.8533,-0.8533,1.7067,-0.8533,-0.8533,1.7067,-0.8533,-0.8533,1.7067,-0.8533,-0.8533,0.8533]);
k=full(k);
f1=k/m;
f=full(f1);
[x,d]=eig(f);
w=sqrt(d);
h1=0.05; %阻尼比
h2=0.05; %阻尼比
time=1:1:10
plot(time,x(1));
a1=(2*w(1,1)*w(2,2)*(h1*w(2,2)-h2*w(1,1)))/(w(2,2)^2-w(1,1)^2);
a2=2*(h2*w(2,2)-h1*w(1,1))/(w(2,2)^2-w(1,1)^2);
c=a1*m+a2*k;
c=full(c);
dzbo; %地震波的调用
ag1=a(:,2);
xs=70/max(abs(ag1));
ag=ag1*0.01*xs;
ag=ag';
seita=1.4;
dt=0.02;
T=10.0;
N1=ones(1,10);
N=N1';
n=T/dt+1;
k1=k+6*m/(seita^2*dt^2)+3*c/(seita*dt);
u=zeros(10,n);
v=zeros(10,n);
ju=zeros(10,n);
for j=2:n
G1=-m*N*seita*(ag(j)-ag(j-1));
G2=m*(6*v(:,j-1)/seita/dt+3*ju(:,j-1));
G3=c*(3*v(:,j-1)+seita*dt*ju(:,j-1)/2);
dp=G1+G2+G3;
du1=inv(k1)*dp;
du=du1/seita;
dju1=6.0*du1/(seita^2*dt^2)-6.0*v(:,j-1)/(seita*dt)-3.0*ju(:,j-1);
dju=dju1/seita;
du=dt*v(:,j-1)+dt^2*ju(:,j-1)/2.0+dt^2*dju;
dv=dt*ju(:,j-1)+dt*dju/2.0;
u(:,j)=u(:,j-1)+du;
v(:,j)=v(:,j-1)+dv;
ju(:,j)=-N*ag(j)-inv(m)*k*u(:,j)-inv(m)*c*v(:,j);
end
t=0:0.02:10;
plot(t,u(10,:));
不知哪位高手可帮忙解决下！急！
谢谢！！！
[ 本帖最后由 jlde_639 于 2007-9-25 10:39 编辑 ] |
|
发表于 2007-9-25 10:24
提升卡
变色卡
显身卡
楼主
发表于 2007-9-26 16:07
