Pythonで学ぶ数値積分

数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使って数値積分の台形公式を用い、積分の値を数値的に求めてみたいと思います。Python（SciPy）に組み込まれている数値積分の方法もあわせて紹介し、その値を比較してみます。

積分

積分

関数 $y=f(x)$ を区間 $a\leq x \leq b$ で積分するというのは、 $xy$ -平面上に関数 $y=f(x)$ のグラフを描いたとき、区間 $a\leq x \leq b$ 内で、グラフと $x$ 軸とが囲む領域の（符号付きの）面積を求めることに相当します。この面積を

$\quad \displaystyle{\int_a^bf(x)dx}$

と表します。

台形公式

関数 $y=f(x)$ の積分

$\quad \displaystyle{S[a,b]=\int^b_a f(x)dx}$

を考えましょう。積分区間 $[a,b]$ を微小区間 $h$ によって $N$ 分割します：

$\quad \displaystyle{x_n=nh+a\quad(n=0,1,\dots,N)\qquad\qquad h=\frac{b-a}{N}}$

このとき、微小区間 $[x_n,\,x_{n+1}]$ がつくる台形の面積

$\quad \displaystyle{\Delta s_n=\frac{f(x_n)+f(x_{n+1})}{2}h}$

によって各微小区間からの積分への寄与を近似し、これをもって数値的な積分の近似値を与える方法を台形公式といいます。つまり、積分の近似値を

$\quad \displaystyle{\begin{aligned}S[a,b]=\int^b_a f(x)dx\ &\sim \ \sum_{n=0}^N \frac{f(x_n)+f(x_{n+1})}{2}h \\ &=h \left(\frac{f(x_0)}{2}+f(x_1)+f(x_2)+\dots+f(x_{N-1})+\frac{f(x_N)}{2}\right)\end{aligned}}$

によって与えます。これをPythonで実装してみましょう。

裳華房の『数値計算』（柳田・中木・三村）では、数値計算の基本的な手法を紹介しています：

<br />

積分を数値的にいかにして求めるかを数学的背景を基礎に解説していて、原理の面から理解したいという人にはうってつけの本だと思います。

例：１次関数

まず具体例として１次関数

$\quad \displaystyle{f(x)=x}$

の区間 $0\leq x\leq 1$ での積分を考えると、台形公式は厳密に正しい答えを与えるはずなので、実装したプログラムが数値的に正しいかどうかの動作確認として使うことができるでしょう：

import numpy as np

N = 10**6    # 分割数
xmin = 0    # 積分区間の下端
xmax = 1    # 積分区間の上端
h = (xmax - xmin)/N    # 微小区間の幅
p = np.linspace( xmin, xmax, N+1)   # 積分区間を N 等分する
f = p
S = h*(np.sum(f)-f[0]/2-f[N]/2)
S

から得られる答えは

$\quad \displaystyle{0.5}$

で、正しい答えを与えています。

例：円周率

次に具体例として

$\quad \displaystyle{f(x)=\sqrt{1-x^2}}$

を考えてみましょう。この関数を区間 $0\leq x\leq 1$ で積分すると、半径１の円の４分の１の面積に等しいはずです。というのは、原点を中心とする半径１の円は方程式

$\quad \displaystyle{x^2+y^2=1}$

によって定義されますから、この方程式を $y$ について解き、正の符号をとったものが、ここで考えている関数だからです。１次関数の場合と同様にして

import numpy as np

N = 10**6
xmin = 0
xmax = 1
h = (xmax - xmin)/N
p = np.linspace( xmin, xmax, N+1)
f = np.sqrt(1-p**2)
S = h*(np.sum(f)-f[0]/2-f[N]/2)
4*S    # 円周率と直接比較できるように４倍しています

から

$\quad \displaystyle{3.1415926524138134}$

という数値を得ます。

数学定数としての円周率

Numpyには重要な数学定数のひとつである円周率の値が組み込まれています：

import numpy as np

np.pi

から、その値は（近似値として）

$\quad \displaystyle{\pi=3.141592653589793}$

ですので、上の積分から求めた数値は、９桁ほど合致していることがわかります。分割数 $N$ を大きくすれば、それに応じて精度が上がります。

組み込み関数（数値積分）との比較

さいごにSciPyライブラリで実装されている組み込み関数としての数値積分と比較してみます。ただし、SciPyライブラリの実装では、台形公式とは異なる方法で積分の近似値を与えています。さて、SciPyでは上の関数 $y=\sqrt{1-x^2}$ の数値積分を以下のようにして求めることができます：

import numpy as np
from scipy import integrate

y = lambda x: 4*np.sqrt(1-x**2)   # 円周率と比較するために４倍しています
integrate.quad(y, 0, 1)

から、その値が

$\quad \displaystyle{3.1415926535897922\pm 3.533564552071766\times 10^{-10}}$

と計算されます。かなりの高精度で円周率と一致していることがわかります。

キーワード：Python、積分、数値計算、台形公式

pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts