利用Python實(shí)現(xiàn)數(shù)值積分的方法
1. 栗子
為了加深大家的印象,首先我們來(lái)看個(gè)例子:
圖示如下:
2. 矩形計(jì)算面積
我們知道,在數(shù)學(xué)中,積分運(yùn)算表示上述曲線和x軸圍成的封閉區(qū)域的面積,為此,我們?cè)跀?shù)值預(yù)算中,來(lái)近似計(jì)算上述區(qū)域的面積,最直觀的想法就是拆成一個(gè)個(gè)小的矩形來(lái)計(jì)算對(duì)應(yīng)的面積。
2.1 左側(cè)邊長(zhǎng)計(jì)算面積
為了計(jì)算每個(gè)小矩形的面積,設(shè)計(jì)到邊長(zhǎng)高的選擇,這里我么以左側(cè)函數(shù)取值作為對(duì)應(yīng)矩形的高來(lái)計(jì)算相應(yīng)的小矩形的面積,圖示如下:
對(duì)應(yīng)的代碼如下:
import numpy as np x = np.linspace(0, 3, 1001) f = lambda x: x**3 - 4*x**2 + 4*x + 2 a = 0.5 b = 2.5 Ax = np.linspace(a, b, 101) Ay = f(Ax) def defInt_left(f, a, b, N): ? ? # left-hand point ? ? result = 0; FX = []; Xn = [] ? ? dx = abs(b - a)/N ? ? while a < b: ? ? ? ? result += f(a)*dx ? ? ? ? FX += [f(a)] ? ? ? ? Xn += [a] ? ? ? ? a += dx ? ? return result, FX, Xn, dx N = 4 I_left, FX, Xn, dx = defInt_left(f, a, b, N) print(I_left)
上述代碼中,我們將橫坐標(biāo)拆分為4小份,也就是拆分成4個(gè)小矩形,然后使用函數(shù)左側(cè)的點(diǎn)坐位小矩形的高,上述代碼的運(yùn)行結(jié)果如下:
5.25
2.2 右側(cè)邊長(zhǎng)計(jì)算面積
這里和上述原理類似,只不過(guò)每個(gè)小矩形的高采用右側(cè)邊長(zhǎng)函數(shù)取值來(lái)近似計(jì)算,圖例如下:
樣例代碼如下:
def defInt_right(f, a, b, N): ? ? # right-hand point ? ? result = 0; FX = []; Xn = [] ? ? dx = abs(b - a)/N ? ? while a < b: ? ? ? ? result += f(a + dx)*dx ? ? ? ? FX += [f(a + dx)] ? ? ? ? Xn += [a] ? ? ? ? a += dx ? ? return result, FX, Xn, dx N = 4 I_right, FX, Xn, dx = defInt_right(f, a, b, N) print(I_right)
運(yùn)行結(jié)果如下:
5.0
2.3 中值邊長(zhǎng)計(jì)算面積
看了上述兩種近似計(jì)算方式,有同學(xué)就說(shuō)有取左側(cè)點(diǎn)算出來(lái)面積大的,有取右側(cè)點(diǎn)算出來(lái)面積小的,那干脆折中一下,我們來(lái)以中值坐位矩形的高來(lái)計(jì)算對(duì)應(yīng)的面積。圖例如下:
代碼實(shí)現(xiàn)如下:
def defInt_middle(f, a, b, N): ? ? # middle point ? ? result = 0; FX = []; Xn = [] ? ? dx = abs(b - a)/N ? ? while a < b: ? ? ? ? result += f(a + dx/2)*dx ? ? ? ? FX += [f(a + dx/2)] ? ? ? ? Xn += [a] ? ? ? ? a += dx ? ? return result, FX, Xn, dx N = 4 I_mid, FX, Xn, dx = defInt_middle(f, a, b, N) print(I_mid)
運(yùn)行結(jié)果如下:
5.0625
3. 梯形計(jì)算面積
讀到這里的同學(xué)可能會(huì)思考,既然可以將封閉區(qū)域劃分成一個(gè)個(gè)的小矩形,那當(dāng)然也可以將其劃分成梯形來(lái)近似計(jì)算相應(yīng)的面積,圖例如下:
樣例代碼如下:
def defInt_trapezoid(f, a, b, N): ? ? # trapezoidal rule ? ? result = 0; FXa, FXb = [], []; Xn = [] ? ? dx = abs(b - a)/N ? ? while a < b: ? ? ? ? result += (f(a) + f(a + dx))*dx/2 ? ? ? ? FXa += [f(a)]; FXb += [f(a + dx)] ? ? ? ? Xn += [a] ? ? ? ? a += dx ? ? return result, FXa, FXb, Xn, dx N = 4 I_trap, FXa, FXb, Xn, dx = defInt_trapezoid(f, a, b, N) print(I_trap)
運(yùn)行結(jié)果如下:
5.125
4. 真值比對(duì)
最后,我們來(lái)針對(duì)不同的N來(lái)講封閉區(qū)域劃分成對(duì)應(yīng)的小份,分別針對(duì)性的計(jì)算上述四種方式的積分值,樣例代碼如下:
Nx = range(1, 11) I1, I2, I3, I4 = [], [], [], [] for Ni in Nx: ? ? i1, *_ = defInt_left(f, a, b, Ni); I1 += [i1]; ? ? i2, *_ = defInt_right(f, a, b, Ni); I2 += [i2]; ? ? i3, *_ = defInt_middle(f, a, b, Ni); I3 += [i3]; ? ? i4, *_ = defInt_trapezoid(f, a, b, Ni); I4 += [i4];
最后將其與真值進(jìn)行對(duì)比,如下:
可以看出,隨著劃分區(qū)域的增多,采用梯形計(jì)算面積方式最逼近真值。
5. 總結(jié)
本文重點(diǎn)介紹了使用不同面積劃分方法來(lái)近似計(jì)算積分取值的原理和相應(yīng)的代碼實(shí)現(xiàn),其中采用梯形計(jì)算面積的方式隨著劃分子區(qū)域數(shù)目的增加最接近真值,推薦大家使用。
到此這篇關(guān)于利用Python實(shí)現(xiàn)數(shù)值積分的方法的文章就介紹到這了,更多相關(guān)Python實(shí)現(xiàn)數(shù)值積分內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家!
- python計(jì)算分段函數(shù)值的方法
- Python中的四種交換數(shù)值的方法解析
- python中利用numpy.array()實(shí)現(xiàn)倆個(gè)數(shù)值列表的對(duì)應(yīng)相加方法
- Python產(chǎn)生一個(gè)數(shù)值范圍內(nèi)的不重復(fù)的隨機(jī)數(shù)的實(shí)現(xiàn)方法
- python畫(huà)圖——實(shí)現(xiàn)在圖上標(biāo)注上具體數(shù)值的方法
- python畫(huà)柱狀圖--不同顏色并顯示數(shù)值的方法
- Python數(shù)值求解微分方程方法(歐拉法,隱式歐拉)