Python - 用 scipy 求解伯努利梁方程 [英] Python - solving Bernoulli's beam equation with scipy
问题描述
回答问题的过程已经在下面链接的问题中开始，但该主题专门关于集成功能，已回答.所以我添加了一个新问题.
近似解决方案(种类):下面的代码是由 willcrack 制作的.形状看起来比上一个问题中的好，但值仍然不正常.
from scipy import 集成将 numpy 导入为 np导入 matplotlib.pyplot 作为 plt# 梁参数L = 100w = 10小时 = 10我 = (w*h**3)/12E = 200000F = 100# 积分参数a = 0.0b = L# 定义梁方程def d2y_dx2(x,y=None):返回 (-F*x)/(E*I)# 定义积分1 - 斜率定义斜率(x):斜率_res = np.zeros_like(x)对于 enumerate(x) 中的 i,val:y,err = 积分.quad(f,a,val)斜率_res[i]=y返回斜率_res# 定义积分1 - 偏转定义 defl(x):defl_res = np.zeros_like(x)对于 enumerate(x) 中的 i,val:y, err =integrate.dblquad(d2y_dx2,0,val, lambda x: 0, lambda x: val)defl_res[i]=y返回 defl_res# 阴谋图, (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, 斜率(t))ax3.plot(t1, defl(t1))plt.show()您正在对一个微分方程进行积分，您在循环中计算定积分的方法可以说是次优的.
Scipy 中的标准方法是使用
OP 报告了一个问题，理解 solve_ivp(lambda x, Y: [Y[1], x], ....
我们有一个标准形式的一阶 ODE 系统
y₁' = f₁(x, y₁(x), …, yₙ(x))… = …yₙ' = f₁(x, y₁(x), …, yₙ(x))可以写成，用大写字母表示向量
Y' = F(x, Y(x))求解微分方程组solve_ipv 正是需要这个F(x, Y) 函数.
代替 lambda 表达式，您可以编写如下函数定义，这可能更容易解释
def F(x, Y):# 解包函数值的向量y = y[0]h = Y[1]# 计算导数dy_over_dx = hdh_over_dx = x# 返回导数向量返回 [dy_over_dx, dh_over_dx]s = solve_ivp(F, …)在答案中简洁(过于简洁?)表示为 lambda x,Y:[Y[1],x] ...
The process of answering the question has already started in the question on the link bellow, but that topic was specifically about integrating a function, which was answered. So I added a new question.
Python - Integrating a function and plotting results
THE PROBLEM: how to solve a beam equation y''(x) = M(x) / (E*I) using scipy integrate.
SOLUTION, courtesy of gboffi:
#---------- DESCRIPTION
# cantilever beam with point load P at the free end
# original beam equation: y''(x) = M(x)/(E*I)
# moment equation: M(x) = -P*x
# x goes from the free end to the clamped end
# we have a second order diff eq: y''(x) = x
# we implement a new function:
# h = y',
# h' = y'' = M(x) = x
# we get a system of two ODE of first order
# y' = h
# h' = x
# we write the equations in vector form
# Y' = F(x, Y(x)) = F(x,Y)
# we define a function that returns the original values
#----------- CODE
from __future__ import division
from numpy import linspace
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# Exact solution, E*Iy = const, y1 = y', y0 = y,
w = 10 #beam cross sec width (mm)
h = 10 #beam cross sec height (mm)
Iy = (w*h**3)/12 #cross sec moment of inertia (mm^4)
E = 200000 #steel elast modul (N/mm^2)
L = 100 #beam length(mm)
P = 100 #point load (N)
x = linspace(0, L, 51)
y1 = (-P/(2*E*Iy))*x**2+(P*L**2)/(2*E*Iy)
y0 = (-P/(6*E*Iy))*x**3+((P*L**2)/(2*E*Iy))*x-(2*P*L**3)/(6*E*Iy)
# Define the vector function for E=const for integration
def F(x,Y):
#unpack the vector function
y = Y[0]
h = Y[1]
#compute the derivatives
dy_dx = h
dh_dx = (-P/(E*Iy))*x
#return the vector of derivatives values
return [dy_dx, dh_dx]
# Numerical solution
s = solve_ivp(
F, # Y[0]=y0, Y[1]=y1, dy0dx=y1, dy1dx=x
[L, 0.0], # interval of integration (NB: reversed, because...)
[0.0, 0.0], # initial conditions (at the 1st point of integ interval)
t_eval=linspace(L, 0, 101) # where we want the solution to be known
)
# Plotting
fig, (ax1, ax2) = plt.subplots(2)
ax1.plot(x, y0, label="Exact y")
ax2.plot(x, y1, label="Exact y'")
ax1.plot(s.t[::2], s.y[0][::2], label="Numeric y", linestyle='', marker='.')
ax2.plot(s.t[::2], s.y[1][::2], label="Numeric y'", linestyle='', marker='.')
plt.show()
EXACT SOLUTION: exact solution is made by integrating the beam equation twice using definite integrals and use the boundary conditions to define the integral constants. Everything is explained in the wiki link above. Below is the code to plot the y''(x), y'(x) (slope) and y(x) (deflection). The diagram is turned around, the free end of the beam is at 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()
APPROXIMATE SOLUTION (KIND OF): the code below was made by willcrack. The shape looks better than in the previous question but the values are still not ok.
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()
You are integrating a differential equation, your approach of computing in a loop the definite integrals is, let's say, sub-optimal.
The standard approach in Scipy is the use of scipy.integrate.solve_ivp, that uses a suitable integration method (by default, Runge-Kutta 45) to provide the solution in terms of a special object.
As usual in the field of numerical integration of ordinary differential equations, the method is limited to a system of 1st order differential equation, but your 2nd degree equation can be transformed to a system of 1st degree equations introducing an helper function
Y" = M ⇒ {y' = h, h' = M}
While this sounds complicated, its implementation is quite simple
In [51]: #########################################################################
...: # L, EJ = 1.0
...: #########################################################################
...: # exact solution
...: from numpy import linspace
...: x = linspace(0, 1, 51)
...: y1, y0 = (x**2-1)/2, (x**3-3*x+2)/6
...: #########################################################################
...: # numerical solution
...: from scipy.integrate import solve_ivp
...: s = solve_ivp(
...: lambda x, Y: [Y[1], x], # Y[0]=y0, Y[1]=y1, dy0dx=y1, dy1dx=x
...: [1.0, 0.0], # interval of integration (NB: reversed, because...)
...: [0.0, 0.0], # initial conditions (at the 1st point of integ interval)
...: t_eval=np.linspace(1, 0, 101) # where we want the solution to be known
...: )
...: #########################################################################
...: # plotting
...: from matplotlib.pyplot import grid, legend, plot
...: plot(x, y0, label="Exact y")
...: plot(x, y1, label="Exact y'")
...: plot(s.t[::2], s.y[0][::2], label="Numeric y", linestyle='', marker='.')
...: plot(s.t[::2], s.y[1][::2], label="Numeric y'", linestyle='', marker='.')
...: legend() ; grid() ;
In [52]:
The OP reported an issue understanding solve_ivp(lambda x, Y: [Y[1], x], ....
We have a system of 1st order ODEs in normal form
y₁' = f₁(x, y₁(x), …, yₙ(x))
… = …
yₙ' = f₁(x, y₁(x), …, yₙ(x))
that can be written, using capital letters to signify vector quantities
Y' = F(x, Y(x))
to solve the system of differential equations solve_ipv needs exactly this F(x, Y) function.
In place of the lambda expression one could write a function definition like the following, that is possibly more self explanatory
def F(x, Y):
# unpack the vector of function values
y = Y[0]
h = Y[1]
# compute the derivatives
dy_over_dx = h
dh_over_dx = x
# return the vector of derivatives values
return [dy_over_dx, dh_over_dx]
s = solve_ivp(F, …)
that in the answer was succinctly (too much succinctly?) was expressed as lambda x,Y:[Y[1],x] …
这篇关于Python - 用 scipy 求解伯努利梁方程的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
