Python-関数の統合と結果のプロット

Question

ベルヌーイのビーム方程式を数値的に解き、結果をプロットしようとしています。方程式の一次導関数は勾配であり、二次導関数はたわみです。私は問題に段階的に取り組み、最初に関数をプロットしてからそれを積分し、積分結果を同じ図にプロットしました。

これまでの私のコードは次のとおりです。問題は、integrate.quadが単一の値を返し、そこから複数の値を取得しようとしているという事実にあると思います。誰かがそれにアプローチする方法を知っていますか？

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

編集：ベローはwillcrackの答えで修正されたコードです。このコードは機能するようになりましたが、結果は正しくありません。下部に、正しいビーム方程式の解析解を使用して結果をプロットするためのコードを追加しました。

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

以下の分析ソリューション、コード、および結果。偏向されたビームの形状は向きを変え、ビームの端は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

多分あなたはこのようなことを試すことができます

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

結果：

Answer 2

斜面の解決策を見つけたと思います。もう1つは後で試してみます。これがアップデートです。

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

新しい結果：

まだ最後のものと格闘しています...