数値積分の種類

ある関数の数値積分を考えます。
解析的に積分を行う事をまず考えましょう。1回の関数の評価で、全桁一致するからです。

  1. 解析的に計算
  2. どんな値でも関数値が得られる場合
  3. 特定の位置の値だけ関数値が得られる場合

解析的に計算


例えば,すぐには思いつかなそうな積分
\(
\displaystyle \int \sqrt{x}e^{-x^2} dx
\)

http://www.wolframalpha.com/input/?i=integral+sqrt(x)*exp(-x%5E2)
でも、wolfram alphaで計算させますと、不完全ガンマ関数を用いることで計算することが出来ることが分かります。

また、
森口繁一, 宇田川銈久, 一松信著『数学公式Ⅰ』、『数学公式Ⅱ』、『数学公式Ⅲ』岩波書店(1960)
にwolframでも答えてくれないような一般的な関数の積分も載っていますので、こちらも参照しましょう。

残念ながら解析的に積分出来ない場合、数値積分を行います。
ほとんどの数値積分法は、関数\(f(x)\)の積分を
\(
\displaystyle \int f(x) dx \approx \sum_{i} w_i f(x_i)
\)

と右辺の形で近似します。
ここで右辺に現れる、
あらかじめ決められた分点の位置\(x_i\)と
あらかじめ決められた重み\(w_i\)
をどのように決めるのか?で数値積分の良さが決まります。
上の例にあてはまらないものとして乱数を利用するモンテカルロ積分が挙げられます。

数値積分は大きく2つに分けられます。

1. どんな位置の値でも関数値が得られる場合
2. 特定の位置の値だけ関数値が得られる場合

です。

どんな位置の値でも関数値が得られる場合


とても簡単に実装、1-2桁程度合っていれば良い場合


モンテカルロ積分が簡単で、被積分関数が非連続性を持っていても、積分が存在すれば求めることが出来ます。
欠点としては、再現性が乱数に依存してしまう点でしょう。シード値を記録して置けばよいですが…。
良い乱数の発生方法として、メルセンヌツイスタと呼ばれる乱数の発生方法が良いそうです。

一様乱数を用いて1次元積分を行う場合、
\(
\displaystyle \int_a^b f(x)dx \approx \frac{b-a}{N}\sum_{i=1}^N f(x_i)+O(\frac{1}{\sqrt{N}})
\)

として計算できます。ここで\(x_i\)は位置\([a,b]\)に渡って発生される一様乱数です。
コードでは積分値をsとすると

s=0d0
do i=1,N
   x=([0,1)の範囲で発生する乱数)*(b-a)+a
   s=s+f(x)
enddo
s=s*(b-a)/dble(N)

で実装できます。

多重積分の場合、通常の積分法の収束が指数関数的に遅くなるので、モンテカルロ積分は良く使われます。
詳細はモンテカルロ積分をご覧ください。

モンテカルロ法 (Monte Carlo method)
確率密度関数からモンテカルロ積分まで

簡単に実装、2-4桁程度合っていれば良い場合


台形則による積分が良いでしょう。
区間\([a,b]\)で積分したい場合、区間を\(N+1\)点で分割すると
・等間隔の分点の場合
\(
\displaystyle \int_a^b f(x) dx \approx \frac{b-a}{N}\left(\frac{1}{2}[f(a)+f(b)]+\sum_{i=1}^{N-1} f(i\frac{b-a}{N}+a)\right)
\)

・非等間隔の分点の場合
\(
\displaystyle \int_a^b f(x) dx \approx
\frac{1}{2}\left[(x_1-x_0)f(x_0)+(x_N-x_{N-1})f(x_N)
+\sum_{i=1}^{N-1} (x_{i+1}-x_{i-1})f(x_i)\right]
\)

と求める事が出来ます。ここで、\(x_0=a,~x_N=b\)です。

s=0d0
do i=1,N
   c=1d0
   if(i.eq.1.or.i.eq.N)c=0.5d0
   x=(i-1)*(b-a)/dble(N-1)+a
   s=s+c*f(x)
enddo
s=s*(b-a)/(N-1)

3次元での台形則は以下のようなプログラムで実装できます。

s=0d0
do ix=1,Nx
   cx=1d0
   if(ix.eq.1.or.ix.eq.Nx)cx=0.5d0
   x=(ix-1)*(xb-xa)/dble(Nx-1)+xa
   do iy=1,Ny
      cy=1d0
      if(iy.eq.1.or.iy.eq.Ny)cy=0.5d0
      y=(iy-1)*(yb-ya)/dble(Ny-1)+ya
      do iz=1,Nz
         cz=1d0
         if(iz.eq.1.or.iz.eq.Nz)cz=0.5d0
         z=(iz-1)*(zb-za)/dble(Nz-1)+za
         s=s+cx*cy*cz*f(x,y,z)
      enddo
   enddo
enddo
s=s*((xb-xa)/dble(Nx-1))*((yb-ya)/dble(Ny-1))*((zb-za)/dble(Nz-1))

3,4桁以上の計算精度が欲しい時、台形則では間に合いません。
刻み幅を小さくすると精度は確実に上がりますが、刻み幅を1/10すると精度が1桁上がるだけなので、なかなか収束しません。
しかし、台形則も馬鹿にはできません。例えば、周期関数の1周期に渡る積分や、両端で定数に近づいていく積分の場合に、台形則はとんでもない精度をたたき出します(たしかガウス求積法と同じ精度まで上がると思いましたが、詳しい文献が見つけられないので参照はありません)。
この考えを元に考え出されたのが二重指数関数型数値積分です。端点特異性がある場合に特に強いため、学んでおいて損は無いでしょう。


スポンサーリンク


積分する関数に性質が見当たらないが、高い精度が欲しい場合


関数の形によって分点の取り方を自動的に決めてくれる適応型の自動積分が良いです。この種の積分で有名なのは、

NetlibにあるQUADPACK
http://www.netlib.org/quadpack/
最速・高精度の数値積分, QUADPACKプログラム(1/2次元, 実/複素積分) -シキノート)、
もしくは
名古屋大学大型計算機センター、Numpacにあるaqnn9d
http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/index.html
を参照してください。

QUADPACKはガウス-ルジャンドル求積法とガウス-クロンロッド求積法による適応型自動積分であり、
aqnn9dはニュートン・コーツ9点則に基づく適応型自動積分です。
おおよそ単精度ではaqnn9dの方が関数の評価回数が少なくて済み、倍精度以上ではQUADPACKの方が少なくて済みます。
この単精度/倍精度の違いは積分アルゴリズムそのものの特性です。

自動積分について、頒布されているfortranコードの良し悪しを最速・高精度の数値積分で比較したので、詳しいことを知りたい方はそちらをご覧ください。

また、刻み幅制御されたルンゲクッタ法による積分も良いでしょう。
定積分
\(
\displaystyle g(x)= \int_a^x f(x’)dx’
\)

を行うには、初期条件\(g(a)=0\)の下で、
\(
\displaystyle \frac{d g(x’)}{dx’} = f(x’)
\)

を点\(x\)まで上の微分方程式を解けば良いのです。

滑らかな関数(※1)であることが分かっている場合


積分区間内に特異性(※2)は無いとします。
この場合はチェビシェフ多項式による展開に基づくクレンショウ-カーティス(Clenshaw–Curtis)積分が優秀です。
定積分を
\(
\begin{align}
\displaystyle \int_{-1}^1 f(x) dx &\sim \sum_{k=0}^{n} \omega_k f(x_k)\\
x_k&=\cos\left(\frac{k\pi}{N}\right),\;\;\;k=0,1,\cdots,N\\
\omega_0 &=\omega_N =\frac{1}{N^2 -1} \\
\omega_s &=\omega_{N-s} = \frac{4}{n^2}\sum^{n/2}_{j=0}{\prime\prime }\frac{1}{1-4j^2}\cos\left(\frac{2j\pi s}{n}\right),\;\;\;s=0,1,\cdots,\frac{N}{2}
\end{align}
\)
ここで、\(\sum”\)は和の最初と最後の項に\(\frac{1}{2}\)を掛けることを意味します。

クレーンショー・カーチス数値積分則 Ooura’s Mathematical Software Packages
または
Clenshaw–Curtis求積法 シキノート
最速・高精度の数値積分
をご覧下さい。

端点に特異性がある場合


計算区間内の途中に特異点は無く、計算区間の端点に特異性がある場合、変数変換型の数値積分である二重指数関数型数値積分が有効です。
具体的には例えば
\(
\displaystyle \int_0^1 \sqrt{x} dx
\)


\(
\displaystyle \int_{-1}^{1} x^{-0.8} dx
\)

の事です。これらの積分を行う際に二重指数関数型数値積分は非常に少ない分点数で非常に高精度の結果を返します。

重みと分点位置は以下の通り決められます。積分を
\(
\begin{align}
\int_{-1}^{1} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\tanh\left(\frac{\pi}{2}\sinh(ih)\right), \\
w_i&=\frac{h\frac{\pi}{2}\cosh(ih)}{\cosh^2\left(\frac{\pi}{2}\sinh(ih)\right)}
\end{align}
\)
として近似します。\(N−,N+\)は離散化誤差と打ち切り誤差が等しくなるように決められます。
\(N−,N+\)はプログラム上では、\(x_i\)がコンピュータの扱える桁数を超えず、\(w_i\)がアンダーフローを起こさない範囲で決められます。

詳しくは二重指数関数型数値積分 -シキノート
をご覧ください。
プログラム自体は大浦様による二重指数関数型数値積分公式の方が優秀なので、こちらを参照すると良いでしょう。

被積分関数が(特定の関数関数)×(多項式)の場合


被積分関数\(f(x)\)が(重み関数\(\omega(x)\))×(多項式\(g(x)\))で記述される特別な場合、ガウス求積法による数値積分が有効です。

多項式\(g(x)\)とは\(x^n\)(nは整数)の足し算で書ける関数です。例えば\(x^4+3x^2+7x\)とかのことです。\(x^2+2\sqrt{x}\)とか、\(x\)の0以上の整数乗ではない項があってはいけません。

重み関数\(\omega(x)\)とは、特定の形をした関数の事です。例えば

\(\omega(x)=1~\to\) ガウス-ルジャンドル求積法、ガウス-ロバート求積法、ガウス-radau求積法
\(\omega(x)=x^n e^{-x}~\to\) ガウス-ラゲール求積法
\(\omega(x)=e^{-x^2}~\to\) ガウス-エルミート求積法
\(\omega(x)=\frac{1}{\sqrt{1-x^2}}~\to\) ガウス-チェビシェフ(第一種)求積法
\(\omega(x)=\sqrt{1-x^2}~\to\) ガウス-チェビシェフ(第二種)求積法
\(\omega(x)=(1-x)^\alpha (1+x)^\beta~\to\) ガウス-ヤコビ求積法

の形を数値積分したい場合、ガウス求積法を考えてみてはどうでしょうか。
上では省略しましたが、ガウス求積法を用いる際には積分区間も決まっています。しかし、これは簡単な変数変換によって書き換えることが出来ます。

例えば、ガウス-ルジャンドル求積法は区間\([-1,1]\)出なければ使うことはできませんが、変数変換
\(
\displaystyle \int_a^b f(x)dx=\frac{b-a}{2}\int_{-1}^{1}f(\frac{b-a}{2}+\frac{a+b}{2}) dx
\)

の変換が出来るので、左辺を計算しようと思ったら右辺を計算しにいけばいいのです。

ガウス-ルジャンドル求積法 http://mathworld.wolfram.com/Legendre-GaussQuadrature.html
ガウス-radau求積法 http://mathworld.wolfram.com/RadauQuadrature.html
ガウス-ロバート求積法 http://mathworld.wolfram.com/LobattoQuadrature.html
ガウス-ラゲール求積法 http://mathworld.wolfram.com/Laguerre-GaussQuadrature.html
ガウス-エルミート求積法 http://mathworld.wolfram.com/Hermite-GaussQuadrature.html
ガウス-チェビシェフ求積法 http://mathworld.wolfram.com/ChebyshevQuadrature.html
ガウス-ヤコビ求積法 http://mathworld.wolfram.com/Jacobi-GaussQuadrature.html

二重積分を行いたい場合


一番理解しやすい方法は直積を考えて計算する方法です。
「二重積分」を「一重積分した結果を一重積分する」という問題に焼直します。
すなわち、
\(
\displaystyle \int dx \int dy f(x,y)
\)

という問題は、
\(
\begin{align}
\int dx \int dy f(x,y)=\int dy\left[\int f(x;y) dx\right]
\end{align}
\)

と考えるのです。ここで、\(f(x;y)\)のセミコロンの意味は、

『\(f\)は\(x\)と\(y\)の変数だよ。だけど、今変数として注目したいのは\(x\)についてだけ。だから\(y\)は一応書いておくけど、今あらわな変数として扱わないよ。』

ということを示しています。
さて、計算を進めましょう。\([~~]\)内は今、固定したある\(y\)について\(x\)のみの変数なので、求積法の形にして
\(
\begin{align}
\int dy\left[\int f(x;y) dx\right]&\approx \int dy\left[\sum_i w^{(x)}_i f(x_i,y)\right]\\
&= \sum_i w^{(x)}_i \int f(x_i,y)dy \\
&\approx \sum_i w^{(x)}_i \sum_j w^{(y)}_j f(x_i,y_j)
\end{align}
\)

となります。ここで、\(w^{(x)}_i\)は\(x\)方向の\(i\)番目の分点の重みを表します。
これは、例えば被積分関数が\(x\)方向については端点特異性を持っていて、\(y\)方向については簡単な多項式で描けている、といった場合に別々の求積法を使うことで計算量を減らし、高精度の結果が得られるということを示しています。

表面積分を行いたい場合


表面積分
\(
\displaystyle \int_0^\pi d\theta \int_0^{2\pi} d\varphi f(\theta,\varphi) \sin\theta \approx 4\pi \sum_{i} w_i f(\theta_i, \varphi_i)
\)

を行う場合、Lebedev求積法かGaussian-quadrature(ガウス-ルジャンドル求積法と台形則の組み合わせ)が良いです。
lebedev求積法は、球面調和関数の直交性を利用した求積法で、\(f(\theta,\varphi)\)が最大\(l=131\)までの球面調和関数で展開されるのであれば、厳密な値を返すという積分法です。

fortranのプログラムはhttp://www.ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/Lebedev-Laikov.F
にあります。このプログラムと

program main
  implicit none
  integer::i,N
  double precision,allocatable::x(:),y(:),z(:),w(:)
  double precision::s
  double precision,parameter::pi=dacos(-1d0)
  double precision,external::f
 
  N=146
  allocate(x(1:N),y(1:N),z(1:N),w(1:N))
  call LD0146(x,y,z,w,N)
  s=0d0
  do i=1,N
     s=s+w(i)*f(x(i),y(i),z(i))
  enddo
  write(6,*)4d0*pi*s
 
  stop
end program main


function f(x,y,z)
  implicit none
  double precision::x,y,z,f
 
  f=x*x+y*y+z*z
 
  return
end function f

を一緒にコンパイルしてください。プログラムでは\(\theta,\varphi\)での指定ではなく、\(x,y,z\)のデカルト座標系のものになっています。動かすと半径1の球面の立体角に関する積分、すなわち\(4\pi\cdot 1^2\)の値を返します。

3. Lebedev angular quadrature -Methods of Numerical Integration.

Gaussian-quadratureは、\(\theta\)方向についてはガウス-ルジャンドル求積法、\(\varphi\)方向については台形則により積分する、という方法です。良く使われるのでgaussian-quadratureと名前がついているようです。
\(\varphi\)方向で台形則が使われる理由は、\(\varphi\)方向について周期的であるため、台形則は非常に高精度の結果を与えるためです。
\(\theta\)方向でガウス-ルジャンドル求積法が用いられるのは、球面調和関数の\(m=0\)がルジャンドル多項式で展開されるためです。他の\(m\ne 0\)についてはルジャンドル多項式は完全系なので、その性質から計算される太郎という思惑なのだと思います。

lebedev求積法による計算は最大\(l=131\)の求積点と重みが知られていますが、これ以上の\(l\)に関してのlebedev求積点、重みに関する論文は見つけられませんでした。
もしも\(l=131\)以上が出てくることが分かっている場合、現状ではGaussian-quadratureによる求積が良いでしょう。

多次元(6,7次元以上)の積分を行いたい場合


この場合はモンテカルロ積分が良いです。

一様乱数によるD次元の積分は、
\(
\begin{align}
\displaystyle &\int_{x^{(1)}_{a}}^{x^{(1)}_{b}}\int_{x^{(2)}_{a}}^{x^{(2)}_{b}}\cdots\int_{x^{(D)}_{a}}^{x^{(D)}_{b}}
f(x^{(1)},x^{(2)},\cdots,x^{(D)})dx^{(D)}\cdots dx^{(2)} dx^{(1)}\\
&\approx
\frac{(x^{(1)}_{b}-x^{(1)}_{a})(x^{(2)}_{b}-x^{(2)}_{a})\cdots(x^{(N)}_{b}-x^{(N)}_{a})}{N}
\sum_{i=1}^N f(x^{(1)}_i,x^{(2)}_i,\cdots,x^{(D)}_i)
+O(\frac{1}{\sqrt{N}})
\end{align}
\)

と求める事が出来ます。

モンテカルロ積分の誤差はサンプル数にのみに依存します。
\(O(\frac{1}{\sqrt{N}})\)なので、サンプル数\(N\)を100倍に増やすと積分精度が\(\frac{1}{\sqrt{100}}=\frac{1}{10}\)、つまり1桁上がる事を意味しています。

一般的な数値積分法、例えば台形則を用いて1次元積分を\(N\)点で積分するとします。
すると、2次元の積分では\(N^2\)点、D次元では\(N^D\)点被積分関数を評価する必要があります。
台形則の組み合わせによる1次元当たりの積分の精度は\(N\cdot O(h^2)\sim O(\frac{1}{N})\)で決まるので、
1桁、積分精度を上げたいとすれば被積分関数の評価回数は\(10^{D}\)倍に増やさなければなりません。

しかしモンテカルロ積分では積分の精度は次元に依らずサンプル数にのみで決定されます。これは、
1桁、積分精度を上げたいとすれば被積分関数の評価回数は\(100\)倍に増やさなければならないことを意味しています。

台形則と比べるとおおよそ3重積分を解こうとする際に、モンテカルロ法の方が有利になる、ということが分かります。

スポンサーリンク

特定の\(x\)だけ関数値\(f(x)\)が得られる場合


”特定の\(x\)”が等間隔であったり、特に整然としていないと仮定します。
この場合は、

台形則による積分
もしくは
3次スプライン補間による積分

がどんな分点数、分点の間隔であっても、プログラムの変更をほとんどしないで計算できるので良いでしょう。

もしも分点が等間隔の場合で、データ配列の大きさが\(1,2,…,2N+1, (Nは整数)\)の場合はシンプソン則による計算でも良いでしょう。
また、分点が等間隔で、データの配列の大きさが\(1,2,…, 2^N+1(Nは整数)\)の場合はロンバーグ積分法が使えます。

3次スプライン補間による積分は
スプライン補間(1, 2次元) -シキノートに置いてあるc3spline_integralを用いると良いでしょう。

注釈

(※1)”滑らか”とは、
積分区間内の至る所で被積分関数の\(x\)のn階微分の右からの極限と左からの極限が一致すること、を意味します。
(※2)”特異性”とは、
積分区間のどこか一点でも被積分関数の\(x\)のn階微分の右からの極限と左からの極限が一致しない点がある、もしくは積分区間の端点で被積分関数の\(x\)のn階微分が無限大に発散する、ことを意味しています。

モンテカルロ積分

モンテカルロ積分は乱数を利用した積分方法です。

モンテカルロ積分の精度を1桁上げるためには評価点の数を100倍に増やす必要があります。
ですが、この性質は次元数に依らないため、多重積分(5~8次元以上など)を行うときに有利な性質を持ちます。

また、とにかく実装が簡単なので非常に使い易いです。

1次元のモンテカルロ積分


1次元のモンテカルロ積分は定積分を
\(
\displaystyle \int_a^b f(x)dx \approx \frac{1}{N}\sum_{i=1}^N\frac{f(x_i)}{p(x_i)}+O(\frac{1}{\sqrt{N}})
\)

として求めます。
ここで、\(x_i\)は\(p(x)\)に従って\(i\)番目に発生させた乱数で、
\(p(x)\)は区間\([a,b]\)内の乱数の分布を表す関数で、
\(
\displaystyle \int_a^b p(x) dx =1
\)

を満たします。一様乱数であれば
\(
\displaystyle p(x)=\frac{1}{b-a}
\)

です。モンテカルロ積分は一様乱数のみで出来るのではなく、任意の乱数分布に対して適用することが出来ます。

一様乱数によるでD次元の積分は、
\(
\begin{align}
&\int_{x^{(1)}_{a}}^{x^{(1)}_{b}}\int_{x^{(2)}_{a}}^{x^{(2)}_{b}}\cdots\int_{x^{(D)}_{a}}^{x^{(D)}_{b}}
f(x^{(1)},x^{(2)},\cdots,x^{(D)})dx^{(D)}\cdots dx^{(2)} dx^{(1)} \\
&\hspace{1em}\approx
\frac{(x^{(1)}_{b}-x^{(1)}_{a})(x^{(2)}_{b}-x^{(2)}_{a})\cdots(x^{(D)}_{b}-x^{(D)}_{a})}{N}
\sum_{i=1}^N f(x^{(1)}_i,x^{(2)}_i,\cdots,x^{(D)}_i)
+O(\frac{1}{\sqrt{N}})
\end{align}
\)

と求める事が出来ます。

モンテカルロ積分の誤差はサンプル数にのみに依存します。
モンテカルロ積分の誤差は\(O(\frac{1}{\sqrt{N}})\)なので、サンプル数\(N\)を100倍に増やすと積分精度が\(\frac{1}{\sqrt{100}}=\frac{1}{10}\)、つまり1桁上がる事を意味しています。

一般的な数値積分法、例えば台形則を用いて1次元積分を\(N\)点で積分するとします。
すると、2次元の積分では\(N^2\)点、D次元では\(N^D\)点被積分関数を評価する必要があります。
台形則を組み合わせて積分する場合、精度は\(N\cdot O(h^2)\sim O(\frac{1}{N})\)で決まるので、
1桁、積分精度を上げたいとすれば被積分関数の評価回数は\(10^{D}\)倍に増やさなければなりません。

しかしモンテカルロ積分では積分の精度は次元に依らずサンプル数にのみで決定されます。これは、
1桁、積分精度を上げたいとすれば被積分関数の評価回数は\(100\)倍に増やさなければならないことを意味しています。

台形則と比べると収束速度は、おおよそ3重積分を解こうとする際に、モンテカルロ法の方が有利になる、ということが分かります。

(注意)
3次元の定積分を考えます。
この時、台形則では100点で真値と1桁、モンテカルロ積分では10000点で真値と2桁合っているとします。
この状態から精度を1桁上げようとすると関数の評価回数は
台形 \(100\times 10^3= 100,000 \)点
モンテカルロ積分 \(10000\times 100\approx 1,000,000\)点
となり、この場合では台形の方が評価回数が少なくて済みます。ずっと繰り返せばモンテカルロ積分の方が少なくなりますが、現実的な評価回数ではなくなると思います。
なので、評価回数が増えていく割合が変わるのは3重積分あたりですが、最終的に評価回数が少なくなるのはもっと多次元の積分を実行する時かもしれません。

台形則と比較しましたが、台形則の代わりにガウス求積法を用いた場合、もう少し高次元で逆転します。

実用的には6~10次元以上でモンテカルロ積分が使われるようです。

乱数について


モンテカルロ積分は乱数を利用して積分を行うため、乱数の性質が重要になってきます。
私は乱数について詳しくありませんが、調べた限りメルセンヌツイスタという乱数発生法が良い乱数を発生するようです。
メルセンヌツイスタの詳細はMersenne Twister Home Page にあるのでご参照ください。
fortran90コードはhttps://gist.github.com/ykonishi/5569005にあります。

fortran90コード

1次元の積分
\(
\displaystyle \int_0^2 \sin x dx \approx 1.416146836547142\cdots
\)

を実際にメルセンヌツイスタを用いて、計算します。
関数の評価点数は10000点です。

下のコードを実行すると

> gfortran -fno-range-check mt19937.f90 main2.f90
> ./a.out
       10000   1.4167388988307803        1.4161468365471424  
>

という出力を得ます。1列目が関数の評価数、2列目がモンテカルロ積分による積分値、3列目が真値です。
1万点の評価で、ただのsin関数を計算しているのに4桁程度というのは経験的に非常に良くない精度です。
しかも、乱数を用いているので、4桁一致していると判断するのは危険です。
もっと乱数を発生させて、値が変わらないのかチェックを欠かすことはできません。

モンテカルロ積分は以下のプログラムで可能です。

program main
  use mt19937
  implicit none
  integer::i,j,N
  double precision::a,b,s,x,y,z,ans
  double precision,external::f

  call sgrnd(1)
 
  a=0d0
  b=2d0
  N=10000
  ans=cos(a)-cos(b)
     
  s=0d0
  do i=1,N
     x=grnd()*(b-a)+a
     s=s+f(x)
  enddo
  s=s*(b-a)/dble(N)

  write(6,*)N,s,ans
 
  stop
end program main

function f(x)
  implicit none
  double precision::f,x
 
  f=sin(x)
 
  return
end function f

コンパイルは上のファイルをmain.f90と名付けた場合

gfortran mt19937.f90 main.f90

とします。もしもエラーが出てしまう場合

gfortran -fno-range-check mt19937.f90 main.f90

とします。
このエラーはgfortranコンパイラ(ver4.8.4で確認)の問題です。
fortranの整数は
-2147483648から2147483647までで定義されているのにもかかわらず、
gfortran(ver4.8.4)は
-2147483647から2147483647の範囲でなければエラーを出すようです。
しかし、実際に整数に-2147483648を代入し、このエラーを無視する上のオプションを付けて実行すると、ちゃんと値が入ります。なので、コンパイラに問題があると判断しました。

3次元積分


3次元の定積分を2つ考えます。
\(
\displaystyle \int_0^1 dz \int_0^1 dy \int_0^1 dx \sin(x)\sin(y)\sin(z) \approx 0.097144222323873819\cdots
\)

\(
\displaystyle \int_0^1 dz \int_0^1 dy \int_0^1 dx 100 e^{x}\sin(30y)z^5 \approx 0.80735242505763782\cdots
\)

比較する積分法はモンテカルロ積分、台形則、quadpackによる自動積分です。
quadpackは15-31点ガウス求積法を用います。3次元関数のfortranによるquadpackコードはこちら

結果は以下の図の方になりました。

quadpackについては余りにも精度が違い過ぎるので図には重ねていません。
それぞれ、関数の評価回数と積分値を載せました。下線は真値と一致している桁です。

quadpack, key=3
3375points, 0.097 144 222 323 873 861
50625points,0.807 352 425 057 637 71

左の図で、台形則では一次元方向の分点数が10倍(3次元なのでグラフ上では1000倍)になると積分の精度が約1桁上昇し、モンテカルロ積分では点数が100倍増えると積分の精度が約1桁上昇しているのが見て取れます。
右の図ではあたかも台形則の方が傾きが急峻そうに見えますが、まだ安定している積分になってなく、傾きが緩やかになっていきそうなことが分かります。

QUADPACKは流石と言いますが、次元が増えると辛いですがまだまだモンテカルロ積分を使わなくてよさそうです。関数の評価点数が5桁なのにほぼ全桁一致しています。

用いたメインのプログラムは以下のものです。これと上で紹介した
メルセンヌツイスタのfortran90プログラム,3次元quadpackプログラムを一緒にコンパイルしてください。

参考文献


モンテカルロ法 (Monte Carlo method)
確率密度関数からモンテカルロ積分まで
メルセンヌツイスタMersenne Twister Home Page
メルセンヌツイスタのfortran90コードhttps://gist.github.com/ykonishi/5569005

スプライン補間(1, 2次元)

3次スプライン補間に関する簡単な説明とfortran90による計算コードです。

概要

・与えられたデータ点すべてを通るように補間
・与えられたデータ点間を3次関数で近似する
・与えられたデータ点間の隣り合う補間関数(スプライン関数)の一階微分、二階微分が連続
・与えられたデータ点の両端は二階微分がゼロと仮定

1次元の場合


区間\([x_i, x_{i+1}]\)の間を3次関数
\(
S_i(x)=a_i(x-x_i)^3+b_i(x-x_i)^2+c_i(x-x_i)+d_i,~(i=0,1,\cdots N-1)
\)

で補間します。補間は与えられたデータ点間の隣り合う区間の補間関数の一階微分と二階微分が一致するように決めます。
すなわち
\(
\displaystyle \frac{d S_i(x_i)}{dx}=\frac{d S_{i-1}(x_i)}{dx} \\
\displaystyle \frac{d^2 S_i(x_i)}{dx^2}=\frac{d^2 S_{i-1}(x_i)}{dx^2} \\
\)

という式が成り立っているとします。与えられたデータ点の端、点\(x_0\)と\(x_N\)では自然境界条件と呼ばれる境界条件
\(
\displaystyle \frac{d^2 S_0(x_0)}{dx^2}=\frac{d^2 S_{N-1}(x_N)}{dx^2}=0
\)

を課します…この境界条件を課す理由は良く分かりません。[1]に端ではない領域で補間が滑らかになるようにこう決める、とあります。良く分かりません。きっと解く都合上こう決めると簡単に解けるからでしょう。実際そうですし…

詳細は[1]が詳しいのでそちらを参照してください。
係数は以下の通り与えられます。
\(
\begin{align}
a_i&=\frac{v_{i+1}-v_i}{6(x_{i+1}-x_i)} \\
b_i&=\frac{v_i}{2}\\
c_i&=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-\frac{1}{6}(x_{i+1}-x_i)(2v_i+v_{i+1}) \\
d_i&=y_i
\end{align}
\)
ここで、\(v_i\)は以下の方程式の解です。
\(
\scriptsize
\begin{eqnarray}
\left(
\begin{array}{ccccc}
2(h_0+h_1)& h_1 & 0 & \cdots & 0 \\
h_1 & 2(h_1+h_2) & h_2 & \cdots & 0\\
0 & h_2 & 2(h_2+h_3) & \cdots & 0\\
\ldots &&&& \\
0& \cdots & 0 & h_{N-2} & 2(h_{N-2}+h_{N-1})
\end{array}
\right)
\left(
\begin{array}{c}
v_1 \\
v_2 \\
v_3 \\
\ldots \\
v_{N-1}
\end{array}
\right)=\left(
\begin{array}{c}
w_1 \\
w_2 \\
w_3 \\
\ldots \\
w_{N-1}
\end{array}
\right)
\end{eqnarray}
\normalsize
\)
ここで、
\(
\begin{align}
h_i&=x_{i+1}-x_i,~ (i=0,1,\cdots,N-1)\\
w_i&=6\left(\frac{y_{i+1}-y_{i}}{h_i}-\frac{y_i-y_{i-1}}{h_{i-1}}\right),~ (i=1,2,\cdots,N-1)
\end{align}
\)
です。

計算コード

Fortranコードを以下に置いておきます。
注意しておきますが、[1]を元にしたコードではありません。[2]を元にしたコードです。

・配列 fdata(0:N),xdata(0:N) は与えられたデータで、既知の値です。
・ c3spline1p は任意の点 x の補間値 f 、一階微分 df 、二階微分 df2 を fdata(0:N),xdata(0:N) を元に計算します。
・ c3spline は任意の点配列 xa(0:M) の補間値 fa(0:M) を fdata(0:N),xdata(0:N) を元に計算します。
・ c3spline_integral はスプライン補間された関数の積分 s を計算します。

c3spline1p を用いて複数の点を計算することもできますが、c3spline に比べて遅いです。
計算時間は  c3spline1p : c3spline = 1 : 0.03です。

2次元の補間


2次元の補間を考えます…が、少し横着します。
「横着」は2次元のスプラインを調べるのが面倒で、思いついてしまったのでこれでいいだろう、という意味です。

格子状に与えられたデータ列\(f(x_i,y_j), i=0,\cdots,N_x, j=0,\cdots,N_y\)があり、点\((x,y)\)での補間された値が欲しいとします。
点\((x,y)\)は、それぞれ
\(
\begin{align}
& x_i\le x \lt x_{i+1} \\
& y_j\le y \lt y_{j+1}
\end{align}
\)

に存在するとします。
この時、考え方は以下の通りです。

まず、\(j\)を固定して\(x\)上の補間値を\(N_y\)個(点\(f(x,y_0),f(x,y_1),f(x,y_2),\cdots,f(x,y_{N_y})\))を得ます。
その後、補間された\(N_y\)個の補間値に対して補間を行い、点\((x,y)\)に補間値\(f(x,y)\)を得ることが出来ます。

ここで問題なのが、\(x\)方向、\(y\)方向のどちらを先に補間するのか?
そしてその2つの結果は変わるのか?ということです。

直感的に考えて良さそうなのは2つの平均を取ってしまう事でしょう。しかし、単純に考えて計算時間が2倍になります。
直ぐには分からなそうなので、プログラムではそれを指定できる形にする、にとどめましょう。

偏微分を求められます。最後に補間した方向の偏微分を得ることが出来ます。

c3spline2d1pは元となるデータ配列を入れると、点(x,y)の補間された値を返します。
c3spline2dは元となるデータ配列を入れると、配列で指定された格子状の点を返します。

上記コードを動かすと、以下の画像を作ることが出来ます。

赤い点は補間に用いたオリジナルのデータ列を表し、黒の線は補間された関数を表します。また、底面のカラーマップは本当の関数との相対誤差を表しています。xが大きいところの端で補間した結果の関数が本当の関数とあまり一致していません。

これは、元にするデータ点の両端において、元の関数の二階微分がゼロという仮定が満たされていないためです。

上の図を再現するgnuplotのコードは以下のものです。

set palette defined ( -1 '#ffffff', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set ticslevel 0
set view 46,40,1,1

set xl "x"
set yl "y"
set cbr[0:1]
unset key
splot "interpolate_data1.d" u 1:2:6 w pm3d at b
replot "interpolate_data1.d" u 1:2:3 w l lw 1
replot "original_data.d" u 1:2:3 w p pt 7 ps 0.5

参考文献

[1]https://www.mk-mode.com/rails/docs/INTERPOLATION_SPLINE.pdf

[2]http://150.19.250.16/MULTIMEDIA/numeanal/node16.html

アイキャッチ画像のフォント:キルゴhttp://getsuren.com/killgoU.html