单自由度体系分析是结构动力学和地震工程学习的启蒙知识，也是复杂结构动力分析的基石，具有重要的作用和研究价值。本文将探讨如何用Python的科学计算模块来实现单自由度体系的动力响应分析。

经典的弹性有阻尼单自由度体系在任意激励P(t)作用下的运动方程可写为

或

求解此运动方程的经典方法有杜哈梅积分法、Newmark-β法、中心差分法等。借助Python求解此方程最快捷的方法就是将方程降阶为状态空间形式调用scipy.integrate.odeint()函数。

上述方程的状态空间形式为

采用odeint()进行单自由度体系共振简谐激励和地震激励作用下响应求解的完整代码为

import scipy as sp
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def pd(tn,*p_args):
   """ 等时间间隔离散信号a对应的连续化函数（线性插值方法实现） """
   dt = p_args[0] # 离散信号时间间隔
   a  = p_args[1] # 离散信号数据

   ind = int(sp.floor(tn/dt))
   if (ind+1)>=len(a):
       return 0.0
   else:
       al = a[ind]; ar = a[ind + 1];
       k = (ar - al) / dt
       return al+k*(tn-ind*dt)

def solve_sdof_resonance(omg = 1.0*2.0*sp.pi, zeta = 0.02):

   y0 = sp.asarray([0.0,0.0]) # 初始条件
   p = lambda t: sp.sin(omg*t) # 共振简谐激励

   f_sdof = lambda y,t: sp.asarray([ y[1], -p(t)-2.0*zeta*omg*y[1]-omg*omg*y[0] ])

   t = sp.linspace(0, 30, 1001)
   y = odeint(f_sdof, y0, t)

   # 绘图:
   plt.figure("Resonance Response",(12,4))
   plt.plot(t,y[:,0])
   plt.grid(True)
   plt.show()

def solve_sdof_eqwave(omg = 1.0*2.0*sp.pi, zeta = 0.02):

   y0 = sp.asarray([0.0,0.0]) # 初始条件

   a = sp.loadtxt("EQ-S-3.txt") # 读取地震波
   dt = 0.005 # 时间间隔
   n = len(a)
   p = lambda t: pd(t,dt,a) # 共振简谐激励

   f_sdof = lambda y,t: sp.asarray([ y[1], -p(t)-2.0*zeta*omg*y[1]-omg*omg*y[0] ])

   t = sp.linspace(0.0,dt*(n-1),n)
   y = odeint(f_sdof, y0, t)

   # 绘图:
   plt.figure("Seismic Response",(12,8))
   plt.subplot(2,1,1)
   plt.plot(t,a,label=r"输入地震波加速度时程")
   plt.grid(True)
   plt.legend()
   plt.subplot(2,1,2)
   plt.plot(t,y[:,0],label=r"SDOF体系位移响应时程")
   plt.grid(True)
   plt.legend()
   plt.show()

if __name__ == '__main__':
   solve_sdof_resonance()
   solve_sdof_eqwave()

共振响应曲线

地震激励及响应曲线