Python - การแก้สมการลำแสงของ Bernoulli ด้วย scipy

Question

ขั้นตอนการตอบคำถามได้เริ่มต้นแล้วในคำถามบนลิงค์ร้อง แต่หัวข้อนั้นเกี่ยวกับการรวมฟังก์ชันซึ่งได้รับคำตอบโดยเฉพาะ ผมจึงเพิ่มคำถามใหม่

Python - การรวมฟังก์ชันและการพล็อตผลลัพธ์

ปัญหา: วิธีแก้สมการคาน y '' (x) = M (x) / (E * I) โดยใช้ scipy integrate

โซลูชันได้รับความอนุเคราะห์จาก 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()

Answer 1

คุณกำลังรวมสมการเชิงอนุพันธ์แนวทางการคำนวณของคุณในลูปอินทิกรัลที่แน่นอนคือสมมุติว่าย่อยที่เหมาะสมที่สุด

แนวทางมาตรฐานใน 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]: 

---

---

solve_ivp(lambda x, Y: [Y[1], x], ...สหกรณ์รายงานปัญหาการทำความเข้าใจ

เรามีระบบ ODE ลำดับที่ 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]...