Mr.Bo，爱好科技，航天迷一枚，工科研究僧 阅读原文 这个问题其实也等价于在计算机中如何表示 ，对于半径为 1 的单位圆来说，它的面积就是 。 那么估计单位圆的面积就转换为了估计 ，那么先缩小本问题的目标，如何计算单位圆的面积？ （1）最基础的方法就是使用几何概率，即蒙特卡罗方法： 即撒点，假设向正方形中均匀撒入 n 个点，共有 m 个点落入圆中，那么落入圆中点的概率就为： 这样很容易就得到了圆的面积为 ，发现没有，成功绕开了 算出了圆的面积，并且在单位圆情况下能够成功估计 ，对于概率来说，实验次数越多越准，也就意味着撒点数量越多，圆的面积越精准，也意味着 的精度越高。 （2）第二种方法，就是离散，这是非常重要的思想 首先就是微积分的思想，即将圆面积拆分成很多微小面积之和，通过对微小面积求和来逼近真实面积。很典型的例子就是将正方形划分为一个一个的小正方形格子，统计落在圆内的所有正方形小格子的比例。高赞回答举得一系列例子，其实也就是这个思想的体现。在高中其实就已经讲过，曲边梯形面积的求法，这就是微积分最核心的思想，无论理论多么高大上，理论基础都是这个：化曲为直，无限逼近。 图源百度 其次就是级数的思想，即将对单位圆面积的求解划分为一小块一小块面积的总和，很典型的一个级数就是： 包括传统的割圆法等等，体现的核心思想就是两个字：逼近，无非就在于逼近的精度如何罢了。 在完成单位圆面积的估算之后，其他所有圆的面积，无非就是单位圆面积的倍数了，这仅仅与圆的直径 / 半径比值相关，那么就简单通过直尺来得到圆的面积，绕开了 。 在实际使用中，根据所需的精度对 进行估计即可，比如 c 语言的 math 库函数里直接对 值定义，也算典型的直接拿来用（狗头） #if defined(_USE_MATH_DEFINES) && !defined(_MATH_DEFINES_DEFINED) #define _MATH_DEFINES_DEFINED #define M_E 2.71828182845904523536 #define M_LOG2E 1.44269504088896340736 #define M_LOG10E 0.434294481903251827651 #define M_LN2 0.693147180559945309417 #define M_LN10 2.30258509299404568402 #define M_PI 3.14159265358979323846 #define M_PI_2 1.57079632679489661923 #define M_PI_4 0.785398163397448309616 #define M_1_PI 0.318309886183790671538 #define M_2_PI 0.636619772367581343076 #define M_2_SQRTPI 1.12837916709551257390 #define M_SQRT2 1.41421356237309504880 #define M_SQRT1_2 0.707106781186547524401 #endif /* _USE_MATH_DEFINES */ 那么蒙特卡罗方法的精度如何？写个代码便知： sample_points = 10.^(2:1:8); estimate_pi = zeros(length(sample_points),1); for i = 1:length(sample_points) n= sample_points(i); x = rand(n,1); rng('shuffle'); y = rand(n,1); rng('shuffle'); s = sum(x.^2+y.^2 <= 1); estimate_pi(i) = s/n*4; end err=abs(estimate_pi-pi) plot(err,'r','LineWidth',1); xlabel('取值点（10^x）'); ylabel('估计误差'); 下图展示了在 点的情况下，估计值为 3.14123220000000，这显然与真实值差距较大，当然这是由于 MATLAB 的随机函数不是真随机导致的，但是蒙特卡罗方法是随机策略，每次运行的值都不同，因此很少在计算机中使用，常用的方法一般为数值积分方法。 再用级数 来试试，代码为： n=1000; y=0; num=0; for i=1:1:n num=num+(1/(i*i)); end pi=sqrt(6*num); disp(pi) 结果为 3.140638056205995，依旧不太理想。 如果采用对 进行数值积分，代码为： a=0;b=1;s=0;n=1000;i=0; h=(b-a)/n; for i=0:n-1 xi=a+i*h; yi=1/(1+(xi)^2); xj=a+(i+1)*h; yj=1/(1+(xj)^2); s=s+(yi+yj)*h/2; end estimate_pi=vpa(4*s,30); 那么可以得到 为 3.141592486923，显然精度已经很高了. 如果还想要更高的精度，就需要迭代速度更快，收敛速度更快的方法，感兴趣可以去查一下，在这里就不再过多赘述了。 阅读原文
