Python - integracja funkcji i wykreślanie wyników

Question

Próbuję numerycznie rozwiązać równanie belki Bernoulliego i wykreślić wyniki. Pierwsza pochodna równania to nachylenie, a druga to ugięcie. Podchodziłem do problemu krok po kroku, najpierw wykreślam funkcję, a następnie integrujemy ją i wykreślam wyniki całkowania na tym samym diagramie.

Jak dotąd mój kod jest poniżej. Podejrzewam, że problem polega na tym, że funkcja integrate.quad zwraca jedną wartość i próbuję uzyskać z niej wiele wartości. Czy ktoś wie, jak do tego podejść?

from scipy import integrate
import numpy as np

from pylab import *

# Beam parameters
L = 100
w = 10
h = 10
I = (w*h**3)/12
E = 200000
F = 100

def d2y_dx2(x):
    return (-F*x)/(E*I)

a = 0.0
b = L

res, err = integrate.quad(d2y_dx2, a, b)

t = np.linspace(a,b,100)

ax = subplot(111)
ax.plot(t, d2y_dx2(t))
ax.plot(t, res(t))

show()

EDYCJA: Poniżej jest zmodyfikowany kod z odpowiedzią willcracka. Ten kod teraz działa, ale wyniki są nieprawidłowe. Na dole dodałem kod do wykreślania wyników za pomocą analitycznych rozwiązań równania belki, które są poprawne.

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)


def something(x):
    return integrate.quad(d2y_dx2)[0]
    
# Define the integration1 - slope
def slope(t):
    slope_res = []
    for x in t:
        res1, err = integrate.quad(d2y_dx2, a, b)
        slope_res.append(res1)
    return slope_res

# Define the integration1 - deflection
def defl(t1):
    defl_res = []
    for t in t1:
        res2, err = integrate.dblquad(d2y_dx2,a,b, lambda x: a, lambda x: b)
        defl_res.append(res2)
    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()

Poniżej przedstawiono rozwiązanie analityczne, kod i wyniki. Kształt odchylonej belki jest odwrócony, koniec belki znajduje się w punkcie 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()

Answer 1

Może możesz spróbować czegoś takiego

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)


def something(x):
    return integrate.quad(d2y_dx2)[0]
    
# Define the integration1 - slope
def slope(t):
    slope_res = []
    for x in t:
        res1, err = integrate.quad(d2y_dx2, a, b)
        slope_res.append(res1)
    return slope_res

# Define the integration1 - deflection
def defl(t1):
    defl_res = []
    for t in t1:
        res2, err = integrate.dblquad(d2y_dx2,a,b, lambda x: a, lambda x: b)
        defl_res.append(res2)
    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()

Wynik:

Answer 2

Myślę, że znalazłem rozwiązanie na stok. Spróbuję później. Oto aktualizacja.

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

Nowy wynik:

Wciąż zmagam się z ostatnim ...