Python - 积分函数并绘制结果
Python - Integrating a function and plotting results
我正在尝试对伯努利梁方程进行数值求解并绘制结果。方程的一阶导数是斜率，二阶导数是挠度。我一步步接近这个问题，先绘制函数，然后对它进行积分，并将积分结果绘制在同一张图上。
到目前为止我的代码如下。我怀疑问题在于 integrate.quad returns 是单个值，而我正试图从中获取多个值。有谁知道如何处理它？
from scipy import integrate
import numpy as np
from pylab import *
# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100
def d2y_dx2(x):
return (-F*x)/(E*I)
a = 0.0
b = L
res, err = integrate.quad(d2y_dx2, a, b)
t = np.linspace(a,b,100)
ax = subplot(111)
ax.plot(t, d2y_dx2(t))
ax.plot(t, res(t))
show()
编辑：波纹管是修改后的代码，包含 willcrack 的答案。此代码现在有效，但结果不正确。在底部，我添加了使用正确的光束方程解析解绘制结果的代码。
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100
# Integration parameters
a = 0.0
b = L
# Define the beam equation
def d2y_dx2(x,y=None):
return (-F*x)/(E*I)
def something(x):
return integrate.quad(d2y_dx2)[0]
# Define the integration1 - slope
def slope(t):
slope_res = []
for x in t:
res1, err = integrate.quad(d2y_dx2, a, b)
slope_res.append(res1)
return slope_res
# Define the integration1 - deflection
def defl(t1):
defl_res = []
for t in t1:
res2, err = integrate.dblquad(d2y_dx2,a,b, lambda x: a, lambda x: b)
defl_res.append(res2)
return defl_res
# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
t = np.linspace(a,b,100)
t1 = np.linspace(a,b,100)
ax1.plot(t, d2y_dx2(t))
ax2.plot(t, slope(t))
ax3.plot(t1, defl(t1))
plt.show()
解析解、代码和结果如下。偏转光束的形状转过来，光束的末端在 x = 0.
from __future__ import division #to enable normal floating division
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
w = 10 #beam cross sec width (mm)
h = 10 #beam cross sec height (mm)
I = (w*h**3)/12 #cross sec moment of inertia (mm^4)
I1 = (w*h**3)/12
E = 200000 #steel elast modul (N/mm^2)
L = 100 #beam length(mm)
F = 100 #force (N)
# Define equations
def d2y_dx2(x):
return (-F*x)/(E*I)
def dy_dx(x):
return (1/(E*I))*(-0.5*F*x**2 + 0.5*F*L**2)
def y(x):
return (1/(E*I))*(-(1/6)*F*(x**3) + (1/2)*F*(L**2)*x - (1/3)*F*(L**3))
# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
a = 0
b = L
x = np.linspace(a,b,100)
ax1.plot(x, d2y_dx2(x))
ax2.plot(x, dy_dx(x))
ax3.plot(x, y(x))
plt.show()
也许你可以尝试这样的事情
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100
# Integration parameters
a = 0.0
b = L
# Define the beam equation
def d2y_dx2(x,y=None):
return (-F*x)/(E*I)
def something(x):
return integrate.quad(d2y_dx2)[0]
# Define the integration1 - slope
def slope(t):
slope_res = []
for x in t:
res1, err = integrate.quad(d2y_dx2, a, b)
slope_res.append(res1)
return slope_res
# Define the integration1 - deflection
def defl(t1):
defl_res = []
for t in t1:
res2, err = integrate.dblquad(d2y_dx2,a,b, lambda x: a, lambda x: b)
defl_res.append(res2)
return defl_res
# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
t = np.linspace(a,b,100)
t1 = np.linspace(a,b,100)
ax1.plot(t, d2y_dx2(t))
ax2.plot(t, slope(t))
ax3.plot(t1, defl(t1))
plt.show()
结果：
我想我找到了斜坡的解决方案。
稍后我会尝试另一个。
这是更新。
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100
# Integration parameters
a = 0.0
b = L
# Define the beam equation
def d2y_dx2(x,y=None):
return (-F*x)/(E*I)
# Define the integration1 - slope
def slope(x):
slope_res = np.zeros_like(x)
for i,val in enumerate(x):
y,err = integrate.quad(f,a,val)
slope_res[i]=y
return slope_res
# Define the integration1 - deflection
def defl(x):
defl_res = np.zeros_like(x)
for i,val in enumerate(x):
y, err = integrate.dblquad(d2y_dx2,0,val, lambda x: 0, lambda x: val)
defl_res[i]=y
return defl_res
# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
t = np.linspace(a,b,100)
t1 = np.linspace(a,b,100)
ax1.plot(t, d2y_dx2(t))
ax2.plot(t, slope(t))
ax3.plot(t1, defl(t1))
plt.show()
新结果：
还在为最后一个而苦苦挣扎...
我正在尝试对伯努利梁方程进行数值求解并绘制结果。方程的一阶导数是斜率，二阶导数是挠度。我一步步接近这个问题，先绘制函数，然后对它进行积分，并将积分结果绘制在同一张图上。
到目前为止我的代码如下。我怀疑问题在于 integrate.quad returns 是单个值，而我正试图从中获取多个值。有谁知道如何处理它？
from scipy import integrate
import numpy as np
from pylab import *
# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100
def d2y_dx2(x):
return (-F*x)/(E*I)
a = 0.0
b = L
res, err = integrate.quad(d2y_dx2, a, b)
t = np.linspace(a,b,100)
ax = subplot(111)
ax.plot(t, d2y_dx2(t))
ax.plot(t, res(t))
show()
编辑：波纹管是修改后的代码，包含 willcrack 的答案。此代码现在有效，但结果不正确。在底部，我添加了使用正确的光束方程解析解绘制结果的代码。
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100
# Integration parameters
a = 0.0
b = L
# Define the beam equation
def d2y_dx2(x,y=None):
return (-F*x)/(E*I)
def something(x):
return integrate.quad(d2y_dx2)[0]
# Define the integration1 - slope
def slope(t):
slope_res = []
for x in t:
res1, err = integrate.quad(d2y_dx2, a, b)
slope_res.append(res1)
return slope_res
# Define the integration1 - deflection
def defl(t1):
defl_res = []
for t in t1:
res2, err = integrate.dblquad(d2y_dx2,a,b, lambda x: a, lambda x: b)
defl_res.append(res2)
return defl_res
# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
t = np.linspace(a,b,100)
t1 = np.linspace(a,b,100)
ax1.plot(t, d2y_dx2(t))
ax2.plot(t, slope(t))
ax3.plot(t1, defl(t1))
plt.show()
解析解、代码和结果如下。偏转光束的形状转过来，光束的末端在 x = 0.
from __future__ import division #to enable normal floating division
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
w = 10 #beam cross sec width (mm)
h = 10 #beam cross sec height (mm)
I = (w*h**3)/12 #cross sec moment of inertia (mm^4)
I1 = (w*h**3)/12
E = 200000 #steel elast modul (N/mm^2)
L = 100 #beam length(mm)
F = 100 #force (N)
# Define equations
def d2y_dx2(x):
return (-F*x)/(E*I)
def dy_dx(x):
return (1/(E*I))*(-0.5*F*x**2 + 0.5*F*L**2)
def y(x):
return (1/(E*I))*(-(1/6)*F*(x**3) + (1/2)*F*(L**2)*x - (1/3)*F*(L**3))
# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
a = 0
b = L
x = np.linspace(a,b,100)
ax1.plot(x, d2y_dx2(x))
ax2.plot(x, dy_dx(x))
ax3.plot(x, y(x))
plt.show()
也许你可以尝试这样的事情
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100
# Integration parameters
a = 0.0
b = L
# Define the beam equation
def d2y_dx2(x,y=None):
return (-F*x)/(E*I)
def something(x):
return integrate.quad(d2y_dx2)[0]
# Define the integration1 - slope
def slope(t):
slope_res = []
for x in t:
res1, err = integrate.quad(d2y_dx2, a, b)
slope_res.append(res1)
return slope_res
# Define the integration1 - deflection
def defl(t1):
defl_res = []
for t in t1:
res2, err = integrate.dblquad(d2y_dx2,a,b, lambda x: a, lambda x: b)
defl_res.append(res2)
return defl_res
# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
t = np.linspace(a,b,100)
t1 = np.linspace(a,b,100)
ax1.plot(t, d2y_dx2(t))
ax2.plot(t, slope(t))
ax3.plot(t1, defl(t1))
plt.show()
结果：
我想我找到了斜坡的解决方案。 稍后我会尝试另一个。 这是更新。
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100
# Integration parameters
a = 0.0
b = L
# Define the beam equation
def d2y_dx2(x,y=None):
return (-F*x)/(E*I)
# Define the integration1 - slope
def slope(x):
slope_res = np.zeros_like(x)
for i,val in enumerate(x):
y,err = integrate.quad(f,a,val)
slope_res[i]=y
return slope_res
# Define the integration1 - deflection
def defl(x):
defl_res = np.zeros_like(x)
for i,val in enumerate(x):
y, err = integrate.dblquad(d2y_dx2,0,val, lambda x: 0, lambda x: val)
defl_res[i]=y
return defl_res
# Plot
fig, (ax1, ax2, ax3) = plt.subplots(3)
t = np.linspace(a,b,100)
t1 = np.linspace(a,b,100)
ax1.plot(t, d2y_dx2(t))
ax2.plot(t, slope(t))
ax3.plot(t1, defl(t1))
plt.show()
新结果：
还在为最后一个而苦苦挣扎...