數學傳播

本文謹獻給作者的恩師 Darío Castellanos 教授(生於1937年12月4日;逝於1995年11月23日)。 作者永遠懷念他的教誨。

要獲得近似圓周率的公式並不是困難的事情。舉例來說, 我們已有下列級數, $$ \sin x = x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \cdots, $$ $$ \cos x = 1 - \frac{x^{2}}{2} + \frac{x^{4}}{24} - \cdots, $$ $$ \frac{1}{2} \sin {2x} = \sin {x} \cos {x} = x - \frac{2}{3}x^{3} + \frac{2}{15}x^{5} - \cdots $$ 如果只考慮前兩式, 然後忽略五次或五次以上的項, 我們便可得到 $$ 3 \sin {x} - x \cos {x} \approx 2x, $$ 因而 \begin{equation} %(1) x \approx \frac{3\sin {x}}{2+\cos {x}}, \end{equation} 這是首先由Nicolaus Cusanus給出的公式; 以後荷蘭的數學家兼物理學家Willebrord Snellius 也獨立地發現了此公式。 Snellius亦即發現光學上著名的折射原理(Snell's law)的Snell。 兩個名字的差異起因於拉丁文和英文的拼法不同。

現在,如果考慮所有三式, 然後忽略七次或七次以上的項, 我們便可得到 $$ 14 \sin {x} - 6x \cos {x} + \sin {x} \cos {x} \approx 9x $$ 或 \begin{equation} %(2) x \approx \sin {x} \frac{14+\cos{x}}{9+6\cos{x}}, \end{equation} 這是首先由數學巨人Isaac Newton所發現。

由這兩條公式所得到的近似值會越來越好, 如果所取的角度越來越小。若取角度為 $\displaystyle\frac{\pi}{12}$, 公式(1)給了圓周率值 $\pi\approx 3.141509994$, 但是Newton的公式給了 $\pi\approx 3.141592169$。 當然,真正的數值是3.141592653589793238...

在此得到(1)和(2)的方法是由Newton所導出 。現在就讓我們遵循著大師的路徑來擴展他的過程。

首先考慮下列方程式 \begin{equation} %(3) x \approx A_{1} \sin {x} - A_{2} x \cos {x} + A_{3} \sin {2x} - A_{4} x \cos {2x} + \cdots + A_{2s-1} \sin {sx} - A_{2s} x \cos {sx}, \end{equation} 在此 $A_{k}$ 為有待決定的常數。接著利用適當的Maclaurin級數展開右手邊的項,我們可以發現下列系統 $$ \begin{array}{rcrcrcrcrcrcrcr} A_{1} & - & A_{2} & + & 2A_{3} & - & A_{4} & + & \cdots & + & sA_{2s-1} & - & A_{2s} & = & 1 \\ A_{1} & - & 3A_{2} & + & 2^{3}A_{3} & - & 3\cdot 2^{2}A_{4} & + & \cdots & + & s^{3}A_{2s-1} & - & 3s^{2}A_{2s} & = & 0 \\ & & & & & & \cdots\cdots & & & & & & & & \\ A_{1} & - & (4s\!-\!1)A_{2} & + & 2^{4s\!-\!1}A_{3} & - & (4s\!-\!1)2^{4s\!-\!2}A_{4} & + & \cdots & + & s^{4s\!-\!1}A_{2s\!-\!1} & - & (4s\!-\!1)s^{4s\!-\!2}A_{2s} & = & 0 \end{array} $$ 有 $2s$ 條方程式及 $2s$ 未知數可用來解 $A_{k}$。解(3)中的 $x$ 我們發現 $$ x \approx \frac{A_{1}\sin{x}+A_{3}\sin{2x}+A_{5}\sin{3x}+\cdots+A_{2s-1}\sin {sx}} {1+A_{2}\cos{x}+A_{4}\cos{2x}+A_{6}\cos{3x}+\cdots+A_{2s}\cos {sx}}. $$ 取 $s=1$, 則 $A_{1}=\displaystyle\frac{3}{2}$, $A_{2}=\displaystyle\frac{1}{2}$, 我們得到了公式(1)。

為了列出其他 $s$ 值所對應的公式,我們利用加法公式 $\sin {2x}=2\sin {x}\cos {x}$, $\cos {2x}=\cos^{2}x-\sin^{2}x$, 等等, 以避免三角函數多倍角度上複雜的計算。 取 $s=2$, 得到公式 \begin{equation} %(4) x \approx \frac{5}{3}\sin {x}\frac{16+5\cos {x}}{17+16\cos {x} + 2\cos^{2}x}, \end{equation} 取 $s=3$, 得到公式 \begin{equation} %(5) x \approx \frac{7}{5}\sin {x}\frac{92+66\cos {x}+7\cos^{2}x} {82+111\cos {x}+36\cos^{2}x+2\cos^{3}x}, \end{equation} 取 $s=4$, 得到公式 \begin{equation} %(6) x \approx \frac{1}{35}\sin {x}\frac{91648+103511\cos {x}+28544\cos^{2}x+1522\cos^{3}x} {1667+2944\cos {x}+1560\cos^{2}x+256\cos^{3}x+8\cos^{4}x}. \end{equation} 讓 $x=\displaystyle\frac{\pi}{6}$, 我們得到一些有趣的圓周率近似值。公式(4), (5)和(6)分別給出 \begin{eqnarray} 5 \Big (\frac{32+5\sqrt{3}}{37+16\sqrt{3}}\Big ) &=& 3.141592229\ldots, \\ %(7) \frac{21}{5} \Big (\frac{389+132\sqrt{3}}{436+225\sqrt{3}}\Big ) &=& 3.14159265346\ldots, \\ %(8) \frac{3}{70} \Big (\frac{45, 224+209, 305\sqrt{3}}{5, 683+3, 136\sqrt{3}}\Big ) &=& 3.141592653589754\ldots; %(9) \end{eqnarray} 再說一次,真正的數值當然是 3.141592653589793238...

這裡有一個有趣的地方值得注意。(8)和(9)改善了下列Ramanujan的公式 \begin{equation} %(10) \frac{63}{25} \Big (\frac{17+15\sqrt{5}}{7+15\sqrt{5}}\Big ) = 3.14159265380\ldots \end{equation} 雖然以相似大小的係數觀點來看, (10)確實比(7)來得近似圓周率, 但是我們要強調的是Ramanujan導出公式(10)所使用的方法是奠基於橢圓模函數的理論 (elliptic modular functions), 它所需要考慮到的和牽連到的內容和事情要比我們在此使用的簡單過程要複雜和困難得多太多了。

這裡還有一個問題尚待回答。我們如何知道上述系統經常有解? 也就是說, 上述系統的係數矩陣經常是非奇異的(non-singular)? 這個問題的答案很有趣, 但是卻不太明顯。它等同於去找到 $n$ 個函數使得它們的線性組合 $k(x)$ 和某一被給的函數 $f(x)$ 在下列的觀點下一致: $f^{(j)}(0)=k^{(j)}(0)$, $j=0, 1, \ldots, n-1$。

為了證明這是經常可能的, 讓我們考慮一具有 $n$ 個獨立解的 $n$ 階齊次線性微分方程。 讓 $k(x)$ 為此 $n$ 個函數的線性組合, 且具有下列的初始條件: $f^{(j)}(0)\!=\!k^{(j)}(0)$, $j\!=\!0, 1, \ldots, n-1$。 由於存在和唯一定理, 我們知道如此的 $k$ 值是存在和唯一的。 更甚的是, 若 $n=2s$ 且 $s$ 為此微分方程獨立解的個數且為奇數, 則假使 $f$ 為奇函數(偶函數), 它便可以被解釋為 $s$ 個奇函數(偶函數)的線性組合。 特別地, 下列微分方程 $$ (D^{2}+1^{2})^{2} (D^{2}+2^{2})^{2} \cdots (D^{2}+s^{2})^{2}y = 0 $$ 為 $4s$ 階而且有 $2s$ 個奇函數解如 $\sin {x}$, $x\cos {x}$, $\sin {2x}$, $x\cos {2x}$, $\ldots$, $\sin {sx}$, $x\cos {sx}$.

因此,選擇奇函數 $y=x$ 我們可以寫成 $$ x \approx A_{1}\sin {x} - A_{2}x\cos {x} + A_{3}\sin {2x} - A_{4}x\cos {2x} + \cdots + A_{2s-1}\sin {sx} - A_{2s}x\cos {sx}, $$ 此即(3)。

但是, 如果我們用矩陣來表達這個系統並且用行列式來證明這個矩陣經常有唯一解, 則可能會更加有趣。定義一 $(2s)\times(2s)$ 行列式 $$ H_{s}(u_{1}, \ldots, u_{s}) = \left | \begin{array}{ccccccccccccc} u_{1} & & -1 & & u_{2} & & -1 & & \cdots & & u_{s} & & -1 \\ u_{1}^{3} & & -3u_{1}^{2} & & u_{2}^{3} & & -3u_{2}^{2} & & \cdots & & u_{s}^{3} & & -3u_{s}^{2} \\ u_{1}^{5} & & -5u_{1}^{4} & & u_{2}^{5} & & -5u_{2}^{4} & & \cdots & & u_{s}^{5} & & -5u_{s}^{4} \\ \vdots & & \vdots & & \vdots & & \vdots & & \vdots & & \vdots & & \vdots \\ u_{1}^{4s-1} & & (1-4s)u_{1}^{4s-2} & & u_{2}^{4s-1} & & (1-4s)u_{2}^{4s-2} & & \cdots & & u_{s}^{4s-1} & & (1-4s)u_{s}^{4s-2} \end{array} \right| $$ 我們可以證明 $$ H_{s}(u_{1}, \ldots, u_{s}) = -2u_{1}^{3} (u_{2}^{2}-u_{1}^{2})^{4} \cdots (u_{s}^{2}-u_{1}^{2})^{4} H_{s-1}(u_{2}, \ldots, u_{s}), $$ 而且 $$ H_{s}(u_{1}, \ldots, u_{s}) = (-2)^{s} (u_{1}\cdots u_{s})^{3} \prod_{1\le i\lt j\le s}(u_{j}^{2} - u_{i}^{2})^{4}. $$ 對所有的 $j$ 讓 $u_{j}=j$, 我們看到了這個係數矩陣是非奇異的。 此種型式的行列式我們叫它 confluent Vandermonde 行列式。

現在讓我們考慮下列方程式 \begin{equation} %(11) x \approx A_{1}\sin {x} - A_{2}\sin {2x} + A_{3}\sin {3x} - \cdots + (-1)^{s+1}A_{s}\sin {sx} \end{equation} 在等號右手邊的項, 使用適當的正弦級數經過Maclaurin級數展開後,我們可以發現下列系統 $$ \begin{array}{rcrcrcrcrcr} A_{1} & - & 2A_{2} & + & 3A_{3} & - & \cdots & + & (-1)^{s-1}sA_{s} & = & 1 \\ A_{1} & - & 2^{3}A_{2} & + & 3^{3}A_{3} & - & \cdots & + & (-1)^{s-1}s^{3}A_{s} & = & 0 \\ {\vdots} & & {\vdots} & & {\vdots} & & {\vdots} & & {\vdots} & & {\vdots} \\ A_{1} & - & 2^{2s-1}A_{2} & + & 3^{2s-1}A_{3} & - & \cdots & + & (-1)^{s-1}s^{2s-1}A_{s} & = & 0 \end{array} $$ 有 $2s$ 條方程式及 $2s$ 未知數可用來解 $A_{k}$。

欲證明這個係數矩陣是非奇異的, 我們可以考慮微分方程 $(D^{2}+1^{2})(D^{2}+2^{2})\cdots(D^{2}+s^{2})y=0$。 它有 $2s$ 階並且有 $s$ 個奇函數解如 $\sin {x}$, $\sin {2x}$, $\cdots$, $\sin {sx}$。

可是, 就如前面已敘述過的一樣, 我們可以使用行列式來證明這個矩陣經常有唯一解。定義 $$ G(v_{1}, \ldots, v_{s}) = \left | \begin{array}{ccccc} v_{1} & v_{2} & v_{3} & \cdots&v_{s} \\ v_{1}^{3} & v_{2}^{3} & v_{3}^{3} & \cdots & v_{s}^{3} \\ v_{1}^{5} & v_{2}^{5} & v_{3}^{5} & \cdots & v_{s}^{5} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ v_{1}^{2s-1} & v_{2}^{2s-1} & v_{3}^{2s-1} & \cdots & v_{s}^{2s-1} \end{array} \right | $$ 我們可以證明 $$ G(v_{1}, \ldots, v_{s}) = v_{1} \cdots v_{s} \prod_{1\le i\lt j\le s}(v_{j}^{2}-v_{i}^{2}). $$ 由於這個行列式本質上也是屬於Vandermonde行列式, 而且對所有的 $j$ 讓 $v_{j}=(-1)^{j-1}j$, 我們看到了這個係數矩陣是非奇異的。

現在假設 $P$ 是近似於圓周率的一個數值。讓 $x=P-\pi$ 表示近似值和真實值之間的誤差, 同時注意到 $\sin {(k(P-\pi))}=(-1)^{k}\sin {kP}$。因此, 在(11)式中,我們用 $P-\pi$ 取代 $x$, 得到

\begin{eqnarray*} %(11.1) P - \pi \approx \sum_{k=1}^{s} (-1)^{k-1}A_{k}\sin (k(P-\pi)).\tag*{(11.1)} \end{eqnarray*} 緊接著, 使用 $\sin {(k(P-\pi))}=(-1)^{k}\sin {kP}$ 在(11.1)中, 我們得到

\begin{equation} %(11.3) P - \pi \approx -\sum_{k=1}^{s} A_{k}\sin {kP},\tag*{(11.3)} \end{equation} 經過移項後, 我們得到了

\begin{equation} %(11.5) \pi\approx P+A_{1}\sin {P}+A_{2}\sin {2P}+A_{3}\sin {3P}+\cdots+A_{s}\sin {sP}.\tag*{(11.5)} \end{equation} 解上面的系統後並令 $s=1$ 我們發現

\begin{equation} %(12) \pi \approx P + \sin {P}, \end{equation}

這是由D. Shanks所給出的。讀者知道他是誰嗎? 他在1961年的七月和作者的恩師J. W. Wrench, Jr., 利用位於紐約IBM資料處理中心的一台超級電腦IBM 7090, 以C. Störmer發現的公式 $$ \frac{\pi}{4} = 6\tan^{-1} \Big (\frac{1}{8}\Big ) + 2\tan^{-1}\Big (\frac{1}{57}\Big ) + \tan^{-1} \Big (\frac{1}{239}\Big ). $$ 作為計算主程式, 而用Gauss發現的公式 $$ \frac{\pi}{4} = 12\tan^{-1} \Big (\frac{1}{18}\Big ) + 8\tan^{-1} \Big (\frac{1}{57}\Big ) - 5\tan^{-1} \Big (\frac{1}{239}\Big ). $$ 作為核對副程式, 使用8小時43分鐘的時間計算了100,265位圓周率。因此, 間接來說, 他也是作者的師祖輩呢, 哈哈。

假設 $P$ 為圓周率的 $n$ 位小數近似值, 則(12)式將可得到 $3n$ 位小數近似值。 這是因為我們注意到(12)式也可以寫成 $\pi+x-(x-\displaystyle\frac{x^{3}}{3!}-\frac{x^{5}}{5!}-\cdots) =\pi+\frac{x^3}{3!}-\frac{x^5}{5!}+\cdots$, 則誤差 $x$ 可以減少至 $\displaystyle\frac{x^{3}}{3!}$。

解上面的系統後並令 $s$ 為更大的正整數,則我們就可得到更佳的近似值: \begin{equation} %(13) \pi \approx P + \frac{4}{3}\sin {P} + \frac{1}{6}\sin {2P} \end{equation} 得到 $5n$ 位小數, \begin{equation} %(14) \pi \approx P + \frac{3}{2}\sin {P} + \frac{3}{10}\sin {2P} + \frac{1}{30}\sin {3P} \end{equation} 得到 $7n$ 位小數, \begin{equation} %(15) \pi \approx P + \frac{8}{5}\sin {P} + \frac{2}{5}\sin {2P} + \frac{8}{105}\sin {3P} + \frac{1}{140}\sin {4P} \end{equation} 得到9n位小數, 等等。例如, 讓 $P=3.1$, (12)給了 $\pi=3.1415806$, (13)給了 $\pi=3.1415926494$, (14)給了 $\pi=3.141592653588$, (15)給了 $\pi=3.1415926535897922$。

若將上面所敘述過的概念和下列的Fourier展開作一比較, 我們會發現另一個有趣之處: $$ x = 2\sum_{n=1}^\infty \frac{(-1)^{n-1}\sin {nx}}{n}, \qquad -\pi \lt x \lt \pi. $$ 用 $P-\pi$ 取代 $x$, 而且如同前面所敘述過的, 利用 $\sin(k(P-\pi))=(-1)^{k}\sin {kP}$, 得到 \begin{equation} %(16) \pi = P + 2\sum_{n=1}^\infty \frac{\sin {nP}}{n}, \qquad 0 \lt P \lt 2\pi. \end{equation}

在研究過Fourier級數後, 我們知道, 當一個系統 $S=\{\phi_{0}, \phi_{1}, \phi_{2}, \ldots\}$ 在 $[a, b]$ 上為單範正交(orthonormal)時, 則在此系統的元素所有可能線性組合之中, 若以最小均方誤差的觀點來看, 會有一個在 $[a, b]$ 上為黎曼可積分的(Riemann-integrable)函數 $f$ 其Fourier級數的部份和 $s_{n}$ 將會產生對於 $f$ 的最佳可能近似值。 從(13)式到(16)式都是由正弦函數所線性組合而成的, 而且對於奇函數 $f(x)=\pi$, $0\lt x\lt \pi$, 這些周期為 $2\pi$ 的正弦函數會產生越來越好的近似值。

作者在此希望對其恩師, 臺灣大學數學系程舜仁教授, 在Vandermonde行列式方面的教導表達內心由衷的感謝。

---本文作者現就讀於國立臺灣大學數學研究所碩士班三年級---