一维动态有限差分matlab程序

reference:

http://wenku.baidu.com/link?url=KE3pz7qs2J5gxE2h-CR7rl0uIg3NyozFZVdOMTnCcOkamQxfIoHhATt6yEbWvR1qDbUm_FonuQh5CBSmtS9jhX85-WmbRuvknz3t2LHatI3

网上很多关于一个节点的一维动态有限差分matlab程序，但是如果要做多个节点呢（也即是一维），就没有多少文章，今天改了个代码，看起来还行，但是节点不能太多(<100)

fd_dyna.m

function [u,v,ac]=fd_dyna(m,c,k,u0,v0,all_time,P0,dt)
 %结构动力学 中心差分法
 
  %m=质量；k=刚度；c=阻尼；u0=初始位移；v0=初始速度；all_time=所用时间；P0=荷载幅值；dt=时间步长；
  %u=位移；v=速度；ac=加速度；ek=等效刚度；p=荷载；ep=等效荷载；t=时间；
% method 2 for input (jimmy comment)%
% m=input('输入质量矩阵M ：');
% c=input('输入阻尼矩阵C：'); 
% k=input('输入刚度矩阵K：'); 
% u0=input('输入初始位移 u0： ');
% v0=input('输入初始速度 v0： ');
% all_time=input('输入模拟时间 time：');
% dt=input('输入时间步长dt ：');
% P0=input('输入load ：');
% [mm,nn]=size(k);
% n=all_time/dt;     %计算步数 
% u=zeros(mm,floor(n)+1);      %设定存储位移矩阵
% v=zeros(mm,floor(n)+1); %设定存储速度矩阵
% ac=zeros(mm,floor(n)+1);%设定存储加速度矩阵 
% P=zeros(mm,floor(n)+1);%设定存储荷载矩阵
%jimmy comment%




%  1d bar from heat.m (jimmy add for)%
q=0.065;
Z=1;                 % the length of bar
K=2.;               % conductivity
n=50;
delz=Z/(n-1);
z=0:delz:Z;
a=ones(1,n)*-2;
b=ones(1,n-1);
A=diag(a,0)+diag(b,-1)+diag(b,1);
A(1,1) = 1;         % adjust the A matrix for the top BC and for the gradient BC at the Bottom
A(1,2) = 0;
A(2,1) = -1;         % gradient BC
A(2,2) = 1;


 p=zeros(n,n);;      


C=zeros(n,1);       % make the BC vector
C(1) = 0;       %BC
C(2) = delz*q/K;
%  1d bar from heat.m%


 %A0=input('请按格式和顺序输入初始矩阵，如A0=[m,k,c,u0,v0,all_time,P0,dt]=[9240 1460000 6410 0.05 0 30 7300*3 0.05]=');   
 Init=[1 1 1 0.05 0 30 100*3 0.005]
 m=Init(1,1);k=Init(1,2);c=Init(1,3);u0=Init(1,4);v0=Init(1,5);all_time=Init(1,6);P0=Init(1,7);dt=Init(1,8); 
  
  if dt>2*sqrt(m/k)                %判断时间步长是否满足稳定条件    
      disp('不满足稳定条件:dt<=Tn/pi,请重新输入符合稳定条件的时间步长dt')    
      return     
  elseif 0<dt<=2*sqrt(m/k)    
      disp('满足稳定条件为：dt<=Tn/pi') 
  end 
 t=[0:dt:all_time];                             %将时间分步，采用等时间步长；
  [mm,nn]=size(t);                               %计算t的向量长度，得出步数；
  u=zeros(n,nn);                              %设定存储u的矩阵；
  v=zeros(n,nn);                              %设定存储v的矩阵；
  ac=zeros(n,nn);   
  u(1,2)=u0;                                     %赋值向量第2项为u0； 
  v(1,2)=v0;                                     %赋值向量第2项为v0；
%   ac(:,:,2)=(P0-c*v(:,2)-k*u(:,:,2))./m;              %求出初始加速度ac0；
  ac(1,2)=(P0-c*v(1,2)-k*u(1,2))/m;              %求出初始加速度ac0；
  u(1,1)=u(1,2)-dt*v(1,2)+((dt)^2)*ac(1,2)/2;    %计算初始条件u-1项；
  ek=m/(dt^2)+c/(2*dt);                          %计算等效刚度； 
  % a=k-(2*m)/(dt^2); 
  a=A-diag(ones(1,n)*(2*m)/(dt^2),0); 
  b=diag(ones(1,n)*(m/(dt^2)-c/(2*dt)),0);         %计算方程系数； 
  p(1,2)=P0*sin(0);                              %给出初始荷载条件；
  ep(:,2)=p(:,2)-a*u(:,2)-b*u(:,1);              %计算初始等效荷载；
  u(:,3)=ep(:,2)/ek;                             %计算位移u1=u(:,3) 
  for i=3:nn                                     %从第二项开始进行中心差分法计算；   
     
        p(1,i)=P0*sin(.5*pi*(i-2)*dt)*exp(-10*(i-2));              %给出荷载条件，按照简谐荷载计算；  
    
      ep(:,i)=p(:,i)-a*u(:,i)-b*u(:,i-1);        %计算等效荷载；     %-----------------------得出所需要结果----------------------------------%  
      u(:,i+1)=ep(:,i)/ek;                       %计算位移量；   
      v(:,i)=(u(:,i+1)-u(:,i-1))/(2*dt);          %计算速度量；  
      ac(:,i)=(u(:,i+1)-2*u(:,i)+u(:,i-1))/(dt^2);%计算加速度量； 
  end 
  
  % z axis response (jimmy add)  %
  %  z=z';u=u(1:n,:);v=v(1:n,:);ac=ac(1:n,:);p=p(1:n,:);ep=(1:n,:);
  %------------------------绘制位移、速度、加速度z-axis曲线-----------------------% 
  %plot(t,u,'b-o'),hold on,plot(t,v,'g--p'),hold on,plot(t,ac,'r:x'),grid on,xlabel('时间(s)'),ylabel('位移(m)速度(m/s)加速度(m/s^2)'),title('顶层u,v,ac的时程曲线'); 
 for ii=1:100:all_time/dt
 % subplot(10,1,1);
  plot(z,u(:,ii),'b-');
  %grid,xlabel('时间(s)');
  %ylabel('位移(m)');
  %title('位移u的时程曲线');legend('位移u')
 
  %subplot(3,1,2),plot(z,v(:,1),'k'),grid,xlabel('时间(s)'),ylabel('速度(m/s)'),title('速度v的时程曲线');legend('速度v') 
  %subplot(3,1,3),plot(z,ac(:,1),'r'),grid,xlabel('时间(s)'),ylabel('加速度(m/s^2)') ,title('加速度ac的时程曲线');legend('加速度ac') 
 end 
  % 1 node time response( jimmy comment)
%   t=t(:,1:end-1);u=u(:,2:end-1);v=v(:,2:end);ac=ac(:,2:end);p=p(:,2:end);ep=ep(:,2:end);
%   %------------------------绘制位移、速度、加速度时程曲线-----------------------% 
%   %plot(t,u,'b-o'),hold on,plot(t,v,'g--p'),hold on,plot(t,ac,'r:x'),grid on,xlabel('时间(s)'),ylabel('位移(m)速度(m/s)加速度(m/s^2)'),title('顶层u,v,ac的时程曲线'); 
%   subplot(3,1,1),plot(t,u,'b-'),grid,xlabel('时间(s)'),ylabel('位移(m)'),title('位移u的时程曲线');legend('位移u')
%   subplot(3,1,2),plot(t,v,'k'),grid,xlabel('时间(s)'),ylabel('速度(m/s)'),title('速度v的时程曲线');legend('速度v') 
%   subplot(3,1,3),plot(t,ac,'r'),grid,xlabel('时间(s)'),ylabel('加速度(m/s^2)'),title('加速度ac的时程曲线');legend('加速度ac')