程序实现蒙特卡洛算法计算PI值和积分

原创已于 2022-05-03 22:14:21 修改 · 2.4k 阅读

11 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #概率论

于 2022-01-02 16:19:54 首次发布

嵌入式系统同时被 3 个专栏收录

355 篇文章

订阅专栏

人工智能

155 篇文章

订阅专栏

算法

102 篇文章

订阅专栏

本文介绍了蒙特卡罗方法，一种使用随机数进行数值计算的技术。通过随机采样，该方法可以估算PI值和进行积分计算。示例中展示了如何用C语言和Python实现估算PI的程序，以及利用蒙特卡罗方法进行积分的近似计算。随着样本数量增加，结果趋向准确。此外，文章还探讨了PI的数学特性，包括其无理性和超越性。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数或伪随机数来解决很多计算问题的方法。与它对应的是确定性算法，蒙特卡洛是摩纳哥的一座城市，注意是摩纳哥不是有卡萨布兰卡的摩洛哥，蒙特卡洛以赌博闻名，这大概也是其名称被用来作为一种概率统计算法的名称的缘由把，毕竟都与概率现象有关。

由于很多情况下计算实际值会消耗大量时间，因此蒙特卡洛方法的思想是通过随即采样来对实际值进行估算，并且可以得到良好的估算结果。它的一个简单示例是用蒙特卡洛方法估算PI值，如下图：

下图中，是一个内切于正方形的单位圆，如果按照概率计算，随机生成的坐标点落在圆形的概率与落在整个正房形区域内的面积比为：

$Prob=\frac{\pi r^2}{(2r)^2}=\frac{\pi}{4} \Rightarrow \pi= 4\cdot Prob$

我们编程实现：

#include <stdio.h>
#include <stdlib.h>

int main(void)
{
    int tries = 0;
    int success = 0;

    for(tries = 0; tries < 80000000; tries ++)	
    {
        float x=rand()/(RAND_MAX+1.0);
        float y=rand()/(RAND_MAX+1.0);
        
        if((x*x + y*y) <= 1)
        {
            success ++;
        }
    }

    double pi = 4 * (double)success / (double)tries;
    printf("pi = %f\n", pi);

    return 0;
}

蒙特卡洛求积分

为了理解蒙特卡洛积分的含义，这里有一个可视化图，对于从区间[a, b]中均匀选取的每一个x值，我们计算f(x)。用(b-a)乘以f(x)得到矩形的面积，如上图所示。通过计算这些面积的平均值，我们可以得到f在区间[a, b]下曲线的近似面积。其理论基础是大数定律。

手搓积分验证大数定律和中心极限定理_papaofdoudou的博客-优快云博客

辛钦大数定律的Python实践_papaofdoudou的博客-优快云博客

说白了，就是随着样本增大，样本均值无限趋近于真是的均值。

import random

a = 0;
b = 1;

N=200000000;

def f(x):
    return 2*x**5;

def uniform():
    return random.uniform(a, b)

def compute():
    M=f(uniform())
    # for k in range(N):
        # M += ((f(uniform()) - M) / (k + 1))
    # return M*(b-a)
    sum = 0
    res = 0
    for k in range(N):
        sum += f(uniform())
    res = sum / k;
    return res * (b-a);


integral = compute()

print(integral)

运行程序得到近似结果：

$I = \int_{0}^{1}2x^5dx=\frac{1}{3}= 0.333333\cdots$

理论值和蒙特卡洛算法得到的值相差无几。

PI的精确值

PI的精确值是多少呢，很遗憾，关于该取值的消息相当糟糕，由于PI是无理数，因此，我们不可能将其表示为两个整数的比值，特别是，想要将圆的直径和圆周都表示为同一个计量单位的整数倍，是绝对不可能的。

虽然 $\sqrt{2}$ 也是无理数，但是我们至少可以这样描述它，即“其平方为2的数”，换句话说，我们可以用整数的算术来表达 $\sqrt{2}$ 所满足的关系，也就是，它是这样一个数，满足 $x^2=2$ ,虽然我们也不知道它的取值到底是多少，但我们知道它的性质。结果表明，PI极为不同，它不仅不能够用分数表示，事实上，他也不能满足任何的代数关系，PI有什么用呢？除了用来表示圆周率之外，其实它并没有什么别的作用。PI就是PI，像PI这样的数，被称为超越数，超越数有很多，除了PI，还有e,超越数根本就超出了代数所具有的表达能力。真的很神奇，人类居然还能够知道像超越数这样的数！

PI还有其他的表示方法，比如下面的公式,它是莱布尼兹发现的：

$\frac{\pi}{4}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\frac{1}{11}+\cdots$

用程序表达：

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv)
{
	double sum = 0;
	int i = 0;
	int b = 1;

	for(i = 0; i < 200000000; i ++)
	{
		sum += (double)1/(double)(2*i + 1) * b;
		b = -b;
	}

	printf("%s line %d, pi is %f.\n", __func__, __LINE__, 4*sum);

	return 0;
}