Python - Bernoulli'nin ışın denklemini scipy ile çözme

Nov 26 2020

Soruyu cevaplama süreci, aşağıdaki bağlantıdaki soruda çoktan başladı, ancak bu konu özellikle yanıtlanan bir işlevi entegre etmekle ilgiliydi. Bu yüzden yeni bir soru ekledim.

Python - Bir işlevi entegre etmek ve sonuçları çizmek

SORUN: scipy integratı kullanarak bir y '' (x) = M (x) / (E * I) kiriş denklemi nasıl çözülür .

ÇÖZÜM, gboffi'nin izniyle:

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

KESİN ÇÖZÜM: Kesin çözüm, kiriş denkleminin belirli integraller kullanılarak iki kez entegre edilmesi ve integral sabitleri tanımlamak için sınır koşullarının kullanılmasıyla yapılır. Her şey yukarıdaki wiki bağlantısında açıklanmıştır. Aşağıda, y '' (x), y '(x) (eğim) ve y (x) (sapma) değerlerini çizecek kod bulunmaktadır. Diyagram ters çevrilir, ışının serbest ucu x = 0'dadır.

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

YAKLAŞIK ÇÖZÜM (TÜR): Aşağıdaki kod willcrack tarafından yapılmıştır. Şekil, bir önceki sorudan daha iyi görünüyor, ancak değerler hala iyi değil.

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

Yanıtlar

2 gboffi Nov 26 2020 at 09:45

Diferansiyel bir denklemi entegre ediyorsunuz, hesaplama yaklaşımınız bir döngüde belirli integraller diyelim, optimalin altında.

Scipy'deki standart yaklaşım scipy.integrate.solve_ivp, çözümü özel bir nesne açısından sağlamak için uygun bir entegrasyon yöntemi (varsayılan olarak, Runge-Kutta 45) kullanan kullanılmasıdır.

Adi diferansiyel denklemlerin sayısal entegrasyonu alanında her zaman olduğu gibi, yöntem 1. mertebeden diferansiyel denklem sistemiyle sınırlıdır, ancak 2. derece denkleminiz bir yardımcı fonksiyon sunan 1. derece denklemler sistemine dönüştürülebilir.