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:モンテカルロ法の積分教えてください。
2重積分は、2重ループで計算できますよ。
普通の積分はとけますか?
小さくきって、四角形を作っていき、その合計を出す方法。
2重積分は、1方向をとめて、1方向ずつ、加算していきます。
ちょっと出かける用事がありますので、帰宅後、お答えします。
すみません。
普通の積分はとけますか?
小さくきって、四角形を作っていき、その合計を出す方法。
2重積分は、1方向をとめて、1方向ずつ、加算していきます。
ちょっと出かける用事がありますので、帰宅後、お答えします。
すみません。
Re:モンテカルロ法の積分教えてください。
ある枠内に一様に乱数の点を発生させて、その図形の内部に入る割合によって面積を求める方法なんです。
普通の積分は何とかできます。
忙しい中すいません。
普通の積分は何とかできます。
忙しい中すいません。
Re:モンテカルロ法の積分教えてください。
う、、モンテカルロ法昔やったんですけど、
すっかり忘れてしまいました・・。
いまや部分積分すら危うい予感・・;
調べてみて
http://momonga.t.u-tokyo.ac.jp/~ooura/intde-j.html
この辺が近いかな?と思うんですけど、どうでしょうか。
力不足でごめんなさい(_ _||)
すっかり忘れてしまいました・・。
いまや部分積分すら危うい予感・・;
調べてみて
http://momonga.t.u-tokyo.ac.jp/~ooura/intde-j.html
この辺が近いかな?と思うんですけど、どうでしょうか。
力不足でごめんなさい(_ _||)
Re:モンテカルロ法の積分教えてください。
初めまして、kuzといいます。
最近Cから離れてしまっていて、あっている自信は
ありません(というか間違えてる可能性が高いです)
しかもスマートじゃありません・・・・
ただ、お答えがない様ですので投稿いたします。
参考程度に見ていただければ幸いです
まず、Rの両辺を36で割ると
X^2/3^2 + Y^2/2^2 <=1 となるので
これは楕円の方程式をあらわします。
ここから-3<=X<=3,-2<=Y<=2
であることがわかります。
この範囲内の乱数の組み合わせで
Rを満たすものの組をf(x)に代入した結果の
総計を組数で割った値が近似値となります。
乱数を出力する所が見苦しいですが、以下ソースです
最近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:モンテカルロ法の積分教えてください。
> 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π になるはずです。
> 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:モンテカルロ法の積分教えてください。
すいません説明不足でした。
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で少数の乱数を発生させる方法が他に考えつかなかったので、
こんな、見苦しいことをしています、申し訳ありません
負についは確かに不要ですね。範囲内の乱数という意味で
負もあるほうが良いかと思ったのですが、蛇足でした。
まだ、わかりにくい点多々ありそうですが、
あくまで参考程度ということで、ご容赦ください
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で少数の乱数を発生させる方法が他に考えつかなかったので、
こんな、見苦しいことをしています、申し訳ありません
負についは確かに不要ですね。範囲内の乱数という意味で
負もあるほうが良いかと思ったのですが、蛇足でした。
まだ、わかりにくい点多々ありそうですが、
あくまで参考程度ということで、ご容赦ください