精细积分法的算例
对于一个单自由度的体系，质量，刚度，阻尼比为。初始状态下体系的处于静止状态，当时刻作用一个大小为的简谐荷载，作用时间为，求体系的位移时程，并将结果与采用解析积分法作比较。积分的基本原理是将任意载荷视为一系列脉冲激励的叠加，系统的响应可以表示成对这些脉冲激励的响应函数的叠加。在单自由度系统中，积分的形式为：其中，是系统在时间的位移响应，是外载荷函数，是质量，是无阻尼自然频率，为有阻尼时结构的频率，是阻尼比.importscipy.integrateasintegrateimportmatplotlib.pyplotaspltimportnumpyasnpdefcompute_t(n,H,eta):##计算矩阵TN=20tau=eta/(2**N)H_tau1=H*tauH_tau2=H@H_tau1*tauH_tau3=H@H_tau2*tauH_tau4=H@H_tau3*tauT_a=H_tau1+H_tau2/2.0+H_tau3/6.0+H_tau4/24.0foriinrange(N):T_a=2.0*T_a+T_a@T_aT=np.eye(n)+T_areturnTdefprecise_integration_solver(H,F,ug,step,dt):#精细积分法invH=np.linalg.inv(H)T=compute_t(2,H,dt)x=np.zeros((step,2))xi=np.zeros(2)foriinrange(step-1):r0=-F*(-ug[i])r1=-F*(-ug[i+1]+ug[i])/dtxi=T.dot(xi+invH.dot(r0+invH.dot(r1)))-invH.dot(r0+invH.dot(r1)+dt*r1)x[i,:]=xiu=x[:,0]returnu#杜哈梅积分函数deffunc(tau):returnnp.sin(tau)*np.exp(-zeta*omg*(t-tau))*np.sin(omg*(t-tau))if__name__=="__main__":plt.style.use("ggplot")plt.rcParams['font.sans-serif']=['MiSans']#这个字体能够支持中英文字符plt.rcParams['axes.unicode_minus']=False#正常显示负号omg=1.0#自振频率zeta=0.05#固有阻尼比dt=0.005#时间间隔t0=np.arange(0.0,40.0,dt)step=len(t0)dpm=np.zeros(step)ug=np.sin(t0)#简谐荷载H=np.array([[0.0,1.0],[-omg*omg,-2.0*zeta*omg]])F=np.array([0.0,1.0])u=precise_integration_solver(H,F,ug,step,dt)foriinrange(step):t=t0[i]res=integrate.quad(func,0,t)[0]#调用scipy计算杜哈梅积分dpm[i]=1/omg*resfig,(ax0,ax1)=plt.subplots(nrows=2,ncols=1,sharex=True,)ax0.plot(t0,ug)ax0.set(title='简谐荷载')ax0.grid()ax1.plot(t0,u,label='精细积分法')ax1.plot(t0,dpm,linestyle='--',color='#1661ab',label='杜哈梅积分法')ax1.set(title='位移时程')ax1.grid()ax1.legend(frameon=True,#启用边框framealpha=0.8,#设置透明度fancybox=False,#边框圆角与否shadow=True,#是否启用阴影facecolor='lightgrey',#背景颜色edgecolor='dimgrey',#边框颜色loc='upperleft')fig.savefig("11.png",dpi=300)plt.show()来源：数值分析与有限元编程
