前回のあらまし
前回の記事では三角関数の値を計算し、その結果をpythonのライブラリを利用して描画する練習問題を解いた。いよいよ今回は、その成果を数値積分に応用すべく、まずは数値積分の部分のコードを書き下してみる。
数値積分の台形公式
数値積分のもっとも簡単な形式は、高校の数学の教科書に載っているように
\begin{equation}\frac{x_N-x_0}{N} \sum_{i=0}^Nf(x_i) \rightarrow \int _{a=x_0} ^{b=x_N} f(x)dx \quad (N\rightarrow\infty)\label{def-integral}\end{equation}
という定義式を用いる方式である。しかし、この定義では$N\rightarrow\infty$という極限を想定している点が要注意点である。実際の数値計算では無限個の和は計算できないので、$N$を大きな数にするにせよ、有限の値が選ばれることになる。したがって、その「打ち切り誤差」をなるべく小さくするように公式を工夫する必要がある。
積分の台形公式は、離散化した関数を台形の面積で囲むことで近似する数値積分であるが、その実装は非常に簡単である。つまり、端点（始点と終点の２点のみ）の値を半値に修正するだけでよいのだ！
まずは積分 \begin{equation}\int_0^{\pi/2}\sin(x)dx = \left[-\cos(x)\right]_0^{\pi/2} = -0-(-1) = 1\end{equation} をpythonの数値積分がどの程度の精度で計算できるか確認してみることにする。
また積分のオリジナルの定義式(\ref{def-integral})と台形公式とを用いた数値結果の比較も行ってみたい。
プログラムは次のような感じに書いてみた。
#!/usr/local/bin/python3 # Integral of sin(x) # # Ver. 0 [Oct.4, 2023] import math NMAX=2000 N=10 while True: dx = 0.5*math.pi/N sum = 0.0 for i in range(N+1): x = dx * i sum += math.sin(x) # Edge-point handling for the trapezoidal formula sum -= 0.5*(math.sin(0) + math.sin(N*dx)) print(N, sum * dx) N += 10 if N > NMAX: break
まずは$N=10$から出発し、$N=10, 20, 30,\cdots$と次第に刻み幅を増やしながら、数値積分がどのように変わっていくかを調べている。プログラム後半の
# Edge-point handling for the trapezoidal formula sum -= 0.5*(math.sin(0) + math.sin(N*dx))
の部分が台形公式に基づく修正部分である。この一文をコメントアウトするとオリジナルの積分の定義式(\ref{def-integral})になる。
計算結果は一度外部ファイルに出力し、それをgnuplotで描画した。
もっとも粗い計算($N=10$)の場合でも、すでに10%未満の精度が確保されている。しかし、台形公式の精度は定義式の単なる適用に比べれば圧倒的に良いことがわかるだろう。$N=10$の場合、オリジナルの定義式では誤差8%を切る程度であるが、台形公式はすでに誤差1%を切る精度である(0.3%)。$N=50$では、オリジナルがようやく2%を切るところまで到達したのに対し、台形公式は0.01%を達成している。台形公式は端点２点のみに修正を加えるだけであることを考えれば、それを使わない手はないと誰もが思うだろう。
次は、台形公式自体の収束の程度を見てみよう。スケールに騙されないように、オリジナルの定義式の計算結果を除いて表示する。せっかくだから、matplotlibを使って描画してみた。
この計算に用いたコードは下の通りである。
#!/usr/local/bin/python3 # Integral of sin(x) # # Ver. 1 [Oct.4, 2023] import math import matplotlib.pyplot as pl NMAX=200 Nsp = 10 MMAX = NMAX // Nsp xp = [0.0] * MMAX yp = [0.0] * MMAX N=10 count = 0 while True: dx = 0.5*math.pi/N sum = 0.0 for i in range(N+1): x = dx * i sum += math.sin(x) # Edge-point handling for the trapezoidal formula sum -= 0.5*(math.sin(0) + math.sin(N*dx)) sum *= dx print(N, sum) xp[count] = N yp[count] = sum count += 1 N += Nsp if N > NMAX: break pl.plot(xp,yp,"b--", xp,yp,"ro") pl.show()
このテスト計算の結果をもとに、実際の計算においても台形公式を用いた$N=200$程度の刻み幅の数値計算を実施しようと思う。
割り算と整数
上のプログラムで、配列の次元を指定するために、整数同士の割り算をしている部分が（冒頭近くに）ある。
NMAX=200 Nsp = 10 MMAX = NMAX // Nsp xp = [0.0] * MMAX yp = [0.0] * MMAX
pythonの演算の基本は、どうやら浮動小数らしい。単にNMAX / Nspとやるとエラーが出てしまう。結果が浮動小数になっているから、配列の次元の指定方法に逸脱する、というわけである。FORTRANでもC言語でも、5 / 2 といった演算は自動的に整数型が維持されて2という結果となる。しかし、pythonはデフォルトで浮動小数が返されるため$5/2\rightarrow 2.5$となるのだそうだ（参考文献はこちら）。
pythonでは（上のプログラムリストの３行目にあるように）5 // 2と書くことで2という整数型の返り値を得ることになっているそうである。
面白いことに、たとえ割り切れる演算だとしても、4 / 2 =2.0という風に浮動小数にキャストされてしまう設定だという。したがって、どんな場合でも割り算の商を得るには4 // 2 (=2)と書かねばならない。
いよいよ$\sin (x^2)$の積分へ
これで準備が全て完了した。やっと入試問題に戻ることができる。次の記事でいよいよ \begin{equation}A_k = \int _{\sqrt{k\pi}} ^{\sqrt{(k+1)\pi}} \sin ( x^2 ) dx\end{equation} の計算を実施してみよう。
つづく
