Python复化梯形公式数值积分实现与解析
本文详解Python中复化梯形公式的实现原理与代码。针对高次插值龙格现象，通过将大区间分割为等长小区间降低多项式次数。文章提供求和函数与积分主函数的完整Python代码，解析变量n在公式与代码中的区别，并通过分段绘图直观展示积分值等于梯形面积之和的计算过程。
数值积分方法概述
利用程序进行积分计算的方法众多，本文重点探讨基于牛顿-科特斯公式的数值计算方法。
熟悉插值算法的读者可能会首先想到使用插值函数替代被积函数来求解积分，然而在实际应用场景中，这种方法往往难以奏效。
插值函数通常表现为次数不超过n的多项式。若直接利用插值函数进行积分，会引入高次多项式积分问题，从而将原问题转化为求解n次多项式积分的新难题。随着多项式次数升高，易出现龙格现象（Runge phenomenon），导致误差反而增大，且高次插值求积公式可能变得不稳定，具体原因此处不再赘述。
复化牛顿-科特斯公式原理
牛顿-科特斯公式通过将该问题划分为多个较小的插值区间，从而限制多项式的最高次数，以此解决上述稳定性问题。随后引入参数函数进行推导。
通过将被积函数中带有幂次的项取值范围限制在固定区间内，可将多项式含幂部分的积分转化为固定常数。该常数可通过手工计算得出，并直接代入多项式求积公式中使用。
上式中x的积分区间为[a, b]，步长h = (b - a)/n，此时积分区间转化为[0, n]。需注意的是，该公式表明一个大区间被划分为n个等长的小区间。具体细节可参考任意数值计算相关教材。
此处的n为预先确定的固定值。
由于大插值区间需被划分为多个等长小区间，并在各小区间上分别执行插值与积分操作，因此该形式的牛顿-科特斯公式被称为复化牛顿-科特斯公式。
针对不同的n取值，复化牛顿-科特斯公式有不同的名称。当n=1时，称为复化梯形公式。该公式将每个小区间视为一个梯形（高为h，上底为f(t)，下底为f(t+1)），这与积分无限分割的本质相一致。
当n=2时，复化牛顿-科特斯公式被称为复化辛普森公式（注意：非美国法律界著名的辛普森案）。
Python代码实现
本文主要实现复化梯形公式的代码逻辑。
首先编写一个函数，用于计算节点函数值的求和部分。
"""
@brief: 求和 ∑f(xk) : xk表示等距节点的第k个节点，不包括端点
xk = a + kh (k = 0, 1, 2, ...)
积分区间为[a, b]
@param: xk 积分区间的等分点x坐标集合（不包括端点）
@param: func 求积函数
@return: 返回值为集合的和
"""
def sum_fun_xk(xk, func):
return sum([func(each) for each in xk])随后构建完整的积分计算函数。
"""
@brief: 求func积分 :
@param: a 积分区间左端点
@param: b 积分区间右端点
@param: n 积分分为n等份（复化梯形求积分要求）
@param: func 求积函数
@return: 积分值
"""
def integral(a, b, n, func):
h = (b - a)/float(n)
xk = [a + i*h for i in range(1, n)]
return h/2 * (func(a) + 2 * sum_fun_xk(xk, func) + func(b))整体实现逻辑相当简洁。
试验结果与可视化
接下来进行试验验证。
当大区间被划分为两个小区间时的计算结果。
当大区间被划分为20个小区间时的计算结果。
最终求得的积分值即为图中所有彩色梯形面积的总和。
测试代码与参数说明
以下是完整的测试代码。
if __name__ == "__main__":
func = lambda x: x**2
a, b = 2, 8
n = 20
print integral(a, b, n, func)
''' 画图 '''
import matplotlib.pyplot as plt
plt.figure("play")
ax1 = plt.subplot(111)
plt.sca(ax1)
tmpx = [2 + float(8-2) /50 * each for each in range(50+1)]
plt.plot(tmpx, [func(each) for each in tmpx], linestyle = '-', color='black')
for rang in range(n):
tmpx = [a + float(8-2)/n * rang, a + float(8-2)/n * rang, a + float(8-2)/n * (rang+1), a + float(8-2)/n * (rang+1)]
tmpy = [0, func(tmpx[1]), func(tmpx[2]), 0]
c = ['r', 'y', 'b', 'g']
plt.fill(tmpx, tmpy, color=c[rang%4])
plt.grid(True)
plt.show()需注意，上述代码中的n并非前文公式定义中的n。前文公式中的n指每个具体插值区间（即小区间）内等距插值的节点数，而在复化梯形公式中该值固定为1。
代码中的n指的是将大区间划分为n个小区间的数量。
以上即为复化梯形求积分实例的Python数值计算全部内容，希望能提供参考，欢迎支持。
本站声明
本文内容仅供参考，不代表本站观点。转载或引用请注明出处。如发现内容问题，请通过联系我们页面反馈。
