復化梯形求積分實例——用Python進行數值計算

更新時間：2019年11月20日 10:22:21 作者：行行無別語只道早還鄉(xiāng)

今天小編就為大家分享一篇復化梯形求積分實例——用Python進行數值計算，具有很好的參考價值，希望對大家有所幫助。一起跟隨小編過來看看吧

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

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

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

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

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

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

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

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

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

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

我這篇文章實現(xiàn)的是復化梯形公式：

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

"""
@brief: 求和 ∑f(xk) : xk表示等距節(jié)點的第k個節(jié)點，不包括端點
  xk = a + kh (k = 0, 1, 2, ...)
  積分區(qū)間為[a, b]
   
@param: xk  積分區(qū)間的等分點x坐標集合（不包括端點）
@param: func 求積函數
@return: 返回值為集合的和
"""
def sum_fun_xk(xk, func):
 return sum([func(each) for each in xk])

然后就可以寫整個求積分函數了：

"""
@brief: 求func積分 :
   
@param: a 積分區(qū)間左端點
@param: b 積分區(qū)間右端點
@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))

相當的簡單

試驗：

當把大區(qū)間分為兩個小區(qū)間時：

分為20個小區(qū)間時：

求的積分值就是這些彩色的梯形面積之和。

測試代碼：

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

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

以上這篇復化梯形求積分實例——用Python進行數值計算就是小編分享給大家的全部內容了，希望能給大家一個參考，也希望大家多多支持腳本之家。

您可能感興趣的文章: