数値計算での積分方法、特に等間隔の分点の場合であるニュートン・コーツ積分(とロンバーグ積分)に関する理論とのプログラムを載せます。

1. ニュートンコーツ型の公式
   • 台形則
   • シンプソン積分
   • 高次≠高精度
2. ルンゲ現象
3. ロンバーグ積分
4. サブルーチン”integral”のコードと使い方
   • サブルーチン”integral”を用いた例題
5. サブルーチン”romberg”のコードと使い方
   • サブルーチン”romberg”を用いた例題

等間隔の分点で積分を行う場合、よく使われる方法は台形則、もしくはシンプソン積分です。

台形則

台形則は、分点間を台形近似して面積を求める方法であり、下図のようなイメージで積分を近似する方法です。

積分\(\displaystyle \int_a^b f(x) dx\)を、
\(
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N-1} \frac{f(x_{i+1})+f(x_i)}{2}h
\)
として近似します。この時、誤差のオーダーは\(O(h^3)\)となります。
言い換えると、
台形積分とは、隣り合う分点間を1次関数で近似して求める積分
言えます。

シンプソン積分

シンプソン積分は、隣り合う分点間を2次関数で近似して求める積分です。イメージは

な感じです。数式では
\(
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N/2} \frac{f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})}{3}h
\)
となります。

高次≠高精度

これ以降もあります。分点間を3次、4次、…としていけば高次の積分をすることが可能になります。
この分点間の次数をどんどん上げていって積分をする方法をまとめてニュートンコーツの公式、と呼びます。
例えばの台形則はニュートンコーツの公式の1次に相当し、シンプソン則はニュートンコーツの公式の2次に相当しています。

じゃあ100次の公式作れば精度が凄い公式ができるんじゃないか？となりますが、これは違います。
高次であることと精度が高いことは違うのです。だから使うな、というわけではなくて知っておいてほしいことなんです。
ではなぜか？これはwikipediaのニュートンコーツの公式のページからとってきたものですが、この式に理由があります。

問題となるのが誤差項です。例えば台形近似を見ますと、誤差項は\(h^3 f^{(2)}(\xi)\)とあります。hは刻み幅なので、この項は刻み幅を小さくすれば小さくするほど精度がよくなることを言っています。
が、問題は\(f^{(2)}(\xi)\)の方です。この項は本当の関数が分からないと評価できません。
多項式だったら微分するごとに次数が1つづつ減少するのでこの項は高次になれば小さくなります。じゃあもしも微分するごとに値が増加していったら…？この場合、hを小さくしてもそれを上回るくらい\(f^{(2)}(\xi)\)が大きくなることがあります。これをルンゲ現象と呼びます。

ルンゲ現象

例えば関数
\(
\displaystyle f(x)=\frac{1}{1+25x^2}
\)
を考えます。wikipediaのページを参考にすれば、その導関数と値は、

となります。高次になればなるほど値が増加する関数なのです。
実際には\(h^n f^{(m)}(\xi)\)の積の値が小さいか大きいか？なので、一概には高次は悪い！とは言えません。高次は危ないくらいでしょう。

これを可能な限り回避するために、低次の積分を組み合わせて使う方法がとられます。
こうすることで高次の値を増加を抑えるのです。

そのほかの手法もいろいろあります。分点の位置に縛られなければガウス求積法が最高でしょう。素晴らしい方法です。ぜひ調べてみてください。

[adsense1]

ロンバーグ積分

さて、高次を追い求めるのもいいですが、一線を駕す…とまではいきませんが、ロンバーグ(Romberg)積分というものがあります。
このロンバーグ積分とは、台形則による計算を基本とし、無限小区間での積分結果を補外によって求めよう、というものです。
ここでいう補外の概念は以下のような感じです。

刻み幅1で台形則の結果から積分値S1が分かる
↓
刻み幅0.5で台形則による結果から積分値S2がわかる
↓
刻み幅0.25で台形則による結果から積分値S3がわかる
↓
刻み幅0.175で台形則による結果から積分値S4がわかる

とある程度計算していきます。そうすると刻み幅hを変化させるとそれに応じて積分値がS1→S2→S3→S4…と推移していくわけです。
この推移の仕方を計算すれば、刻み幅無限小の場合の積分値Sが分かる、こういう仕組みになっています。収束加速法と言ったり、リチャードソン補外、なんて言ったりします。
この補外による方法はかなり優れています。先ほどの高次がどうこう、といった問題が発生しないのです。台形則しか使ってないんですから。有限桁計算に適した方法である、とかwikipediaに載ってたりします。

ロンバーグ積分の詳しい計算方法や理論は
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ

日本語では第6章 数値積分と数値微分にロンバーグ積分の記述があります。
英語ではRomberg integralで調べると良いでしょう。参考までに、英語では
Romberg integration exampleや、Romberg Integrationなにかがいいと思います。

2016/03/01
上のRomberg Integrationを参考にして書いたプログラムはこちらです。これは、解析関数fをRomberg積分するプログラムです。
解析関数ではなく、サイズ\(2^{n}+1\)の配列に格納されている場合はこのページの下の方にあるプログラムを参照してください。

で呼び出し、区間a~bの積分\(s=\int_a^b f(x) dx\)を精度preで求めるプログラムです。

サブルーチン”integral”のコードと使い方

下のサブルーチン”integral”は1,2,3次元の短冊近似、台形則、シンプソン積分、シンプソン3/8積分、ブール則(ニュートン・コーツ型積分)をカバーします。
ただし、分点の個数に注意してください。分点の個数を気にせず使える近似方法は、短冊近似、台形近似のみです。

ロンバーグ積分では引数の関係上、ニュートン・コーツ型積分との間をうまく処理できなかったため、別のサブルーチンとして書きます。もしもロンバーグ積分を使いたい場合はもう少し先に進んでください。

まずはニュートン・コーツ型積分です。
下のサブルーチンを使ってください。
使い方は、

という感じです。
実際のプログラムでは、以下のように使用してください。配列yは複素数解列(ただし実軸上の積分)でもokです。

1次元

2次元配列の場合は以下のように指定してください。

3次元の場合はこう指定してください。

積分のためのモジュール

サブルーチン”integral”を用いた例題

実際の使い方。参考にどうぞ。
下のコードは2次元の積分
\(
\displaystyle \int\int xe^{-x^2}e^{-y}dxdy=\frac{e^{-x^2}}{2}e^{-y}+\mbox{const}
\)
を数値積分します。積分範囲を\(x=0\sim 3, \; y=0 \sim 5\)にした場合、その解析解は
\(
\begin{align}
\int_0^3 dx\int_0^5 dy xe^{-x^2}e^{-y}&=\frac{1}{2}(1-e^{-5}-e^{-9}+e^{-14})\\
&=0.4965697373627734784608751894356320535936864348993604\cdots
\end{align}
\)
となります。実行すると

となります。コードは、

3次元の計算の例題です。
下のコードは
\(
\displaystyle \int\int\int x^4y^4z^4 dxdydz=\left(\frac{1}{5}\{0.7^5-(-1)^5\}\right)^3
\)
を数値積分します。積分範囲を\(x,y,z=-1\sim 0.7\)にした場合、その値は
\(
\displaystyle \int_{-1}^{0.7}\int_{-1}^{0.7}\int_{-1}^{0.7} x^4y^4z^4 dxdydz = 0.127496010896795\cdots
\)
となります。実行すると、

と得られます。中身はこうです。分点の数はx,y,z軸各々で違っていて構いません。

ロンバーグ積分のプログラム

ロンバーグ積分(romberg_mod)の例

必要なのは、上で紹介したモジュール “romberg_mod” と下のメインプログラムです。
以下のプログラムは1次元の定積分
\(
\int_1^10 \frac{1}{x^2} dx = 0.9
\)
を分点数\(2^8+1(jx=8)\)個でロンバーグ積分するものです。
精度は\(O(h^{2\cdot min{jm,jx}})\)です。
ここで\(jm\)はモジュールromberg_modの中にパラメータとして宣言されています。

コンパイルは

とでもすればいいでしょう。

実行例(要求精度は\(10^{-8}\))

2次元ではこう。

実行例

[adsense2]

Tweet