from sympy import * import numpy as np import matplotlib.pyplot as plt import math #基本参数 g=9.8 l=60 f=2 Ac=8.15 Ic=3622000*(10e-8) E=2.06*10e3 m=1435.2 #方程运算 H0=(m*g*(l)**2)/(8*f) He=20502900 v1=1-(H0/He) v2=1-H0/(4*He) v3=1-H0/(9*He) p=l**2/(8*f) L=l*(1+(16/3)*(f/l)**2) a21=8*E*Ac*l a22=(m*((np.pi)**2)*(p)**2)*L a2=(a21/a22)*10e4 c2=(2*(np.pi)/l)**2*sqrt((E*Ic*v2)/m) c2=c2*1000#平衡系数 c1=((np.pi/l)**2)*sqrt((E*Ic*v1)/m) c1=c1*1000#平衡系数 c3=((3*np.pi/l)**2)*sqrt((E*Ic*v1)/m) c3=c3*1000 w=symbols('w') G=solve([w**4-(c1**2+c3**2+(10/9)*a2)*w**2+(c1**2*c3**2+a2*(c1**2/9+c3**2))],[w]) print(v1,v2,v3,a21,a22,a2) print(c2,c1,c3) print(G,c2)
以抛物线型拱为例,采用传统的数值计算方法过于繁琐(此题以李国豪版桥梁结构稳定与振动为例)
考虑拱桥的前几阶振型必须先考虑前几阶振型所对应的频率
(以三阶为例,一个为基频另外两个为对称频率)
对于此联立方程不只能求出振型比率
此外思考拱桥抛物线的矢高与拱桥长度的关系
采用最小2乘法拟合分析
根据上式代码计算出的结果
import matplotlib.pyplot as plt import numpy as np from scipy.optimize import leastsq import math plt.rcParams['font.sans-serif']=['SimHei']#显示中文 plt.rcParams['axes.unicode_minus']=False X=np.array([2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10]) Y1=np.array([24.52,24.61,24.68,24.73,24.76,24.788,24.81,24.827,24.842,24.855,24.865,24.875,24.883,24.89,24.896,24.902,24.907]) Y2=np.array([14.75,17.9,21.05,24.17,27.23,30.214,33.086,35.81,38.352,40.662,42.7,44.44,45.877,47.035,47.955,48.683,49.26]) def func1(params,x): A,B=params return A*np.sqrt(1-B*(1/x)) def error1(params,x,y): return func1(params,x) - y def slovePara1(): p0=[1,1] Para = leastsq(error1,p0,args=(X,Y1)) return Para def func2(params,x): k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12=params return 0.235702260395516*np.sqrt((k1+k2*x+k3*x**2+k4*x**3+(k5+k6*x+k7*x**2+k8*x**3+k9*x**4+k10*x**5+k11*x**6)**0.5)/(x*(k12*x**2+1))) def error2(params,x,y): return func2(params,x) - y def slovePara2(): p0=[1,1,1,1,1,1,1,1,1,1,1,1] Para = leastsq(error2,p0,args=(X,Y2)) return Para def solution(): Para = slovePara1() A, B, =Para[0] plt.scatter(X, Y1, color='green',label=u'已知样点',linewidth=0.3) x=np.linspace(2,10,50) y1=A*np.sqrt(1-B*(1/x)) plt.plot(x,y1,color='red',label='w2',linewidth=2) plt.legend() Para = slovePara2() k1, k2, k3,k4,k5,k6,k7,k8,k9,k10,k11,k12 =Para[0] x=np.linspace(2,10,50) y2=0.235702260395516*np.sqrt((k1+k2*x+k3*x**2+k4*x**3+(k5+k6*x+k7*x**2+k8*x**3+k9*x**4+k10*x**5+k11*x**6)**0.5)/(x*(k12*x**2+1))) plt.plot(x,y2,color='blue',label='w1',linewidth=2) plt.legend() plt.scatter(X, Y2, color='green',linewidth=0.3) plt.title(u'w2与w1随f的变化情况') plt.show() solution()
所画图形如下