matlab在地震工程的应用

我是叶子

最近在看徐赵东的书，上面有一个关系用阵型分解法求解地震响应的程序：
%mmain program%
%tish is for solving the dynamic responses %
%this program is for the frame structure%
xs=2*0.287;
ag=dzhbo*0.01*xs;
dt=0.02;
ndzh=400;
cn=3;
m0=[2.762 2.760 2.300]*1e+3;
k0=[2.485 1.921 1.522]*1e+4;
l=diag(ones(cn));
m=diag(m0);
[ik]=matrixju(k0,cn);
[x,d]=eig(ik,m);
d=diag(sqrt(d));
for i=1:cn
     [d1(i),j]=min(d);
     xgd(:,i)=x(:,j);
     d(j)=max(d)+1;
end
w=d1;
x=xgd;
a1=2*w(1)*w(2)*(0.05*w(2)-0.07*w(1))/(w(2)^2-w(1)^2);
a2=2*(0.07*w(2)-0.05*w(1))/(w(2)^2-w(1)^2);
for j=i:cn
     x(:,j)=x(:,j)/x(cn,j);
     znb0(j)=(a1+a2*w(j)^2)/2/w(j);
     zhcan(j)=(x(:,j))'*m*1/((x(:,j))'*m*x(:,j));
     [dlt(j,:),dltacceler(j,:)]=zxzj(znb0(j),w(j),ag);
end
for i=1:cn
     disp1=0;
     accel1=0;
     for j=1:cn
         disp0=zhcan(j)*dlt(j,:)*x(i,j);
         accel0=zhcan(j)*dltacceler(j,:)*x(i,j);
         disp1=disp1+disp0;
         accel1=accel1+accel0;
     end
     disp(i,:)=disp1;
     accel(i,:)=accel1;
end
t=0:dt:ndzh*dt;
subplot(2,2,1)
plot(t,disp(3,:)*1e+3,'k');
subplot(2,2,2)
plot(t,accel(3,:),'k')
     
%this sub-program is for matrix aggregation of system %
function [kcju]=matrixju(korc,cn)
kcju=zeros(cn);
for i=1:cn-1
     kcju(i,i)=korc(i)+korc(i+1);
     kcju(i,i+1)=-korc(i+1);
     kcju(i+1,i)=-korc(i+1);
end
kcju(cn,cn)=korc(cn)
%this sub-program is solving the dynamic response of single degree system%
function[bx,acceler]=zxzj(znb,w,dag)
dt=0.02;
n=400;
x(1)=0;
dx(1)=0;
ddx(1)=0;
s=1+znb*dt*w+w^2*dt*2/6;
for i=1:n
     a(i)=x(i)+dx(i)*dt+ddx(i)*dt^2/3;
     b(i)=dx(i)+ddx(i)*dt/2;
     ddx(i+1)=-1*(dag(i+1)+2*znb*w*b(i)+w^2*a(i))/s;
     dx(i+1)=b(i)+ddx(i+1)*dt/2;
     x(i+1)=a(i)+ddx(i+1)*dt^2/6;
end
bx=x;
acceler=ddx;
·在主程序中
for i=1:cn
     disp1=0;
     accel1=0;
     for j=1:cn
         disp0=zhcan(j)*dlt(j,:)*x(i,j);
         accel0=zhcan(j)*dltacceler(j,:)*x(i,j);
         disp1=disp1+disp0;
         accel1=accel1+accel0;
     end
     disp(i,:)=disp1;
     accel(i,:)=accel1;
end
总感觉 accel0=zhcan(j)*dltacceler(j,:)*x(i,j);这步是不是有问题呀，x(i,j)应该改成x(j,i)
望做过该例子的朋友提提意见

ChaChing

专业方面, 个人水平有限不懂!
不过就编程来说, 可以猜测撰写者是学fortran的
没能利用到matlab的优点

我是叶子

ChaChing 发表于 2013-5-3 21:07

                               
登录/注册后可看大图

专业方面, 个人水平有限不懂!
不过就编程来说, 可以猜测撰写者是学fortran的
没能利用到matlab的优点

谢谢，你觉得那些地方可以优化？

ChaChing

我是叶子 发表于 2013-5-3 21:10

                               
登录/注册后可看大图

谢谢，你觉得那些地方可以优化？

个人水平有限, 仅能建议LZ看下这些帖
From http://home.chinavib.com/blog-63979-18250.html
1.[原创] 提高matlab运行速度的一点心得 http://www.chinavib.com/thread-25305-1-1.html
  [原创] 提高matlab运行速度的一点心得(之二) http://www.chinavib.com/thread-26449-1-1.html
  [原创]提高matlab运行速度和节省空间的一点心得(之三) http://www.chinavib.com/thread-33199-1-1.html
...

我是叶子

非常感谢那你的帮助

江山美人

感觉这个程序所得结果还是有问题