探索圆周率π的计算：从古典算法到现代代码实现
探索圆周率π的计算：从古典算法到现代代码实现
圆周率π，这个神秘而又重要的数学常数，自古以来就吸引着无数数学家和爱好者的目光。它不仅在数学领域有着广泛的应用，在物理、工程等众多学科中也扮演着关键角色。今天，我们就来深入探讨计算π的主要算法，并分别用Matlab和Python实现它们。
一、古典算法：割圆术
割圆术是中国古代数学家刘徽的杰出创造。其原理是通过用圆内接正多边形和外切正多边形来逼近圆。随着正多边形边数的不断增加，它的周长就会越来越接近圆的周长，从而计算出π的近似值。
Matlab实现
function pi_approx = liu_hui_pi(n)
% n为正多边形的边数
theta = 2 * pi / n;
s = sin(theta / 2);
pi_approx = n * s;
end
使用示例：
% 计算正1000边形逼近的π值
approx_pi = liu_hui_pi(1000);
disp(['使用割圆术，1000边形逼近的π值为：', num2str(approx_pi)]);
Python实现
import math
def liu_hui_pi(n):
theta = 2 * math.pi / n
s = math.sin(theta / 2)
pi_approx = n * s
return pi_approx
# 计算正1000边形逼近的π值
approx_pi = liu_hui_pi(1000)
print(f"使用割圆术，1000边形逼近的π值为：{approx_pi}")
二、分析算法
（一）无穷级数算法
- 莱布尼茨公式：
π
/
4
=
1
−
1
/
3
+
1
/
5
−
1
/
7
+
⋯
+
(
−
1
)
n
−
1
/
(
2
n
−
1
)
\pi/4 = 1 - 1/3 + 1/5 - 1/7 + \cdots + (-1)^{n - 1}/(2n - 1)
π/4=1−1/3+1/5−1/7+⋯+(−1)n−1/(2n−1)
这个公式形式简洁，但收敛速度较慢，需要计算很多项才能获得较高精度的π值。
Matlab实现
function pi_approx = leibniz_pi(n)
% n为计算的项数
sum_val = 0;
for k = 1:n
sum_val = sum_val+(-1)^(k - 1)/(2 * k - 1);
end
pi_approx = 4 * sum_val;
end
使用示例：
% 计算10000项的莱布尼茨公式逼近的π值
approx_pi = leibniz_pi(10000);
disp(['使用莱布尼茨公式，10000项逼近的π值为：', num2str(approx_pi)]);
Python实现
def leibniz_pi(n):
sum_val = 0
for k in range(1, n + 1):
sum_val += (-1) ** (k - 1) / (2 * k - 1)
pi_approx = 4 * sum_val
return pi_approx
# 计算10000项的莱布尼茨公式逼近的π值
approx_pi = leibniz_pi(10000)
print(f"使用莱布尼茨公式，10000项逼近的π值为：{approx_pi}")
- 马青公式：
π
=
16
arctan
(
1
/
5
)
−
4
arctan
(
1
/
239
)
\pi = 16\arctan(1/5) - 4\arctan(1/239)
π=16arctan(1/5)−4arctan(1/239)
马青公式利用反正切函数的泰勒展开式，收敛速度相对较快，在早期的π计算中应用广泛。
Matlab实现
function pi_approx = machin_pi()
pi_approx = 16 * atan(1 / 5)-4 * atan(1 / 239);
end
使用示例：
approx_pi = machin_pi();
disp(['使用马青公式计算的π值为：', num2str(approx_pi)]);
Python实现
import math
def machin_pi():
pi_approx = 16 * math.atan(1 / 5)-4 * math.atan(1 / 239)
return pi_approx
approx_pi = machin_pi()
print(f"使用马青公式计算的π值为：{approx_pi}")
（二）迭代算法：高斯 - 勒让德算法
高斯 - 勒让德算法通过一系列迭代公式来逼近π的精确值，每迭代一次，有效数字大致翻倍，收敛速度极快。
迭代公式：
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 b_n}
bn+1=anbn
t
n
+
1
=
t
n
−
p
n
(
a
n
−
a
n
+
1
)
2
t_{n + 1} = t_n - p_n(a_n - a_{n + 1})^2
tn+1=tn−pn(an−an+1)2
p
n
+
1
=
2
p
n
p_{n + 1} = 2p_n
pn+1=2pn
初始值：
a
0
=
1
a_0 = 1
a0=1 ，
b
0
=
1
/
2
b_0 = 1/\sqrt{2}
b0=1/2 ，
t
0
=
1
/
4
t_0 = 1/4
t0=1/4 ，
p
0
=
1
p_0 = 1
p0=1
Matlab实现
function pi_approx = gauss_legendre_pi()
a = 1;
b = 1 / sqrt(2);
t = 1 / 4;
p = 1;
for i = 1:10
an = (a + b) / 2;
b = sqrt(a * b);
t = t - p * (a - an)^2;
a = an;
p = 2 * p;
end
pi_approx = (a + b)^2/(4 * t);
end
使用示例：
approx_pi = gauss_legendre_pi();
disp(['使用高斯-勒让德算法计算的π值为：', num2str(approx_pi)]);
Python实现
import math
def gauss_legendre_pi():
a = 1
b = 1 / math.sqrt(2)
t = 1 / 4
p = 1
for _ in range(10):
an = (a + b) / 2
b = math.sqrt(a * b)
t = t - p * (a - an) ** 2
a = an
p = 2 * p
pi_approx = (a + b) ** 2 / (4 * t)
return pi_approx
approx_pi = gauss_legendre_pi()
print(f"使用高斯-勒让德算法计算的π值为：{approx_pi}")
三、概率算法：蒙特卡洛方法
蒙特卡洛方法基于概率统计的思想。在一个边长为1的正方形及其内切圆中随机投点，根据点在圆内和正方形内的数量比例来估计π的值。如果在正方形内随机投 N N N个点，其中在圆内的点有 M M M个，根据几何概率，π近似为 4 M / N 4M/N 4M/N。
Matlab实现
function pi_approx = monte_carlo_pi(N)
% N为投点数量
x = rand(N, 1);
y = rand(N, 1);
count = sum(x.^2 + y.^2 <= 1);
pi_approx = 4 * count / N;
end
使用示例：
% 投点1000000个
approx_pi = monte_carlo_pi(1000000);
disp(['使用蒙特卡洛方法，1000000个点逼近的π值为：', num2str(approx_pi)]);
Python实现
import random
def monte_carlo_pi(N):
count = 0
for _ in range(N):
x = random.random()
y = random.random()
if x ** 2 + y ** 2 <= 1:
count += 1
pi_approx = 4 * count / N
return pi_approx
# 投点1000000个
approx_pi = monte_carlo_pi(1000000)
print(f"使用蒙特卡洛方法，1000000个点逼近的π值为：{approx_pi}")
四、拉马努金-Ramanujan计算 π \pi π公式
拉马努金（Ramanujan）提出了多个计算圆周率
π
\pi
π的公式，这些公式以其独特的数学美感和快速的收敛速度而闻名。以下是其中一个著名的公式：
[ \frac{1}{\pi} = \frac{2\sqrt{2}}{9801} \sum_{k = 0}^{\infty} \frac{(4k)!(1103 + 26390k)}{(k!)^4 396^{4k}} ]
Matlab实现
function pi_approx = ramanujan_pi()
% 初始化和
sum_val = 0;
% 计算前20项（可根据需要调整项数）
for k = 0:20
numerator = factorial(4 * k) * (1103 + 26390 * k);
denominator = (factorial(k))^4 * 396^(4 * k);
sum_val = sum_val + numerator / denominator;
end
% 根据拉马努金公式计算π
pi_approx = 1 / (sum_val * 2 * sqrt(2) / 9801);
end
使用示例：
approx_pi = ramanujan_pi();
disp(['使用拉马努金公式计算的π值为：', num2str(approx_pi)]);
在这段代码中：
- 首先初始化一个变量
sum_val用于累加每一项的值。 - 通过
for循环计算拉马努金公式中每一项的值，其中factorial函数用于计算阶乘。 - 计算完所有项的和后，根据公式
1 / (sum_val * 2 * sqrt(2) / 9801)计算出π的近似值。
Python实现
import math
def ramanujan_pi():
# 初始化和
sum_val = 0
# 计算前20项（可根据需要调整项数）
for k in range(21):
numerator = math.factorial(4 * k) * (1103 + 26390 * k)
denominator = (math.factorial(k) ** 4) * (396 ** (4 * k))
sum_val += numerator / denominator
# 根据拉马努金公式计算π
pi_approx = 1 / (sum_val * 2 * math.sqrt(2) / 9801)
return pi_approx
approx_pi = ramanujan_pi()
print(f"使用拉马努金公式计算的π值为：{approx_pi}")
在这段Python代码中：
- 同样先初始化变量
sum_val。 - 使用
for循环遍历计算每一项，math.factorial用于计算阶乘。 - 最后根据公式得出π的近似值并返回。
这两种实现方式都只是计算了拉马努金公式的前有限项来得到π的近似值，增加计算的项数可以提高近似的精度，但计算时间也会相应增加。
除上述算法外，还有以下计算圆周率的方法：
- 丘德诺夫斯基算法：这是一种基于超几何级数的算法，公式为 1 π = 12 ∑ k = 0 ∞ ( − 1 ) k ( 6 k ) ! ( 13591409 + 545140134 k ) ( 3 k ) ! ( k ! ) 3 64032 0 3 k + 3 2 \frac{1}{\pi}=12\sum_{k = 0}^{\infty}\frac{(-1)^k(6k)!(13591409 + 545140134k)}{(3k)!(k!)^3640320^{3k + \frac{3}{2}}} π1=12∑k=0∞(3k)!(k!)36403203k+23(−1)k(6k)!(13591409+545140134k)。它收敛速度极快，每计算一项可得到约14位正确小数，常用于计算超高精度的圆周率。
- BBP算法：即贝利-波尔温-普劳夫公式，公式为 π = ∑ k = 0 ∞ ( 1 1 6 k ( 4 8 k + 1 − 2 8 k + 4 − 1 8 k + 5 − 1 8 k + 6 ) ) \pi=\sum_{k = 0}^{\infty}\left(\frac{1}{16^k}\left(\frac{4}{8k + 1}-\frac{2}{8k + 4}-\frac{1}{8k + 5}-\frac{1}{8k + 6}\right)\right) π=∑k=0∞(16k1(8k+14−8k+42−8k+51−8k+61))。其独特之处在于能计算圆周率的任意第n位二进制数字，而无需计算前面的数字，在一些特定场景如验证其他算法计算的某一位正确性时很有用。
五、丘德诺夫斯基算法
Matlab实现
function pi_approx = chudnovsky_algorithm()
% 初始化变量
sum_val = 0;
% 计算前10项（可根据需要调整项数）
for k = 0:10
numerator = (-1)^k * factorial(6 * k) * (13591409 + 545140134 * k);
denominator = (factorial(3 * k)) * (factorial(k)^3) * 640320^(3 * k + 1.5);
sum_val = sum_val + numerator / denominator;
end
% 根据丘德诺夫斯基算法公式计算π
pi_approx = 1 / (12 * sum_val);
end
使用示例：
approx_pi = chudnovsky_algorithm();
disp(['使用丘德诺夫斯基算法计算的π值为：', num2str(approx_pi)]);
Python实现
import math
def chudnovsky_algorithm():
# 初始化变量
sum_val = 0
# 计算前10项（可根据需要调整项数）
for k in range(11):
numerator = (-1) ** k * math.factorial(6 * k) * (13591409 + 545140134 * k)
denominator = (math.factorial(3 * k)) * (math.factorial(k) ** 3) * 640320 ** (3 * k + 1.5)
sum_val += numerator / denominator
# 根据丘德诺夫斯基算法公式计算π
pi_approx = 1 / (12 * sum_val)
return pi_approx
approx_pi = chudnovsky_algorithm()
print(f"使用丘德诺夫斯基算法计算的π值为：{approx_pi}")
六、BBP算法
Matlab实现
function pi_approx = bbp_algorithm()
% 初始化变量
sum_val = 0;
% 计算前1000项（可根据需要调整项数）
for k = 0:1000
term1 = 4 / (8 * k + 1);
term2 = 2 / (8 * k + 4);
term3 = 1 / (8 * k + 5);
term4 = 1 / (8 * k + 6);
sum_val = sum_val+(1/16^k) * (term1 - term2 - term3 - term4);
end
% 根据BBP算法公式计算π
pi_approx = sum_val;
end
使用示例：
approx_pi = bbp_algorithm();
disp(['使用BBP算法计算的π值为：', num2str(approx_pi)]);
Python实现
def bbp_algorithm():
# 初始化变量
sum_val = 0
# 计算前1000项（可根据需要调整项数）
for k in range(1001):
term1 = 4 / (8 * k + 1)
term2 = 2 / (8 * k + 4)
term3 = 1 / (8 * k + 5)
term4 = 1 / (8 * k + 6)
sum_val += (1 / 16 ** k) * (term1 - term2 - term3 - term4)
# 根据BBP算法公式计算π
pi_approx = sum_val
return pi_approx
approx_pi = bbp_algorithm()
print(f"使用BBP算法计算的π值为：{approx_pi}")
通过对这些算法的探讨和代码实现，我们可以看到不同算法在计算π时各有优劣。古典算法展示了古人的智慧，分析算法中的无穷级数和迭代算法则体现了数学的精妙，而概率算法为我们提供了一种独特的计算思路。在实际应用中，可以根据具体需求选择合适的算法来计算π。希望这篇博文能让大家对π的计算有更深入的理解和认识。
扫一扫
2023
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
