利用蒙特卡罗方法求积分

函数y=f(x)=sin(x)/x 在区间[1,2]的积分，refer to http://www.cnblogs.com/leoo2sk/archive/2009/05/29/1491526.html

 

代码如下：

 

#include<iostream> #include<cmath> #include<ctime> #define COUNT 500000 using namespace std; bool InCircle(double x, double y) { if( sin(x) / x >= y ) { return true; } else { return false; } } int main() { double x, y; int a = 1, b = 2; int num = 0; int i; srand((unsigned)time(NULL)); for(i = 0; i < COUNT; i++) { x=rand() * 1.0 / RAND_MAX; y=rand() * 1.0 / RAND_MAX; x = a + (b - a) * x; y = (b - a) * y; if(InCircle(x,y)) { num++; } } double res = (double)num / COUNT; cout << "result is " << res << endl; }

 

结果如下：

 

/home/a/j/nomad2:./a.out result is 0.658302

 

思想：将面积转化为求落在函数下方的点与长度为(b-a)的正方形内的点数之比，类似的还有求圆周率。

 

#include<iostream> #include<cmath> #include<ctime> #define COUNT 500000 using namespace std; bool InCircle(double x,double y) { if((x*x+y*y)<=1)return true; return false; } int main() { double x,y; int num=0; int i; srand((unsigned)time(NULL)); for(i=0;i<COUNT;i++) { x=rand()*1.0/RAND_MAX; y=rand()*1.0/RAND_MAX; if(InCircle(x,y)) num++; } cout<<"PI:"<<(num*4.0)/COUNT<<endl; }

 

运行结果：

 

/home/a/j/nomad2:./a.out PI:3.14107