「Python与地震工程」单自由度体系求解之通用微分方程数值解法（基于scipy.integrate.odeint）
背景图：2008年5月汶川地震
原理
单自由度体系剖析是构造动力学和地震工程学习的启蒙常识，也是简单构造能源剖析的基石，具备重要的作用和钻研价值。本文将探讨如何用Python的科学计算模块来实现单自由度体系的能源响应分析。
经典的弹性有阻尼单自由度体系在任意激励P(t)作用下的静止方程可写为
$$m\ddot{x}\left( t \right) + c\dot{x}\left( t \right) + kx\left( t \right) = P\left( t \right)$$
或
$$\ddot{x}\left( t \right) + 2 \zeta \omega_0 \dot{x}\left( t \right) + \omega_0^2 x\left( t \right) = p\left( t \right)$$
求解此静止方程的经典办法有杜哈梅积分法、Newmark-法、核心差分法等。借助Python求解此方程最快捷（编程者的工作量小，并非最疾速）的办法就是将方程降阶为状态空间模式调用scipy.integrate.odeint()函数。须要留神的是，对于地震激励，须要把离散的地震波“包装”成一个连续函数能力被odeint辨认。
上述方程的状态空间模式为
$$\frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\begin{array}{c} {{y_1}}\\ {{y_2}} \end{array}} \right) = \left( {\begin{array}{c} {{y_2}}\\ {p\left( t \right) - 2\zeta \omega_0 {y_2} - {\omega_0^2}{y_1}} \end{array}} \right)$$
程序代码
采纳odeint()进行单自由度体系共振简谐激励和地震激励作用下响应求解的残缺代码为
import scipy as spfrom scipy.integrate import odeintimport matplotlib.pyplot as pltdef 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 draw_response(title, ta, a, t, u): plt.figure(title,(12,6)) plt.subplot(2,1,1) plt.plot(ta,a,label=r"输出地震波加速度时程") plt.grid(True) plt.legend() plt.xlim(0,t[-1]) plt.subplot(2,1,2) plt.plot(t,u,label=r"SDOF体系位移响应时程") plt.xlabel(r"工夫 (s)") plt.grid(True) plt.legend() plt.xlim(0,t[-1]) plt.show()def solve_sdof_eqwave_odeint(omg0 = 1.0*2.0*sp.pi, zeta = 0.05): y0 = sp.asarray([0.0,0.0]) # 初始条件 acc0 = sp.loadtxt("EQ-S-3.txt") # 读取地震波 dt = 0.005 # 工夫距离 n = len(acc0) t0 = sp.linspace(0.0, dt*(n-1), n) p = lambda t: pd(t, dt, acc0) # 线性间断化地震激励 f_sdof = lambda y,t: sp.asarray([ y[1], -p(t)-2.0*zeta*omg0*y[1]-omg0*omg0*y[0] ]) # 显示地震完结后一段时间内的自在振动衰减状况 ne = round(n*1.2) t = sp.linspace(0.0,dt*(ne-1),ne) y = odeint(f_sdof, y0, t) draw_response("Seismic Response -- scipy.integrate.odeint", t0, acc0, t, y[:,0]) # 绘图if __name__ == '__main__': solve_sdof_resonance() solve_sdof_eqwave()文中所用地震波下载：
EQ-S-3.txt
后果
转载本文请注明出处。
「Python与地震工程」单自由度体系求解之通用微分方程数值解法
科研成果中应用本文代码请援用作者课题组的钻研论文：
作者课题组发表的钻研论文列表（继续更新中……）