Python - решение уравнения пучка Бернулли с помощью scipy
Процесс ответа на вопрос уже начался в вопросе по ссылке ниже, но эта тема была конкретно об интеграции функции, на которую был дан ответ. Итак, я добавил новый вопрос.
Python - интеграция функции и построение графика
ПРОБЛЕМА: как решить уравнение пучка y '' (x) = M (x) / (E * I), используя scipy интегрировать.
РЕШЕНИЕ, любезно предоставлено 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()
ТОЧНОЕ РЕШЕНИЕ: точное решение производится путем двукратного интегрирования уравнения балки с использованием определенных интегралов и использования граничных условий для определения интегральных констант. Все объяснено в вики-ссылке выше. Ниже приведен код для построения графиков y '' (x), y '(x) (наклон) и y (x) (отклонение). Схема развернута, свободный конец балки находится в точке 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()
ПРИБЛИЗИТЕЛЬНОЕ РЕШЕНИЕ (ВИД): приведенный ниже код был создан 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)
# 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()
Ответы
Вы интегрируете дифференциальное уравнение, ваш подход к вычислению определенных интегралов в цикле, скажем так, неоптимален.
Стандартный подход в Scipy - это использование scipy.integrate.solve_ivpподходящего метода интеграции (по умолчанию Runge-Kutta 45) для предоставления решения в виде специального объекта.
Как обычно в области численного интегрирования обыкновенных дифференциальных уравнений, метод ограничен системой дифференциальных уравнений 1-го порядка, но ваше уравнение 2-й степени может быть преобразовано в систему уравнений 1-й степени, вводя вспомогательную функцию
Y" = M ⇒ {y' = h, h' = M}
Хотя это звучит сложно, его реализация довольно проста.
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]:
OP сообщил о проблеме с пониманием solve_ivp(lambda x, Y: [Y[1], x], ....
У нас есть система ОДУ 1-го порядка в нормальной форме
y₁' = f₁(x, y₁(x), …, yₙ(x))
… = …
yₙ' = f₁(x, y₁(x), …, yₙ(x))
которые можно записать заглавными буквами для обозначения векторных величин
Y' = F(x, Y(x))
для решения системы дифференциальных уравнений solve_ipvнужна именно эта F(x, Y)функция.
Вместо лямбда-выражения можно написать определение функции, подобное приведенному ниже, которое, возможно, более понятно само по себе.
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, …)
то, что в ответе было лаконично (слишком лаконично?), было выражено как lambda x,Y:[Y[1],x]...