蒙特卡洛方法求Pi的MATLAB实现和C语言实现

matlab实现

clear
k = 0; %k用于随机点落在1/4圆内计算
for j = 1:1000000 %样本格数
    a = rand(1,2);
    if a(1)^2+a(2)^2<=1 %检验随机数是否满足
    k = k+1;
    end
end
PI = 4*k/j

c实现

/* 利用蒙特卡洛算法近似求圆周率Pi*/ 
/*程序使用：VC++6.0 */ 
#include<stdio.h> 
#include<math.h> 
#include<stdlib.h>  
#define COUNT 800 /*循环取样次数，每次取样范围依次变大*/  
void main() 
{ 
double x,y; 
int num=0; 
int i;  
for(i=0;i<COUNT;i++) 
{
	x=rand()*1.0/RAND_MAX;/*RAND_MAX=32767，包含在<stdio.h>中*/
	y=rand()*1.0/RAND_MAX;
	if((x*x+y*y)<=1)
		num++; /*统计落在四分之一圆之内的点数*/
} 
printf("Pi值等于：%f\n",num*4.0/COUNT); 
}

说明

求pi

在数值积分法中，利用求单位圆的1/4的面积来求得Pi/4从而得到Pi。单位圆的1/4面积是一个扇形，它是边长为1单位正方形的一部分。只要能求出扇形面积S1在正方形面积S中占的比例K=S1/S就立即能得到S1，从而得到Pi的值。怎样求出扇形面积在正方形面积中占的比例K呢？一个办法是在正方形中随机投入很多点，使所投的点落在正方形中每一个位置的机会相等看其中有多少个点落在扇形内。将落在扇形内的点数m与所投点的总数n的比m/n作为k的近似值。P落在扇形内的充要条件是
$$x^2+y^2<=1$$

运行代码

可以看出:随着点数的增加，求得的Pi值渐渐接近真实值。

如果加入程序：srand(time(NULL)); ，同时循环取样次数一定，让取样结果随时间变化，当取样次数为80000000时，可得6次的结果显示：

3.141290 3.141400 3.141268 3.141484 3.141358 3.141462