1. 程式人生 > 程式設計 >復化梯形求積分例項——用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,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個小區間。

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