Python - Integrando uma função e plotando resultados

Nov 24 2020

Estou tentando resolver a equação do feixe de Bernoulli numericamente e plotar os resultados. A primeira derivada da equação é a inclinação e a segunda derivada é a deflexão. Abordei o problema passo a passo, primeiro plote a função e depois integre-a e plote os resultados da integração no mesmo diagrama.

Meu código até agora está abaixo. Suspeito que o problema esteja no fato de que o integrate.quad retorna um único valor e estou tentando obter vários valores dele. Alguém sabe como abordar isso?

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

EDIT: O código abaixo é modificado com a resposta de willcrack Este código agora funciona, mas os resultados não estão corretos. Na parte inferior, adicionei o código para traçar os resultados usando soluções analíticas da equação do feixe que estão corretas.

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

Solução analítica, código e resultados abaixo. A forma do feixe desviado é girada, a extremidade do feixe está em 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()

Respostas

1 willcrack Nov 25 2020 at 12:02

Talvez você possa tentar algo assim

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

Resultado:

1 willcrack Nov 25 2020 at 17:13

Acho que encontrei a solução para a inclinação. Vou tentar o outro mais tarde. Aqui está a atualização.

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

Novo Resultado:

Ainda lutando com o último ...