传递矩阵法matlab程序

tic
%！传递矩阵法！
clc;clear;

%参数输入：
n = 20;%建筑总层数
h = 3.6;%H/n;%层高（m）
H = n*h;%建筑总高（m）
m = 23437*ones（n1）;%各层沿高度的分布质量（kg/m）
M = 84375*ones（n1）;%各层集中质量（kg）
EI = 2.6447*10^12*ones（n1）;%各层抗弯刚度（N/m2）
GA = 4.0387*10^10*ones（n1）;%各层抗剪刚度（N/m2）
v = 0.02;%阻尼
g = 9.8;%重力加速度（m/s2）
k_ = 1.2;%剪应力不均匀系数


%********************************************地震反应计算*******************************************************
%参数输入：
EI_ = EI*（1+v*1i）;
GA_ = GA*（1+v*1i）;
%傅里叶变换
   load a3
   T=max（a（:1））;
   a（:2） = a（:2）/1000;%统一单位m
   Ag=fft（a（:2））;
   N=length（a）;
  %计算omega
   omega=zeros（N1）;
   for j=1:N
     omega（j）=2*pi/T*j;
   end
   %对Ag进行处理
  for j=1:N                
      if j>=2&&j<=N/2
          Ag（j）=2*Ag（j）;
      end
      if j>=N/2+2&&j<=N
          Ag（j）=0;
      end
  end
  
%参数计算:
alpha = k_*h*h*（omega.^2*（m./GA_）‘）;%alpha（:n）表示第n层的所有omega对应的值，下同
beta = h^4*（omega.^2*（m./EI_）‘）;%beta（N*n）
lamda1 = （（alpha.^2/4+beta）.^0.5 - alpha/2）.^0.5;%lamda1（N*n）
lamda2 = （（alpha.^2/4+beta）.^0.5 + alpha/2）.^0.5;%lamda2（N*n）

d = ones（N1）*（h*h*EI_.^-1）‘;%d（N*n）
deta = （lamda1.^2 + lamda2.^2）.^-1;%deta（N*n）
c0 = deta.*（lamda2.^2.*cosh（lamda1） + lamda1.^2.*cos（lamda2））;
c1 = deta.*（lamda2.^2./lamda1.*sinh（lamda1） + lamda1.^2./lamda2.*sin（lamda2））;
c2 = deta.*（cosh（lamda1） - cos（lamda2））;
c3 = deta.*（sinh（lamda1）./lamda1 - sin（lamda2）./lamda2）;

%场传递矩阵tf（4*4*N*n）
tf（11::） = c0 - alpha.*c2;
tf（12::） = -beta.*c3/h;
tf（13::） = -beta.*c2./d;
tf（14::） = beta.*（alpha.*c3 - c1）./d/h;
tf（21::） = -h*（c1 - alpha.*c3）;
tf（22::） = c0;
tf（23::） = h*beta.*c3./d;
tf（24::） = beta.*c2./d;
tf（31::） = -d.*c2;
tf（32::） = d.*c1/h;
tf（33::） = c0;
tf（34::） = beta.*c3/h;
tf（41::） = h*d.*（alpha.*c1 - （beta+alpha.^2）.*c3）./beta;
tf（42::） = d.*c2;
tf（43::） = h*（c1 - alpha.*c3）;
tf（44::） = c0 - alpha.*c2;

%点传递矩阵tp（4*4*N*n）
tp（11::） = ones（Nn）;
tp（22::） = ones（Nn）;
tp（33::） = ones（Nn）;
tp（44::） = ones（Nn）;
tp（12::） = zeros（Nn）; 
tp（12::） = zeros（Nn）;   
tp（13::） = zeros（Nn）;   
tp（14::） = -（omega.^2）*M‘ ;
tp（21::） = zeros（Nn）;   
tp（23::） = zeros（Nn）;
tp（24::） = zeros（Nn）;   
tp（31::） = zeros（Nn）;   
tp（32::） = zeros（Nn）;   
tp（34::） = zeros（Nn）; 
tp（41::） = zeros（Nn）;   
tp（42::） = zeros（Nn）;   
tp（43::） = zeros（Nn）;

%求总传递矩阵Tn
Tn_ = 1;
for i=1:N
    for j=1:n
    Tn_ = tp（::ij）*tf（::ij）*Tn_;
    %Tn（::ij） = Tn_; 
    end
    TN（::i） = Tn_;
    Tn_ = 1;
end

%求地震反应kesi
Ug = -Ag./（omega.^2）;%傅里叶变换后的地面位移
temp = （TN（24:）.*TN（12:） - TN（14:）.*TN（22:））./（TN（22:）.*TN（11:） - TN（12:）.*TN（21:））;
Q_I（:1） = temp（:）.*Ug（:）;
%第一层层底剪力
temp = （TN（14:）.*TN（21:） - TN（24:）.*TN（11:））./（TN（22:）.*TN（11:） - TN（12:）.*TN（21:））;
M_I（:1） = temp（:）.*Ug（:）;
%第一层层底弯矩


kesi_I（1:1） = Q_I（:1）; 
kesi_I（2:1） = M_I（:1）; 
kesi_I（3:1） = zeros（N1）;
k