用程序來(lái)求積分的方法有很多，這篇文章主要是有關(guān)牛頓-科特斯公式。

學(xué)過(guò)插值算法的同學(xué)最容易想到的就是用插值函數(shù)代替被積分函數(shù)來(lái)求積分，但實(shí)際上在大部分場(chǎng)景下這是行不通的。

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

牛頓-科特斯公式解決這一問(wèn)題的辦法是將大的插值區(qū)間分為一堆小的插值區(qū)間，使得多項(xiàng)式的次數(shù)不會(huì)太高。然后通過(guò)引入?yún)?shù)函數(shù)

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

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

n是一個(gè)事先確定好的值。

又因?yàn)橐粋€(gè)大的插值區(qū)間需要被分為等長(zhǎng)的多個(gè)小區(qū)間，并在這些小區(qū)間上分別進(jìn)行插值和積分，因此此時(shí)的牛頓-科特斯公式被稱(chēng)為：復(fù)化牛頓-科特斯公式。

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

當(dāng)n=2時(shí)，復(fù)化牛頓-科特斯公式被稱(chēng)為復(fù)化辛普森公式（非美國(guó)法律界著名的那個(gè)辛普森）。

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

首先寫(xiě)一個(gè)函數(shù)求節(jié)點(diǎn)函數(shù)值求和那部分：

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

然后就可以寫(xiě)整個(gè)求積分函數(shù)了：

"""
@brief: 求func積分 :
   
@param: a 積分區(qū)間左端點(diǎn)
@param: b 積分區(qū)間右端點(diǎn)
@param: n 積分分為n等份（復(fù)化梯形求積分要求）
@param: func 求積函數(shù)
@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))

相當(dāng)?shù)暮?jiǎn)單

試驗(yàn)：

當(dāng)把大區(qū)間分為兩個(gè)小區(qū)間時(shí)：

分為20個(gè)小區(qū)間時(shí)：

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

測(cè)試代碼：

if __name__ == "__main__":
  
 func = lambda x: x**2
 a, b = 2, 8
 n = 20
 print integral(a, b, n, func)
  
 ''' 畫(huà)圖 '''
 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并不是上文開(kāi)篇提到的公式中的n，開(kāi)篇提到的n是指將每一個(gè)具體的插值區(qū)間（也就是小區(qū)間）等距插n個(gè)節(jié)點(diǎn)，復(fù)化梯形公式的n是固定的為1.

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

以上這篇復(fù)化梯形求積分實(shí)例——用Python進(jìn)行數(shù)值計(jì)算就是小編分享給大家的全部?jī)?nèi)容了，希望能給大家一個(gè)參考，也希望大家多多支持腳本之家。

最新評(píng)論