皆さん、こんにちは。
今回は「数値解析」の第2弾として、定積分をコンピュータに計算させる方法「数値積分」について紹介したいと思います。
1. 区分求積法
前回説明した通り、コンピュータでは積分のような極限操作の式を直接計算することができません。なので積分計算も上手く有限の四則演算に置き換えないといけません。
実は、積分を数値解析する最も原始的な方法は、高校数学の時点で登場しています。「区分求積法」です。
a≦x≦bの範囲で関数f(x)を積分すると、曲線y=f(x)と直線x=a, x=bで囲まれる部分の面積が求まりますが、区分求積法ではその面積を、横の長さを分割して長方形に分けることで近似したのでした。
積分区間をN等分して、左から順にx0, x1, ・・・,xNとすると、
区分求積法では次のように定積分の値Iを近似できます。
この分割数NをN→∞にすれば右辺がIに収束します。
左辺の積分区間をよく見る0~1に変換したければ、次のような変数変換をすればOKです。
この変換をすれば、
高校数学でよく見慣れた区分求積法の式になります。
大学入試ではΣの形を積分に直す形で使うことが殆どですが、数値解析では逆に積分をΣの形に直すことで計算をするわけです。
このようにアイデアとしては分かりやすい区分求積法ですが、実際の曲線y=f(x)の形状と長方形の間には大きな隙間やはみ出しがあり、分割数Nを相当大きくしなければ大きな誤差が生じてしまいます。
誤差を小さくするには、少しでもこの隙間やはみ出しを小さくしたいですね。
そうすべく、長方形の代わりに台形を使えないか？
その発想で生まれたのが、次に紹介する「台形公式」です。
2. 台形公式
基本的な発想は区分求積法と同じですが、区分求積法では長方形の面積を足していったところを、代わりに台形の面積を足していくのが「台形公式」です。
長方形を考えたときに比べて、明らかに隙間やはみ出しが減ってることが分かるかと思います。
xk≦x≦xk+1の部分の台形の面積skは、（上底＋下底）×高さ÷2の公式により
と書けるので、それを全部足し上げたものが定積分Iの近似値となります。
例として、定積分
の値を台形公式で求めてみましょう。
この積分はx=tanθと置換することによりπ/4になることが分かります。
台数公式に基づきこの定積分を計算するPythonコードを上げておきます。
厳密解がπ/4なので、積分結果の4倍を出力するように設定してあります。なので、このプログラムを実行した結果が円周率π=3.141592653・・・に近ければ成功というわけです。
分割数を10にして実行すると、
3.330423726463721
となり、整数部分しか一致せずまだまだπとのズレが大きいです。
ところが、分割数を100にして実行すると、
3.1415759869231303
となり、小数点以下4桁まで一致するようになりました。
さらに分割数を1000に増やすと、
3.141592486923127
となり、小数点以下6桁まで一致します。
台形公式⑧の近似精度を見ておきます。
今、定積分IをN分割してできる面積をTkとします。
f(x)の原始関数をF(x)としてTkをx=xkの周りでテイラー展開すると、
とかけます。ここでOはランダウの記号です。
一方、f(x)自体もx=xkの周りでテイラー展開すると、
x=xkにおける微分係数を書き換えることができます。
この結果を⑩に代入すれば、
Tkとskの差が、分割幅Δxの3次式以上だけの多項式になることが分かります。数値解析の分野では、この事実を「3次のオーダーの誤差」と表現します。
⑫でk=0,1,・・・,N-1として和を取れば、
N個分誤差が積み重なることで、誤差の次数が一つ減り2次のオーダーになります。
よって、台形公式による近似計算の誤差は2次のオーダーとなります。
Δxが非常に小さい数であることを考慮すると、誤差の次数のオーダーが高ければ高いほど誤差の大きさが小さい、つまり精度が高いことになります。
⑩の右辺を、2次の項まで切り捨てれば区分求積法の式となり、
区分求積法による近似計算の誤差は1次のオーダーだと分かります。
誤差のオーダーを比較すれば、やはり台形公式による計算の方が精度が高いことが分かりますね。
3. シンプソンの公式
台形公式でも誤差が２次のオーダーという比較的高い精度で定積分を計算できたのですが、それを上回る誤差を4次のオーダーにまでできるアルゴリズムが存在します。
それが、「シンプソンの公式」を用いたものです。
N分割した図形の上側を、定数関数で近似するのが区分求積法、１次式（直線）で近似するのが台形公式でした。
シンプソンの公式では、これを2次式(放物線)で近似してさらなる精度向上を図っています。
さて、大元になるシンプソンの公式を紹介します。
定積分が特定の関数の値の線形和で計算できてしまうという驚きの公式です。
しかも単なる近似ではなく、f(x)が3次式以下なら完全に一致してしまうという。
早速証明していきましょう。
まずは、f(x)が3次以下の多項式のときに完全一致することを示します。これは左辺と右辺を愚直に計算して一致することを確かめればOKです。
次に、f(x)が一般の関数の場合に右辺の形が近似値になっていることを示します。
ここでは、f(x)を、y=f(x)上の異なる3点（両端2点x=aとx=bと、その間にある１点x=c）を通る2次式g(x)で近似します。
g(x)の係数p,q,rは、
㉑～㉓の3つからなる連立１次方程式の解になります。これらを行列の形で書けば、
となります。
今回は、f(x)が2次式g(x)で近似できる、つまりp,q,rが一意に定まることが言えれば十分で、p,q,rの具体的な値の情報は必要ありません。
p,q,rが一意に定まるには、㉔の係数行列が逆行列を持てばいいですね。
連立一次方程式を解き明かそう ～逆行列・ガウスの消去法・クラメールの公式～ - ちょぴん先生の数学部屋 (hatenablog.com)
というわけで行列式を計算すると、
綺麗に因数分解された結果となり、0にならないことが分かります。
よって、㉔を満たすp,q,rがa,b,cを使って一意に定まることが分かり、そんなp,q,rを用いて
と近似できることが分かりました。
㉖で各辺を積分すれば、g(x)は2次式、つまり3次以下の多項式なので、先に証明した通り完全一致するシンプソンの公式が使えます。
ここで、x=a,b,cでy=f(x)とy=g(x)は交わるので
が言えています。よって
としてしまえば、
となって証明完了です。
さて、ここでa=xk、b=xk+1とすれば、先ほど⑨で定義したTkが、シンプソンの公式から
と近似できることになります。
この式で和を計算すれば、定積分Iの近似式が求まります。
シンプソンの公式の近似精度を確かめておきましょう。
㉛の右辺をx=xkの周りでテイラー展開すると、各項が
となるので、
となります。これは4次の項までTkのテイラー展開と一致しているので、
TkとSkの差は5次のオーダーになります。あとはこれをN回分積み重ねればいいので、
結局、シンプソンの公式を用いた近似計算での誤差は4次のオーダーとなります。
これは誤差が２次のオーダーだった台形公式と比較しても圧倒的な高精度ですね。
台形公式のときと同様、
をシンプソンの公式で近似計算してみましょう。
Pythonコードはこちらです。
分割数を10とすると、
3.331925033978532
流石にまだ整数部分しか合いません。
ところが分割数を100とすると、
3.141592653589792
この時点で小数点以下14桁まで再現できました!!
台形公式では同じ分割数100では小数点以下4桁までしか再現できてなかったことと比べると、圧倒的な高精度です。
このようにシンプルなアルゴリズムなのに高精度で定積分を計算できるのが、シンプソンの公式の特徴です。
また、積分の数値解析には「ガウス・ルジャンドルの積分公式」という手法がありますが、これは後日「有限要素法」について紹介する段で詳しく説明したいと思います。
4. 広義積分の場合
さて、ここまで考えた定積分は、積分区間があくまで有限の範囲でした。
では、積分区間に∞が入るようないわゆる広義積分を計算したい場合はどうするか？
一つは、∞などと欲張らずに1000とか10000とか十分に大きい有限の値で∞を置き換えてしまうことです。被積分関数の0への収束スピードが速ければこれでも十分よい近似になります。
もう一つは、置換積分を駆使して有限の積分区間に変換することです。
これには様々な置換の方法が考案されており、例えば次のような例があります。
今回は例として、㊵を採用してf(x)=e^(-x^2)としたもの、つまり「ガウス積分」を計算してみたいと思います。
ガウス積分 ～統計学で最も重要な積分～ - ちょぴん先生の数学部屋 (hatenablog.com)
この形に対してシンプソンの公式を使った近似計算を行います。以下がそのPythonコードです。
ガウス積分の厳密解は√πなので、積分結果の２乗を出力するように設定しました。よって、今回の場合も出力結果がπに近ければ成功です。
分割数100でこのプログラムを実行すると、
3.1415926535897163
小数点以下13桁まで一致しました！
確かにこの方法でガウス積分が近似計算できていることが分かります。
次回は、微分方程式の数値解析について紹介します。
