谁说π难求？盘点圆周率的各种骚操作
引言：π的魅力与困惑
圆周率 π \pi π，一直以来是数学界的明星。从古代的阿基米德到现代计算机科学家，大家都想追求这位“神秘数”的真面目。可是，π 就像个顽皮的孩子，永远不会精确暴露自己。这篇博客将带你领略几种不同的方法来近似 π，从简单到复杂，从经典到另类，让你一步步接近圆周率的真相。
第一种：蒙特卡洛方法 - “用随机点算圆周率”
蒙特卡洛方法是通过随机数来模拟实验的经典方法。在求 π 的过程中，这种方法利用了圆和正方形的面积比值来进行估算。基本思路如下：
- 假设有一个边长为 1 的正方形和一个内接的圆。
- 向正方形中随机撒点，计算这些点中落在圆内的比例。
- 利用比例关系来估算 π 值。
公式如下：
π ≈ 4 × 圆内点数 总点数 \pi \approx 4 \times \frac{\text{圆内点数}}{\text{总点数}} π≈4×总点数圆内点数
MATLAB 代码示例：
matlab">% 参数设置
num_points = 1000000; % 随机点的数量
inside_circle = 0; % 圆内点数计数% 随机点生成
x = rand(1, num_points);
y = rand(1, num_points);% 圆内点数统计
for i = 1:num_pointsif x(i)^2 + y(i)^2 <= 1inside_circle = inside_circle + 1;end
end% 计算 π 的近似值
pi_estimate = 4 * inside_circle / num_points;% 数据可视化
figure;
plot(x, y, '.'); hold on;
theta = linspace(0, 2*pi, 100);
plot(cos(theta), sin(theta), 'r'); % 圆的边界
axis equal;
title(['Monte Carlo Approximation of \pi, Estimated \pi = ', num2str(pi_estimate)]);
xlabel('x');
ylabel('y');
结果分析：
使用蒙特卡洛方法得到的 π 值近似精度取决于随机点的数量。点数越多，估算的 π 值会越接近实际值，但精确度却受限于方法的随机性。
第二种：莱布尼茨公式 - “简单的无穷级数”
莱布尼茨公式是基于无穷级数来逼近 π 的一个简单方法：
π = 4 ∑ k = 0 ∞ ( − 1 ) k 2 k + 1 \pi = 4 \sum_{k=0}^{\infty} \frac{(-1)^k}{2k + 1} π=4k=0∑∞2k+1(−1)k
MATLAB 代码示例：
matlab">% 参数设置
terms = 10000; % 累加项数
pi_estimate = 0; % π 的近似值% 计算 π 的近似值
for k = 0:terms-1pi_estimate = pi_estimate + 4 * ((-1)^k) / (2 * k + 1);
end% 显示结果
disp(['莱布尼茨公式估算的π值为：', num2str(pi_estimate)]);% 数据可视化
n_values = 1:terms;
pi_values = 4 * cumsum((-1).^(n_values - 1) ./ (2 * n_values - 1));
plot(n_values, pi_values);
xlabel('Number of Terms');
ylabel('Estimated \pi');
title('Leibniz Series Approximation of \pi');
grid on;
结果分析：
matlab">莱布尼茨公式估算的π值为：3.1415
虽然莱布尼茨公式的计算方式简单，但其收敛速度极慢。需要大量项数才能获得较高的精度。这种方法更适合作为 π 计算的“入门级”示例。
第三种：高斯-勒让德算法 - “精准的迭代法”
高斯-勒让德算法是一种迭代方法，计算精度极高，非常适合用来计算 π 的高精度值。其迭代公式如下：
a 0 = 1 , b 0 = 1 2 , t 0 = 1 4 , p 0 = 1 a_0 = 1, \quad b_0 = \frac{1}{\sqrt{2}}, \quad t_0 = \frac{1}{4}, \quad p_0 = 1 a0=1,b0=21,t0=41,p0=1
迭代计算：
a n + 1 = a n + b n 2 a_{n+1} = \frac{a_n + b_n}{2} an+1=2an+bn
b n + 1 = a n × b n b_{n+1} = \sqrt{a_n \times b_n} bn+1=an×bn
t n + 1 = t n − p n ⋅ ( a n − a n + 1 ) 2 t_{n+1} = t_n - p_n \cdot (a_n - a_{n+1})^2 tn+1=tn−pn⋅(an−an+1)2
p n + 1 = 2 ⋅ p n p_{n+1} = 2 \cdot p_n pn+1=2⋅pn
最终 π 的近似值为：
π ≈ ( a n + b n ) 2 4 ⋅ t n \pi \approx \frac{(a_n + b_n)^2}{4 \cdot t_n} π≈4⋅tn(an+bn)2
MATLAB 代码示例：
matlab">% 初始值
a = 1;
b = 1 / sqrt(2);
t = 1 / 4;
p = 1;
num_iterations = 10; % 迭代次数% 迭代计算
for n = 1:num_iterationsa_next = (a + b) / 2;b = sqrt(a * b);t = t - p * (a - a_next)^2;p = 2 * p;a = a_next;
end% 计算 π 的近似值
pi_estimate = (a + b)^2 / (4 * t);% 显示结果
disp(['高斯-勒让德算法估算的π值为：', num2str(pi_estimate)]);
结果分析：
matlab">高斯-勒让德算法估算的π值为：3.1416
高斯-勒让德算法以其快速的收敛性闻名，几次迭代即可得到极高精度的 π 值。适合需要精确 π 值的计算场景，但计算过程相对复杂。
第四种：蒲丰投针法古老的“针”与圆周率π的邂逅
第五种：连分数展开 - “π的分数之美”
连分数是一种用分数嵌套来表示数的方法。圆周率 π 的连分数展开近似公式如下：
π = 3 + 1 2 6 + 3 2 6 + 5 2 6 + … \pi = 3 + \frac{1^2}{6 + \frac{3^2}{6 + \frac{5^2}{6 + \dots}}} π=3+6+6+6+…523212
MATLAB 代码示例：
matlab">% 参数设置
terms = 10; % 连分数的项数
pi_estimate = 3;% 计算 π 的近似值
for k = terms:-1:1pi_estimate = 3 + (2 * k - 1)^2 / (6 + pi_estimate);
end% 显示结果
disp(['连分数展开法估算的π值为：', num2str(pi_estimate)]);% 可视化
n_values = 1:terms;
pi_values = zeros(1, terms);
temp_estimate = 3;
for k = 1:termstemp_estimate = 3 + (2 * k - 1)^2 / (6 + temp_estimate);pi_values(k) = temp_estimate;
endplot(1:terms, pi_values, '-o');
xlabel('Number of Terms');
ylabel('Estimated \pi');
title('Continued Fraction Approximation of \pi');
grid on;
结果分析：
matlab">连分数展开法估算的π值为：3.1018
连分数方法提供了一种独特的角度来逼近 π。虽然它的收敛速度较慢，但连分数展开展示了 π 的数学结构之美。
第六种：拉马努金公式 - “近似π的巅峰之作”
印度数学家拉马努金提出了一个以无穷级数形式逼近 π 的公式，其精度和收敛速度非常惊人，广泛应用于高精度计算中。拉马努金的 π 公式如下：
1 π = 2 2 ∑ k = 0 ∞ ( 4 k ) ! ( 1103 + 26390 k ) ( k ! ) 4 39 6 4 k \frac{1}{\pi} = 2\sqrt{2} \sum_{k=0}^{\infty} \frac{(4k)!(1103 + 26390k)}{(k!)^4 396^{4k}} π1=22k=0∑∞(k!)43964k(4k)!(1103+26390k)
这与马青公式有些相似，但通过特殊的处理得到了更快的收敛速度。
MATLAB 代码示例：
matlab">% 参数设置
terms = 5; % 累加项数
pi_inverse = 0;% 计算 π 的倒数
for k = 0:terms-1num = factorial(4 * k) * (1103 + 26390 * k);den = (factorial(k)^4) * 396^(4 * k);pi_inverse = pi_inverse + num / den;
end% 计算 π 的近似值
pi_estimate = 1 / (pi_inverse * (2 * sqrt(2) / 9801));% 显示结果
disp(['拉马努金公式估算的π值为：', num2str(pi_estimate)]);% 可视化
n_values = 1:terms;
pi_values = arrayfun(@(x) 1 / (sum(arrayfun(@(k) factorial(4 * k) * (1103 + 26390 * k) / ((factorial(k)^4) * 396^(4 * k)), 0:x-1)) * (2 * sqrt(2) / 9801)), n_values);
plot(n_values, pi_values, '-o');
xlabel('Number of Terms');
ylabel('Estimated \pi');
title('Ramanujan Series Approximation of \pi');
grid on;
结果分析：
matlab">拉马努金公式估算的π值为：3.1416
拉马努金公式是π计算领域的一个巅峰作品，其收敛速度快得惊人，甚至被用于现代计算机 π 的精确计算中。
总结与反思：π的求值方法之旅
以上我们探索了六种不同方法来逼近圆周率 π，从随机模拟的蒙特卡洛法到收敛神速的拉马努金公式。每种方法展现了数学的不同领域，体现了从几何到级数、从概率到连分数的多样性。
这些方法可以总结为以下几点：
-
简单性与复杂性：例如，蒙特卡洛方法和莱布尼茨公式以其简单的计算方式适合初学者，但精度较低；而高斯-勒让德算法和拉马努金公式虽复杂，但收敛迅速，适合高精度计算。
-
收敛速度：不同方法收敛速度差异巨大。简单的级数求法收敛较慢，而拉马努金公式和马青公式的收敛速度极快。
-
计算成本：一些方法需要大量的计算步骤和迭代（如连分数和莱布尼茨），而另一些方法在相对少的计算量下就能达到高精度（如高斯-勒让德）。
未来展望：
π 的探索仍然没有尽头，未来的数学家可能会发现更多奇妙的计算 π 的方法，也可能会进一步揭示圆周率和宇宙其他恒量之间的关系。
订阅与关注
如果您对 π \pi π感兴趣，欢迎订阅我的博客，获取最新的技术文章和研究动态。也可以通过以下方式与我交流：
- 📧 邮箱：beining_chen@stu.zzu.edu.cn
关于作者
**[Cherngul]**是一位热爱算法与数据科学的技术达人，专注于探索各种优化算法在实际中的应用，思考如何用算法让生活更美好。希望通过分享，让更多人爱上算法的奇妙世界！
标签
- π \pi π
- MATLAB
- 蒙特卡洛方法
- 莱布尼茨公式
- 拉马努金公式
- 马青公式
最后
非常感谢您抽出时间来阅读我的文章！您的意见非常宝贵。文中可能有些地方表达得不够准确或错误，如果您遇到有需要改进或调整的地方。如果有任何问题或建议，请随时告诉我，我会非常感激。再次感谢您的阅读！
版权信息
© 2024 [Cherngul]. 保留所有权利。
