Python - Integrating A Function And Plotting Results

I'm trying to solve Bernoulli's beam equation numerically and plotting the results. First derivative of the equation is the slope and the second derivative is the deflection. I app

Solution 1:

Maybe you can try something like this

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 equationdefd2y_dx2(x,y=None):
    return (-F*x)/(E*I)


defsomething(x):
    return integrate.quad(d2y_dx2)[0]
    
# Define the integration1 - slopedefslope(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 - deflectiondefdefl(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()

Result:

Solution 2:

I think I found the solution for the slope. I'll try the other one later. Here's the update.

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 equationdefd2y_dx2(x,y=None):
    return (-F*x)/(E*I)

    
# Define the integration1 - slopedefslope(x):
    slope_res = np.zeros_like(x)
    for i,val inenumerate(x):
        y,err = integrate.quad(f,a,val)
        slope_res[i]=y
    return slope_res

# Define the integration1 - deflectiondefdefl(x):
    
    defl_res = np.zeros_like(x)
    for i,val inenumerate(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()

New Result:

Still struggling with the last one...