double def_integral1(double(*f)(double), double a, double b, double trial) {
double x, y;
double ya, yb;
double ymax;
double in = 0;
ya = f(a); //e^(-a)=1 , a=0.0
yb = f(b); //e^(-b)=1/e , b=1.0
if (ya < yb) ymax = yb;
else ymax = ya;
// a <= x <= b
// 0 <= y <= ymax
for (double i = 1; i <= trial; i += 1) {
x = a + urnd() * (b - a); // a+乱数*(b-a)
y = urnd() * ymax;
if (x*x+y*y<=1.0) in += 1;[/color] // -e^(-b)-[-e^(-a)]
}
return (in/trial)*4;
}
int main(void) {
double a, b;
srand((unsigned)time(NULL));
long c1, c2;
double d;
a = 0.0;
b = 1.0;
printf("tval = %12.10lf\n", f1i(b) - f1i(a));
c1 = clock();
d = def_integral1(f1, a, b, MAX1);
c2 = clock();
printf("sval_1 = %12.10lf (%ld)\n", d, c2 - c1);
return 0;
モンテカルロ法を用いて定積分をしたい
モンテカルロ法を用いて定積分をしたい
学校の課題で、モンテカルロ法を用いて定積分の値を求めたいのですが、以下のsvalのところが分かりませんでした。理論値が tval =1-1/e になると考えていますが、モンテカルロ法でどうやってネイピア数を求めるのか教えていただけたら幸いです。円周率をモンテカルロ法で求めるやり方は知っています。コメントは自分のメモ程度に使っているので気にしないでいただけたら....
Re: モンテカルロ法を用いて定積分をしたい
タカチャン さんが書きました: ↑4年前学校の課題で、モンテカルロ法を用いて定積分の値を求めたいのですが、以下のsvalのところが分かりませんでした。理論値が tval =1-1/e になると考えていますが、モンテカルロ法でどうやってネイピア数を求めるのか教えていただけたら幸いです。円周率をモンテカルロ法で求めるやり方は知っています。コメントは自分のメモ程度に使っているので気にしないでいただけたら....
コピペに不足部分があったので出し直しました。double def_integral1(double(*f)(double), double a, double b, double trial) { double x, y; double ya, yb; double ymax; double in = 0; ya = f(a); //e^(-a)=1 , a=0.0 yb = f(b); //e^(-b)=1/e , b=1.0 if (ya < yb) ymax = yb; else ymax = ya; // a <= x <= b // 0 <= y <= ymax for (double i = 1; i <= trial; i += 1) { x = a + urnd() * (b - a); // a+乱数*(b-a) y = urnd() * ymax; if (x*x+y*y<=1.0) in += 1;[/color] // -e^(-b)-[-e^(-a)] } return (in/trial)*4; } double f1(double x) { return exp(-x); } double f1i(double x) { // f1の積分 return -exp(-x); } int main(void) { double a, b; srand((unsigned)time(NULL)); long c1, c2; double d; a = 0.0; b = 1.0; printf("tval = %12.10lf\n", f1i(b) - f1i(a)); c1 = clock(); d = def_integral1(f1, a, b, MAX1); c2 = clock(); printf("sval_1 = %12.10lf (%ld)\n", d, c2 - c1); return 0;