2015年8月30日日曜日

二重積分を R と Maxima で解く

二重積分の積分区間が定数でなく,変数の場合,統計解析ソフトRで数値解を求めるには,どうすれば良いか?

最近,Rで重積分というウェブページを見て,このような問題の解説が,ほとんど無いことに気づいた。そのページに書かれていたのは,以下の式である。


二重積分をRとMaximaで解く。






二重積分と言っても驚くことはなく,この場合,通常の単積分を2回繰り返せば良い。

まず, Maxima のプログラム。

f(y):= exp(y^2);
g(x):=integrate(f(y), y, x, 2);
integrate(g(x), x, 0, 2);
float(%);


結果は,以下のとおり。

二重積分を Maxima で解く

26.79907501657212

解析的に解いた結果と数値積分の結果を返してくれる。なお,この数式は, Maxima の出力オプションで,テキスト,画像, LaTeX などある中で, LaTeX を用いて,それを Google chart tools を媒介にして出力したものである。

次に, R による数値解の求め方。

g<- function(x){
  f<- function(y){
  exp(y^2)
  }
  integrate(f, x, 2)$value
}
integrate(Vectorize(g), 0, 2)

注意する点は,二回目の積分,つまり,関数 g の積分で, Vectorize 関数を使い,関数をベクトル化することである。

結果は,
26.79908 with absolute error < 1.3e-12

詳細には検討していないが,実質的に Maxima と同じになる。

参照サイト