蒙特卡洛方法计算圆周率的一个高精度公式

发表于 更新于

通过文章蒙特卡洛方法计算圆周率：积分公式选择的方差分析， 我们可以得到提高蒙特卡洛方法计算精度的方法：

  1. 确保函数的平滑性，这样可以减少方差，于是每次随机点都更加接近平均值。
  2. 减少积分区间，这样同样的随机点数，采样密度更高，有效点更多，所以误差更小。

半角公式的选择

本文就来验证这两条规则，于是仍然选择反正弦三角函数，但是使用半角公式不断减小积分区间, 于是我们选择了\(\frac{\pi}{24576}\)作为积分值，即取积分:

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

其中\(b\)为

\[ b=\frac{1}{2}\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{3}}}}}}}}}}}} \]

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
26
#! /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=1000):
    # 在区间[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 24576*integral_estimate

# 使用Monte Carlo方法估计积分
a, b = 0, np.sqrt(2-np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(2+np.sqrt(3)))))))))))))/2  # 定义积分区间
estimated_integral = monte_carlo_integral(f, a, b)

print("Pi", estimated_integral)

上述程序计算了\(\frac{\pi}{6}\)的\(\frac{1}{2^{12}}\), 随机点只取了1000, 运行结果为 3.1415926538918395, 这个精度已经相当高了，计算量也急剧下降。所以对比之前的文章我们通过半角公式降低的计算量的同时，有效提高了计算精度。这个例子说明，对于同一个物理或数学问题，设计不同的算法，给出不同的公式是相当重要的事情。