挠度曲线微分方程的应用
有限元计算只能得到梁单元两端的节点位移和内力，单元任意截面的位移和内力需要通过插值得到。
梁的挠曲线近似微分方程为
弯矩 、剪力 和载荷集度 之间的微分关系
式(1)对 求导并结合(2)(3)得
由式(5)可以看出，如果在梁的某一区间 上载荷集度 为0，则挠曲线方程为关于 的三次多项式;如果载荷集度 为非零常数，则挠曲线方程为关于 的四次多项式;如果载荷集度 为关于 的一次多项式，则挠曲线方程为关于 的五次多项式。依此类推。而构造梁单元只使用了三次Hermite多项式，因此有限元计算的结果是近似值。
为了使得求解过程一般化，将区间 映射至标准区间 ，映射关系为
其中 ，记
由链式求导法则
则
结合(1)(4)(5)有
q为0的情况
梁在区间 不受分布力作用，即 。将区间 映射到区间 ，如图1所示。在区间 ，挠度 是关于 的三次多项式，在区间 ，挠度 是关于 的三次多项式。记区间 两个节点挠度和转角分别为 和 ， 和 ，注意到 在节点1和2处的值分别是 和 ，利用挠度和转角间的微分关系，构造挠度 的Hermite 插值多项式为
其中
令 ,便可由(8)得到中点的挠度，弯矩和剪力。
q为常数的情况
在区间 ，挠度 是关于 的四次多项式，在区间 ，挠度 是关于 的四次多项式。记区间 两个节点挠度和转角分别为 和 ， 和 ，以及中间点的挠度 ，注意到 在节点1和2处的值分别是 和 ，利用挠度和转角间的微分关系，构造挠度 的Hermite 插值多项式为
其中
q为线性分布的情况
在区间 ，挠度 是关于 的五次多项式，在区间 ，挠度 是关于 的五次多项式。记区间 两个节点挠度和转角分别为 和 ， 和 ，以及中间点的挠度 和转角 ，注意到 在节点1、2和3处的值分别是 和 ，利用挠度和转角间的微分关系，构造挠度 的Hermite 插值多项式为
其中
举例
以长度为1，端部集中力为1的悬臂梁为例，节点位移应该是 ,将梁划分为25个小段，画出弯矩图,剪力图和变形图。
Python代码
import matplotlib.pyplot as plt
import numpy as np
def main():
plt.style.use("ggplot") # 绘图风格
plt.rcParams['font.sans-serif'] = ['MiSans'] # 这个字体能够支持中英文字符
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
EI = 1
l = 1
d = np.array( [0,0,-0.333,-0.5] ) # 节点位移
xi = np.linspace(-1, 1, 25) # 梁划分为25个小段
len = np.size(xi)
# 插值基函数
N1 = 0.25 *(1-xi)**2 * (2+xi)
N2 = 0.25 *(1-xi)**2 * (1+xi)
N3 = 0.25 *(1+xi)**2 * (2-xi)
N4 = 0.25 *(1+xi)**2 * (xi-1)
# 插值基函数的二阶导数
N1_2 = 0.25 * 6 * xi
N2_2 = 0.25 * (-2+6*xi)
N3_2 = -0.25 * 6 * xi
N4_2 = 0.25 *(2+6*xi)
# 插值基函数的三阶导数
N1_3 = 0.25 * 6 *np.ones(len)
N2_3 = 0.25 * 6 *np.ones(len)
N3_3 = -0.25 * 6 *np.ones(len)
N4_3 = 0.25 * 6 *np.ones(len)
v = N1*d[0] + 0.5*l *N2*d[1] + N3*d[2] + 0.5*l *N4*d[3]
v_2 = N1_2 *d[0] + 0.5*l *N2_2 *d[1] + N3_2 *d[2] + 0.5*l *N4_2 *d[3]
v_3 = N1_3*d[0] + 0.5*l *N2_3*d[1] + N3_3*d[2] + 0.5*l *N4_3*d[3]
m = -(2/l)**2 * EI * v_2
fs = (2/l)**3 * EI * v_3
x1 = np.array( xi[0] )
x2 = np.array( xi[-1] )
xx = np.hstack( (x1, xi, x2) )
m1 = np.array( [0] )
mm = np.hstack( (m1, m, m1) )
fs1 = np.array( [0] )
fs2 = np.hstack( (fs1, fs, fs1) )
fig, (ax0,ax1,ax2) = plt.subplots(nrows=3, ncols=1, sharex=True,
figsize=(12, 6))
ax0.plot([xi[0],xi[-1]], [0,0])
ax0.plot(xx, mm)
ax0.set( title='弯矩图')
ax0.grid()
ax1.plot([xi[0],xi[-1]], [0,0])
ax1.plot(xx, fs2)
ax1.set( title='剪力图')
ax1.grid()
ax2.plot([xi[0],xi[-1]], [0,0])
ax2.plot(xi, v)
ax2.set( title='变形图')
ax2.grid()
fig.savefig("test.png")
plt.show()
if __name__ == "__main__":
main() 
以长度为1，均布荷载为8的简直梁为例，节点位移应该是 ,再将跨中线位移-0.104组成插值位移 。将梁划分为25个小段，画出弯矩图,剪力图和变形图。
import matplotlib.pyplot as plt
import numpy as np
def main():
plt.style.use("ggplot") # 绘图风格
plt.rcParams['font.sans-serif'] = ['MiSans'] # 这个字体能够支持中英文字符
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
EI = 1
l = 1
d = np.array( [0,-0.333, 0, 0.333, -0.104] ) # 插值位移
xi = np.linspace(-1, 1, 25) # 梁划分为25个小段
# 插值基函数
N1 = -0.25 * (xi-1)**2 * (2*xi**2 + 3*xi)
N2 = -0.25 *(xi-1)**2 * (xi**2 + xi)
N3 = -0.25 *(1+xi)**2 * (2*xi**2 - 3*xi)
N4 = 0.25 *(1+xi)**2 * (xi**2 - xi)
N5 = (xi-1)**2 * (1+xi)**2
# 插值基函数的二阶导数
N1_2 = -0.25 * (24*xi**2 -6*xi -8)
N2_2 = -0.25 * (12*xi**2 -6*xi -2)
N3_2 = -0.25 * (24*xi**2 +6*xi -8)
N4_2 = 0.25 * (12*xi**2 +6*xi -2)
N5_2 = 12*xi**2 - 4
# 插值基函数的三阶导数
N1_3 = -0.25 * (48*xi -6)
N2_3 = -0.25 * (24*xi -6)
N3_3 = -0.25 * (48*xi +6)
N4_3 = 0.25 * (24*xi +6)
N5_3 = 24*xi
v = N1*d[0] + 0.5*l *N2*d[1] + N3*d[2] + 0.5*l *N4*d[3] +N5*d[4]
v_2 = N1_2 *d[0] + 0.5*l *N2_2 *d[1] + N3_2 *d[2] + 0.5*l *N4_2 *d[3]+N5_2 *d[4]
v_3 = N1_3*d[0] + 0.5*l *N2_3*d[1] + N3_3*d[2] + 0.5*l *N4_3*d[3] +N5_3*d[4]
m = -(2/l)**2 * EI * v_2
fs = (2/l)**3 * EI * v_3
print(m)
x1 = np.array( xi[0] )
x2 = np.array( xi[-1] )
xx = np.hstack( (x1, xi, x2) )
m1 = np.array( [0] )
mm = np.hstack( (m1, m, m1) )
fs1 = np.array( [0] )
fs2 = np.hstack( (fs1, fs, fs1) )
fig, (ax0,ax1,ax2) = plt.subplots(nrows=3, ncols=1, sharex=True,
figsize=(12, 6))
ax0.plot([xi[0],xi[-1]], [0,0])
ax0.plot(xx, mm)
ax0.set( title='弯矩图')
ax0.grid()
ax1.plot([xi[0],xi[-1]], [0,0])
ax1.plot(xx, fs2)
ax1.set( title='剪力图')
ax1.grid()
ax2.plot([xi[0],xi[-1]], [0,0])
ax2.plot(xi, v)
ax2.set( title='变形图')
ax2.grid()
fig.savefig("test1.png",dpi=600)
plt.show()
if __name__ == "__main__":
main()
来源：数值分析与有限元编程 
