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

本文介绍了蒙特卡罗方法,一种使用随机数进行数值计算的技术。通过随机采样,该方法可以估算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;
}

这段程序说明,PI可以表示为无穷级数之和,这个公式至少向我们提供了PI的纯数值表示,这样的表示就是我们所能得到的全部。圆周和直径的比值是PI,然而,对于这样的比值,我们却无能为力,我们所能做的,只能是将它加入从而扩展我们的语言。

总结:

蒙特卡洛方法是通过随机才阳来将一个复杂问题化简,这里用蒙特卡洛方法估算PI值,画了一个正方形和它的内切圆,然后在正方形中随机撒下了一些点,看在圆内的点有多少,在正方形内的点又有多少,由二者的比值即可求得PI的估计值。


结束

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

papaofdoudou

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值