Python: resolviendo la ecuación del haz de Bernoulli con scipy

Nov 26 2020

El proceso de respuesta a la pregunta ya comenzó en la pregunta del enlace a continuación, pero ese tema era específicamente sobre la integración de una función, que fue respondida. Entonces agregué una nueva pregunta.

Python: integración de una función y trazado de resultados

EL PROBLEMA: cómo resolver una ecuación de viga y '' (x) = M (x) / (E * I) usando scipy integration.

SOLUCIÓN, cortesía de 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()

SOLUCIÓN EXACTA: la solución exacta se obtiene integrando la ecuación de la viga dos veces usando integrales definidas y usando las condiciones de contorno para definir las constantes integrales. Todo se explica en el enlace wiki de arriba. A continuación se muestra el código para trazar la y '' (x), y '(x) (pendiente) y y (x) (deflexión). El diagrama se da la vuelta, el extremo libre de la viga está en 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()

SOLUCIÓN APROXIMADA (TIPO DE): el código siguiente fue creado por willcrack. La forma se ve mejor que en la pregunta anterior, pero los valores aún no están bien.

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

Respuestas

2 gboffi Nov 26 2020 at 09:45

Está integrando una ecuación diferencial, su enfoque de calcular en un bucle las integrales definidas es, digamos, subóptimo.

El enfoque estándar en Scipy es el uso de scipy.integrate.solve_ivp, que utiliza un método de integración adecuado (por defecto, Runge-Kutta 45) para proporcionar la solución en términos de un objeto especial.

Como es habitual en el campo de la integración numérica de ecuaciones diferenciales ordinarias, el método está limitado a un sistema de ecuaciones diferenciales de primer orden, pero su ecuación de segundo grado se puede transformar en un sistema de ecuaciones de primer grado introduciendo una función auxiliar.