蒙特卡洛方法计算圆周率

发表于 更新于

使用蒙特卡洛方法求解圆周率的一个经典方法是取一个单位圆和外接正四边形，则正四边形的面积为1, 随机生成点数，统计落在圆内的点的数目m, 然后再以m除以总点数n,获得随机点落在圆内的概率，再乘以正方形的面积就得到了圆周率。 但是本文不计划采用这个方法，我们使用积分来计算，即

\[\pi=2\int_0^1 \frac{1}{\sqrt{1-x^2}} dx\]

按积分中值定理，在\((0,1)\)内必存在一点\(\xi\)满足

\[\pi=2\frac{1}{\sqrt{1-\xi^2}}(1-0)\]

如果我们任意取\((0,1)\)中的一个点，那必然不一定满足上式，但是按照大数定理，当在\((0,1)\)内随机取值的数量达到足够大，并且平均时，按统计规律积分值将接近真实结果，并且随随机取点的数量增大而提高精度。

Python 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#! /usr/bin/env python3
# vim:fenc=utf-8
import numpy as np
# 定义函数f(x)
def f(x):
    return 1/(1-x**2)**0.5

# Monte Carlo方法求积分
def monte_carlo_integral(f, a, b, num_samples=90000000):
    # 在区间[a, b]内随机生成num_samples个样本点
    samples = np.random.uniform(a, b, num_samples)
    
    # 计算这些样本点对应的函数值
    y_values = f(samples)
    
    # 根据平均高度估计积分值，乘以区间的宽度得到近似面积
    integral_estimate = (b - a) * np.mean(y_values)
    
    return 2*integral_estimate

# 使用Monte Carlo方法估计积分
a, b = 0, 1  # 定义积分区间
estimated_integral = monte_carlo_integral(f, a, b)

print("Pi", estimated_integral)

Rust 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
use rand::distributions::{Distribution, Uniform};
use std::f64;

fn f(x: f64) -> f64 {
    1.0 / (1.0 - x.powi(2)).sqrt()
}

fn monte_carlo_integral(a: f64, b: f64, num_samples: usize) -> f64 {
    let between = Uniform::new(a, b);
    let mut rng = rand::thread_rng();
    
    let mut sum = 0.0;
    for _ in 0..num_samples {
        let sample = between.sample(&mut rng);
        sum += f(sample);
    }
    
    let integral_estimate = (b - a) * (sum / num_samples as f64);
    2.0 * integral_estimate
}

fn main() {
    let a = 0.0;
    let b = 1.0;
    let num_samples = 90_000_000;
    
    let estimated_integral = monte_carlo_integral(a, b, num_samples);
    
    println!("Pi: {}", estimated_integral);
}

在这个例子中，我们使用了 rand crate 来生成随机数。因此，在你的 Cargo.toml 文件中需要添加以下依赖：

1
2
[dependencies]
rand = "0.8"

使用cargo build --release编译后再运行代码，可以明显感觉到比Python要快的多。

  • 运行环境：ArchLinux，ThinkPad T490, Intel(R) Core(TM) i5-8265U (8) @ 3.90 GHz，NVIDIA GeForce MX250 [Discrete]，Memory 16G
  • Python 执行总用时1.255s
  • Rust 执行总用时0.466s
  • 它们的精度都能达到3.14xxx, 通过这个例子也显示出了Rust的优越性，在这个例子它需要的时间只是Python的0.37倍。