Python - giải phương trình chùm Bernoulli với scipy

Nov 26 2020

Quá trình trả lời câu hỏi đã bắt đầu trong câu hỏi trên liên kết bên dưới, nhưng chủ đề đó đặc biệt là về tích hợp một hàm, đã được trả lời. Vì vậy, tôi đã thêm một câu hỏi mới.

Python - Tích hợp một hàm và vẽ biểu đồ kết quả

BÀI GIẢI: làm thế nào để giải một phương trình chùm y '' (x) = M (x) / (E * I) bằng tích phân scipy.

GIẢI PHÁP, lịch sự của 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()

GIẢI CHÍNH XÁC: lời giải chính xác được thực hiện bằng cách tích phân phương trình chùm hai lần sử dụng các tích phân xác định và sử dụng các điều kiện biên để xác định các hằng số tích phân. Mọi thứ được giải thích trong liên kết wiki ở trên. Dưới đây là mã để vẽ biểu đồ y '' (x), y '(x) (độ dốc) và y (x) (độ lệch). Sơ đồ quay quanh, đầu tự do của chùm tại 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()

GIẢI PHÁP PHÊ DUYỆT (LOẠI): mã bên dưới được tạo bởi willcrack. Hình dạng có vẻ tốt hơn trong câu hỏi trước nhưng các giá trị vẫn không ổn.

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()

Trả lời

2 gboffi Nov 26 2020 at 09:45

Bạn đang tích phân một phương trình vi phân, cách tiếp cận của bạn tính toán trong một vòng lặp, tích phân xác định, giả sử, là tối ưu phụ.

Cách tiếp cận tiêu chuẩn trong Scipy là sử dụng scipy.integrate.solve_ivp, sử dụng một phương pháp tích hợp phù hợp (theo mặc định, Runge-Kutta 45) để cung cấp giải pháp cho một đối tượng đặc biệt.

Như thường lệ trong lĩnh vực tích phân số của phương trình vi phân thông thường, phương pháp này được giới hạn trong hệ phương trình vi phân bậc 1, nhưng phương trình bậc 2 của bạn có thể được chuyển thành hệ phương trình bậc 1 bằng cách sử dụng hàm trợ giúp