python梯形公式_復化梯形求積分——用Python進行數值計算

用程序來求積分的方法有很多，這篇文章主要是有關牛頓-科特斯公式。

學過插值算法的同學最容易想到的就是用插值函數代替被積分函數來求積分，但實際上在大部分場景下這是行不通的。

插值函數一般是一個不超過n次的多項式，如果用插值函數來求積分的話，就會引進高次多項式求積分的問題。這樣會將原來的求積分問題帶到另一個求積分問題：如何求n次多項式的積分，而且當次數變高時，會出現龍悲歌現象，誤差反而可能會增大，並且高次的插值求積公式有可能會變得不穩定：詳細原因不贅述。

牛頓-科特斯公式解決這一問題的辦法是將大的插值區間分為一堆小的插值區間，使得多項式的次數不會太高。然后通過引入參數函數

將帶有冪的項的取值范圍固定在一個固定范圍內，這樣一來就將多項式帶有冪的部分的求積變為一個固定的常數，只需手工算出來即可。這個常數可以直接帶入多項式求積函數。

上式中x的求積分區間為[a, b]，h = (b - a)/n, 這樣一來積分區間變為[0， n]，需要注意的是從這個公式可以看出一個大的區間被分為n個等長的小區間。 這一部分具體請參見任意一本有關數值計算的書！

n是一個事先確定好的值。

又因為一個大的插值區間需要被分為等長的多個小區間，並在這些小區間上分別進行插值和積分，因此此時的牛頓-科特斯公式被稱為：復化牛頓-科特斯公式。

並且對於n的不同取值牛頓-科特斯有不同的名稱： 當n=1時，叫做復化梯形公式，復化梯形公式也就是將每一個小區間都看為一個梯形(高為h，上底為f(t), 下底為f(t+1))。這與積分的本質：無限分隔  相同。

當n=2時，復化牛頓-科特斯公式被稱為復化辛普森公式(非美國法律界著名的那個辛普森)。

我這篇文章實現的是復化梯形公式：

首先寫一個函數求節點函數值求和那部分：

"""

@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是指將每一個具體的插值區間(也就是小區間)等距插n個節點，復化梯形公式的n是固定的為1.

而代碼中的n指將大區間分為n個小區間。