ページ 1 / 1
モンテカルロ法の積分教えてください。
Posted: 2006年11月30日(木) 23:59
by オス豚
2重積分なんですけど、
∬f(x)dxdy
ただしf(x)=√(9-x^2-9/4*y^2)
R={(x,y); 4x^2+9y^2≦36}
という問題です。Rというのは∬の右下に小さく書いてあるものです。
f(x)はなんとか自分でも作れるんですが、2重積分の計算ができません。
Re:モンテカルロ法の積分教えてください。
Posted: 2006年12月01日(金) 00:21
by 管理人
2重積分は、2重ループで計算できますよ。
普通の積分はとけますか?
小さくきって、四角形を作っていき、その合計を出す方法。
2重積分は、1方向をとめて、1方向ずつ、加算していきます。
ちょっと出かける用事がありますので、帰宅後、お答えします。
すみません。
Re:モンテカルロ法の積分教えてください。
Posted: 2006年12月01日(金) 00:36
by オス豚
ある枠内に一様に乱数の点を発生させて、その図形の内部に入る割合によって面積を求める方法なんです。
普通の積分は何とかできます。
忙しい中すいません。
Re:モンテカルロ法の積分教えてください。
Posted: 2006年12月01日(金) 02:37
by 管理人
う、、モンテカルロ法昔やったんですけど、
すっかり忘れてしまいました・・。
いまや部分積分すら危うい予感・・;
調べてみて
http://momonga.t.u-tokyo.ac.jp/~ooura/intde-j.html
この辺が近いかな?と思うんですけど、どうでしょうか。
力不足でごめんなさい(_ _||)
Re:モンテカルロ法の積分教えてください。
Posted: 2006年12月01日(金) 11:59
by オス豚
何とかできるかがんばって見ます。ありがとうございました。
Re:モンテカルロ法の積分教えてください。
Posted: 2006年12月02日(土) 17:19
by kuz
初めまして、kuzといいます。
最近Cから離れてしまっていて、あっている自信は
ありません(というか間違えてる可能性が高いです)
しかもスマートじゃありません・・・・
ただ、お答えがない様ですので投稿いたします。
参考程度に見ていただければ幸いです
まず、Rの両辺を36で割ると
X^2/3^2 + Y^2/2^2 <=1 となるので
これは楕円の方程式をあらわします。
ここから-3<=X<=3,-2<=Y<=2
であることがわかります。
この範囲内の乱数の組み合わせで
Rを満たすものの組をf(x)に代入した結果の
総計を組数で割った値が近似値となります。
乱数を出力する所が見苦しいですが、以下ソースです
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#define TRY_CNT (10000) //最大試行回数
void main(void){
double randX,randY,ans=0.0f;
int i;
for(i=0;i<TRY_CNT;i++){
srand(time(NULL)); //乱数初期化
for(;;){ //条件に合致するまで
randX = (double)(rand() % 3001);
randY = (double)(rand() % 2001);
if(rand() % 2){
randX = randX / 10000;
}else{
randX = randX / 100000;
}
if(rand() % 2){
randY = randY / 10000;
}else{
randY = randY / 100000;
}
if(rand() % 2){
randX = randX * (-1);
}
if(rand() % 2){
randY = randY * (-1);
}
if(4*randX*randX + 9*randY*randY <= 36){
break;
}
}
ans += sqrt(9-randX*randX - 9/4*randY*randY); //試行回数分集計
}
ans /= TRY_CNT; //試行回数で割る
printf("ans=%f",ans);
}
Re:モンテカルロ法の積分教えてください。
Posted: 2006年12月02日(土) 23:23
by box
> randX = (double)(rand() % 3001);
> randY = (double)(rand() % 2001);
randXは0~3000、randYは0~2000の範囲の値をとります。
次の、10000または100000で割っている理由がよくわかりません。
いずれにせよ、randXは高々0.3、randYは高々0.2です。
その次に符号を変えるロジックがありますが、
後で2乗しているので不要だと思います。
> if(4*randX*randX + 9*randY*randY <= 36){
> break;
randXの2乗は高々0.09、randYの2乗は高々0.04です。
したがって、ここは必ずbreakします。
forループを回る回数は、常に1回です。
> ans += sqrt(9-randX*randX - 9/4*randY*randY); //試行回数分集計
ansの値は、3にひじょうに近くなります。
求める二重積分の値が楕円の面積だとすると、
π*長径*短径=6π になるはずです。
Re:モンテカルロ法の積分教えてください。
Posted: 2006年12月03日(日) 15:42
by kuz
すいません説明不足でした。
R(範囲)が楕円の方程式になっているというこです。
ここでRというのはx,yが満たすべき条件です。
そこで、本来なら適当に乱数を発生させ、
それがRを満たすかどうかを判断し、
満たした場合f(x)に代入し、その総計
(TRY_CNT分)を取得します。そして集計を
TRY_CNTで除算した値が結果となります。
ただ、Rが楕円であることから発生させるべき
乱数の範囲がそれぞれ-3<x<3,-2<y<2ということなのです。
そこでその範囲内の乱数を発生させるために、
randXの場合では2000までの値を1000または10000で
除算することで0.001~2.000での乱数を発生させています。
Cで少数の乱数を発生させる方法が他に考えつかなかったので、
こんな、見苦しいことをしています、申し訳ありません
負についは確かに不要ですね。範囲内の乱数という意味で
負もあるほうが良いかと思ったのですが、蛇足でした。
まだ、わかりにくい点多々ありそうですが、
あくまで参考程度ということで、ご容赦ください
Re:モンテカルロ法の積分教えてください。
Posted: 2006年12月03日(日) 16:26
by オス豚
ありがとうございます。
レベルの高い話しであまりついていけませんけども笑、参考にさせてもらってできました。