1997 STEP I Q6 —— π 的近似分数值 | Colin Young's Personal Site

本篇的问题来自于 1997 STEP I 数学考试中的第6题，要证明的是一个很有意思的关于圆周率 \(\pi\) 近似有理分式表达式的结果。

STEP (Sixth Term Examination Paper) 是英国剑桥考试中心为本科申请者提供的数学专项考试，通常会是剑桥、帝国理工、伦敦大学、华威等名校的数学、计算机、工程等专业的录取要求之一。每年都会有一些非常精彩的考题，我在此挑选一些个人喜欢的考题作分享。

问题概要

考察积分 \(\displaystyle \int_0^1 x^4 (1-x)^4 \, \mathrm{d}x\) 以及 \(\displaystyle \int_0^1 \frac{x^4 (1-x)^4}{1+ x^2} \, \mathrm{d}x\)，证明：

\[\frac{22}{7} - \frac{1}{630} < \pi < \frac{22}{7}\]

分析与解答

出题者的设计给了很有诚意的提示。对答题者而言更多是考察算功，没有特别困难的地方。我们这里也就简单给出推导的关键步骤。

首先

\[\begin{aligned} \int_0^1 x^4 (1-x)^4 \, \mathrm{d}x &= \int_0^1 \left(x^8−4x^7+6x^6−4x^5+x^4 \right)\, \mathrm{d}x\\ & = \frac{1}{9} - \frac{4}{8} + \frac{6}{7} - \frac{4}{6} + \frac{1}{5} \\ & = \frac{1}{630} \end{aligned}\]

利用长除法（long division）或者其他技巧，可以验证

\[x^4 (1-x)^4 = (x^6−4x^5+5x^4−4x^2+4)(x^2+1) - 4\]

于是

\[\begin{aligned} \int_0^1 \frac{x^4 (1-x)^4}{1+ x^2} \, \mathrm{d}x &= \int_0^1 \left(x^6−4x^5+5x^4−4x^2+4 - \frac{4}{x^2+1} \right)\, \mathrm{d}x\\ & = \frac{1}{7} - \frac{4}{6} + 1 - \frac{4}{3} + 4 - 4\arctan1 \\ & = \frac{22}{7} - \pi \end{aligned}\]

在积分区间内，被积函数 \(\frac{x^4 (1-x)^4}{1+ x^2} \geq 0\)，因此有

\[\int_0^1 \frac{x^4 (1-x)^4}{1+ x^2} \, \mathrm{d}x > 0\]

另一方面，对于 \(0\leq x \leq 1\)，有 \(x^2+1 \geq 1\)，于是 \(\frac{x^4 (1-x)^4}{1+ x^2} \leq x^4 (1-x)^4\)，因此又有

\[\int_0^1 \frac{x^4 (1-x)^4}{1+ x^2} \, \mathrm{d}x < \int_0^1 x^4 (1-x)^4 \, \mathrm{d}x\]

结合先前得到的积分结果，可建立不等式关系

\[0 < \frac{22}{7} - \pi < \frac{1}{630}\]

整理得到

\[\frac{22}{7} - \frac{1}{630} < \pi < \frac{22}{7}\]

即可以取 \(\frac{22}{7}\) 作为 \(\pi\) 的近似值，其误差不超过 \(\frac{1}{630}\)。

一些讨论

由这个问题得到启发，我们可以设想如下的思路来得到一系列 \(\pi\) 的近似分式。我们先构造一个辅助多项式函数 \(f(x)\)，通过对它的积分的数值结果的一些估计，得出一个近似等于 \(\pi\) 的有理分式。这个多项式函数 \(f(x)\) 需要具有的性质如下：

在 \(0 \leq x \leq 1\) 的区间内，\(0 \leq f(x) \ll 1\)，如此便有 \(\displaystyle \int_0^1 f(x) \, \mathrm{d}x \approx 0\)，此积分越接近 \(0\)，后面得到的 \(\pi\) 的近似误差就越小
为了使得 \(\displaystyle \int_0^1 \frac{f(x)}{x^2+1} \, \mathrm{d}x\) 能搞出一个带 \(\pi\) 的积分结果，我们要求 \(f(x)\) 除以 \(x^2+1\) 后的余项，只能包含一个常数项 \(c\)，不能含有 \(x\) 的一次项，即 \(f(x)\) 可以被写成 \(f(x)=g(x)(x^2+1) + c\) 的形式，其中 \(g(x)\) 也是一个多项式

如此，在计算 \(\displaystyle \int_0^1 \frac{f(x)}{x^2+1} \, \mathrm{d}x\) 时，我们可以将其写成 \(\displaystyle \int_0^1 \left( g(x) + \frac{c}{x^2+1} \right) \, \mathrm{d}x\)。其中 \(\displaystyle \int_0^1 g(x)\, \mathrm{d}x\) 的积分会给出一个有理分式，\(\displaystyle \int_0^1 \frac{c}{x^2+1}\, \mathrm{d}x\) 则给出一个含有 \(\pi\) 的分式。我们不妨记

\[\int_0^1 \frac{f(x)}{x^2+1} \, \mathrm{d}x = A + B\pi\]

我们将很快看到，这里的 \(A\) 和 \(B\) 会是一组互为异号的有理分式。注意到被积函数 \(f(x)\approx 0\) 的要求，上面的积分结果应为一个也很接近于 \(0\) 的数，于是

\[\pi \approx -\frac{A}{B} = \left| \frac{A}{B} \right|\]

故可以将这里的 \(\left\vert \frac{A}{B} \right\vert\) 作为 \(\pi\) 的近似式。

接下来分析下这个近似的误差如何。\(\displaystyle \int_0^1 f(x) \, \mathrm{d}x\) 的结果显然也是一个有理分式，记作

\[\int_0^1 \frac{f(x)}{x^2+1} \, \mathrm{d}x = C\]

对于 \(0\leq x\leq1\)，有

\[1 \leq x^2+1 \leq 2\]

则

\[\frac{1}{2} f(x) \leq \frac{f(x)}{x^2+1} \leq f(x)\]

于是

\[\frac{1}{2}\int_0^1 f(x) \, \mathrm{d}x < \int_0^1 \frac{f(x)}{x^2+1}\, \mathrm{d}x < \int_0^1 f(x) \, \mathrm{d}x\]

即

\[\frac{1}{2}C < A + B\pi < C\]

考虑到 \(A\) 和 \(B\) 可能的正负号，可以得到

\[\begin{aligned} \frac{1}{2} \left|\frac{C}{B}\right| < \pi - \left| \frac{A}{B} \right| < \left|\frac{C}{B}\right| \qquad \text{if } B>0 \\ \frac{1}{2} \left|\frac{C}{B}\right| < \left| \frac{A}{B} \right| - \pi < \left|\frac{C}{B}\right| \qquad \text{if } B<0 \end{aligned}\]

如果构造的 \(f(x)\) 在积分区间内确实满足 \(\displaystyle \int_0^1 f(x) \, \mathrm{d}x \approx 0\)，则 \(C\) 会是一个很小的正数，于是 \(\displaystyle \pi \approx \left\vert \frac{A}{B} \right\vert\) 的误差将不小于 \(\displaystyle \frac{1}{2 }\left\vert \frac{C}{B} \right\vert\) 但也不会超过 \(\displaystyle \left\vert \frac{C}{B} \right\vert\)。粗暴点将 \(\displaystyle \left\vert \frac{C}{B} \right\vert\) 视为 \(\displaystyle \pi \approx \left\vert \frac{A}{B} \right\vert\) 的误差也大差不差了。

以之前的计算为例，取 \(f(x) = x^4 (1-x)^4\)，它在 \(0\leq x\leq1\) 的区间内确实取值都几乎为 \(0\)，而且计算 \(f(x)\) 除以 \(x^2+1\) 之后的余项也都符合我们对辅助多项式的要求。而我们已经算出

\[A = \frac{22}{7}, \quad B=-1, \quad C=\frac{1}{630}\]

得到了一个误差略小于 \(\frac{1}{630}\) 的结论 \(\pi \approx \frac{22}{7}\)。

如果取更强的辅助函数，我们还可以得到精度更好的 \(\pi\) 的近似值。

照葫芦画瓢，可以验证形如 \(f(x) = x^{4n} (1-x)^{4n}\) 形式的多项式都符合所需的条件。

比如我们试着取 \(f(x) = x^8 (1-x)^8\)，可以算出：

\[\begin{aligned} \int_0^1 \frac{x^8 (1-x)^8}{x^2+1}\, \mathrm{d}x &= -\frac{188684}{15015} + 4\pi \\ \int_0^1 x^8 (1-x)^8 \, \mathrm{d}x &=\frac{1}{218790} \end{aligned}\]

这样得到的近似为：

\[\pi \approx \frac{47171}{15015} = {\color{blue}{3.14159}}17\cdots\]

误差小于 \(\frac{1}{875160} \approx 10^{-6}\)，相当不错了。

最后大力出奇迹一波（当然得借助计算机了），让我们试着取 \(f(x) = x^{20} (1-x)^{20}\)，可以验证：

\[\begin{aligned} \int_0^1 \frac{x^{20} (1-x)^{20}}{x^2+1}\, \mathrm{d}x &= \frac{26856502742629699}{33393321606645} - 256\pi \\ \int_0^1 x^{20} (1-x)^{20} \, \mathrm{d}x &=\frac{1}{5651707681620 } \end{aligned}\]

这样得到的近似为：

\[\pi \approx \frac{26856502742629699}{8548690331301120} = {\color{blue}{3.141592653589793}} 79\cdots\]

误差不超过 \(\frac{1}{1446837166494720} \approx 10^{-16}\)。类比下，用这个近似分式算个地球周长啥的，误差也只是在原子半径的尺度了，嗯，洒家这辈子值了！

几个彩蛋

也许有的读者会好奇，这东西真的有人会愿意耐着性子算？

这里安利一个最近发现的宝藏工具——SymPy，一个支持符号运算的 Python library，实则验（zùo）证（bì）计（tōu）算（lǎn）的利器。

糊了一段能跑的屎山代码：

from sympy import *
import pandas as pd
import math

def display_pis(n=3, precision=10):
    orders = []
    apprx_pis = []
    decimals = []
    max_errors = []
    true_errors = []
    
    for i in range(1,n):
        x = Symbol('x')
        n = 4*i
        
        expr = integrate(x**n * (1-x)**n / (x**2 + 1), (x, 0, 1))
        expr_pi_coeff = expr.as_coefficients_dict()[pi]
        expr_fraction = expr.as_coefficients_dict()[1]
        appr_pi = simplify(-expr_fraction / expr_pi_coeff)
        
        max_error = abs(integrate(x**n * (1-x)**n, (x, 0, 1)) / 
                        expr.as_coefficients_dict()[pi])
        
        orders.append(n)
        apprx_pis.append(appr_pi)
        decimals.append(float(appr_pi))
        max_errors.append(max_error)
        true_errors.append(float(abs(math.pi - appr_pi)))
        
    df_values = pd.DataFrame({'order': orders, 
                              'rational approx. for π': apprx_pis, 
                              'decimal repr.': decimals,
                              'max error est.': max_errors, 
                              'true error': true_errors
                             })
    df_values.index = df_values.index + 1
    pd.set_option("display.precision", precision)
    df_values.style.format({'true error': "{:.2E}".format})
    return df_values

display_pis(n=6, precision=15)

当然，还有更会玩的大佬。我在查素材时，无意发现了 Rutgers 有个数学家，一口气算了 1000 个这玩意儿吧，有兴趣地可以点下面的传送门围观：

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/pi14.txt