はじめに
ウォリス積分(Wallis's integral)
I_n=\int_0^{\pi/2}\sin^n x\,dx
は、円周率の無限積表示(ウォリスの積公式)の導出などに登場する古典的な定積分である。$n$ を増やしていくと $I_n$ は単調に小さくなっていくが、これは0に収束するのか、それとも何らかの値に落ち着くのか。
本稿では、この問いに対して理論(漸化式・漸近評価)と数値計算(シンプソン法)の両面から答える。数値積分については既存の解説記事[1]に譲り、本稿ではシンプソン法の実装と、それを使ってウォリス積分の $n$ 依存性を実際に確認することに主眼を置く。
ウォリス積分の理論
漸化式の導出
部分積分により、$n\ge2$ に対して次の漸化式が成り立つことが知られている。
I_n=\frac{n-1}{n}I_{n-2},\qquad I_0=\frac{\pi}{2},\quad I_1=1
この漸化式から、$n$ が1増えるごとに $I_n$ は係数 $(n-1)/n<1$ 倍されるため、$I_n$ は $n$ について単調減少であることが直ちにわかる。閉じた形では、$n=2m$(偶数)のとき
I_{2m}=\frac{\pi}{2}\cdot\frac{(2m-1)!!}{(2m)!!}
$n=2m+1$(奇数)のとき
I_{2m+1}=\frac{(2m)!!}{(2m+1)!!}
と表される($!!$ は二重階乗)。
nを増やすと0に収束するか
漸化式を繰り返し用いると $I_n=\dfrac{(n-1)!!}{n!!},I_{0\text{ or }1}$ の形になり、各ステップで1未満の係数を掛け続けるので $I_n\to0$ であることは定性的にも予想できる。定量的な収束の速さは、$I_n$ をガンマ関数で表した
I_n=\frac{\sqrt{\pi}}{2}\cdot\frac{\Gamma\!\left(\frac{n+1}{2}\right)}{\Gamma\!\left(\frac{n}{2}+1\right)}
にスターリングの近似を適用することで得られ、$n\to\infty$ のとき
I_n\sim\sqrt{\frac{\pi}{2n}}
という漸近評価が成り立つ。すなわち $I_n$ は0に収束するが、その速さは指数関数的ではなく $O(1/\sqrt{n})$ という比較的緩やかなオーダーである。
nによらない不変量
漸化式 $nI_n=(n-1)I_{n-2}$ の両辺に $I_{n-1}$ を掛けると
n\,I_n\,I_{n-1}=(n-1)\,I_{n-2}\,I_{n-1}
となり、$c_n\equiv n,I_n,I_{n-1}$ とおけば $c_n=c_{n-1}$ が任意の $n\ge2$ で成り立つ。したがって $c_n$ は $n$ によらず一定であり、$c_1=1\cdot I_1 I_0=\pi/2$ より
n\,I_n\,I_{n-1}=\frac{\pi}{2}\qquad(n\ge1)
という関係が得られる。$I_n,I_{n-1}$ はともに $n\to\infty$ で0に近づくが、その積に $n$ を掛けた量は $\pi/2$ という定数に保たれる。これは $I_n\sim\sqrt{\pi/(2n)}$ という漸近評価とも整合する($n\cdot\sqrt{\tfrac{\pi}{2n}}\cdot\sqrt{\tfrac{\pi}{2(n-1)}}\to\tfrac{\pi}{2}$)。
シンプソン法による数値積分
シンプソン法は、区間を偶数個の小区間に分割し、隣接する2区間ごとに被積分関数を放物線で近似して積分する数値積分法である(詳細な導出は文献[1]を参照)。分割数を $2N$、刻み幅を $h=(b-a)/(2N)$ とすると、合成シンプソン公式は
\int_a^b f(x)\,dx\approx\frac{h}{3}\left[f(x_0)+f(x_{2N})+4\sum_{k=1}^{N}f(x_{2k-1})+2\sum_{k=1}^{N-1}f(x_{2k})\right]
で与えられる。本稿ではこれを NumPy で実装し、$f(x)=\sin^n x$ を $[0,\pi/2]$ 上で数値積分することで $I_n$ の近似値を求める。
Python実装
import numpy as np
def simpson(f, a, b, n_intervals):
"""合成シンプソン公式による数値積分(n_intervalsは偶数に丸める)"""
if n_intervals % 2 == 1:
n_intervals += 1
x = np.linspace(a, b, n_intervals + 1)
y = f(x)
h = (b - a) / n_intervals
return h / 3 * (y[0] + y[-1] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-1:2]))
def wallis_exact(n_max):
"""漸化式 I_n=(n-1)/n * I_{n-2}, I_0=π/2, I_1=1 による厳密値"""
I = np.zeros(n_max + 1)
I[0] = np.pi / 2
if n_max >= 1:
I[1] = 1.0
for n in range(2, n_max + 1):
I[n] = (n - 1) / n * I[n - 2]
return I
simpson はシンプソン法そのものを $\sin^n x$ に限らず任意の関数 $f$ に適用できる形で実装し、wallis_exact は前節の漸化式をそのままコードに落として厳密値の系列を求める。両者を $n=0,1,\dots,60$ について計算し、比較する。
n_max = 60
n_intervals = 2000
I_exact = wallis_exact(n_max)
I_simpson = np.array([
simpson(lambda x, n=n: np.sin(x) ** n, 0, np.pi / 2, n_intervals)
for n in range(n_max + 1)
])
error = np.abs(I_simpson - I_exact)
結果
シンプソン法と厳密値の比較
分割数 $2N=2000$ としてシンプソン法で求めた $I_n$ と、漸化式による厳密値を比較すると次の通りである。
| $n$ | シンプソン法 | 厳密値 | |誤差| |
|---|---|---|---|
| 0 | 1.5707963268 | 1.5707963268 | $2.2\times10^{-16}$ |
| 1 | 1.0000000000 | 1.0000000000 | $2.2\times10^{-15}$ |
| 2 | 0.7853981634 | 0.7853981634 | $1.1\times10^{-16}$ |
| 3 | 0.6666666667 | 0.6666666667 | $1.3\times10^{-14}$ |
| 5 | 0.5333333333 | 0.5333333333 | $0$ |
| 10 | 0.3865631585 | 0.3865631585 | $5.6\times10^{-17}$ |
| 20 | 0.2767696821 | 0.2767696821 | $5.6\times10^{-17}$ |
| 30 | 0.2269242244 | 0.2269242244 | $5.6\times10^{-17}$ |
| 50 | 0.1763614288 | 0.1763614288 | $1.1\times10^{-16}$ |
| 60 | 0.1611294174 | 0.1611294174 | $1.1\times10^{-16}$ |
$\sin^n x$ は $[0,\pi/2]$ 上で滑らかな関数であるため、シンプソン法(誤差のオーダーは刻み幅 $h$ に対して $O(h^4)$)は分割数2000程度で倍精度浮動小数点の丸め誤差レベル($10^{-14}$〜$10^{-17}$)まで厳密値に一致しており、両者の系列に有意な差は見られない。
厳密値(丸印)・シンプソン法(バツ印)は完全に重なっており、区別できない。破線は漸近近似 $\sqrt{\pi/(2n)}$。
図から、$I_n$ は $n$ の増加とともに単調に減少し、漸近曲線 $\sqrt{\pi/(2n)}$ にほぼ沿う形で0に近づいていくことが視覚的に確認できる。減少の速さは $1/\sqrt n$ のオーダーであり、指数関数的な減衰と比べるとかなり緩やかである(実際、$n=60$ でもまだ $I_{60}\approx0.161$ であり、0.01を下回るには $n$ をおよそ数千のオーダーまで増やす必要がある)。
シンプソン法の誤差の推移を対数スケールで見ると次のようになる。
誤差は $n$ によらずおおむね $10^{-14}$〜$10^{-17}$ の範囲に収まっており、$n$ が大きくなって被積分関数 $\sin^n x$ が $x=\pi/2$ 付近に鋭く集中するようになっても(分割数2000であれば)シンプソン法の精度が大きく劣化することはなかった。
不変量 $nI_nI_{n-1}=\pi/2$ の数値検証
理論解析で導いた不変量 $c_n=nI_nI_{n-1}$ を、厳密値の系列から実際に計算すると次のようになった。
| $n$ | $n,I_n,I_{n-1}$ |
|---|---|
| 1 | 1.5707963268 |
| 2 | 1.5707963268 |
| 5 | 1.5707963268 |
| 10 | 1.5707963268 |
| 30 | 1.5707963268 |
| 60 | 1.5707963268 |
($\pi/2\approx1.5707963268$)
すべての $n$ に対して $c_n$ は $\pi/2$ に一致しており、理論通り $n$ によらない不変量になっていることが確認できる。
$I_n$ と $I_{n-1}$ はそれぞれ単独では $n\to\infty$ で0に収束していく量だが、$n$ を掛けて積を取るとちょうど収束の遅さ($1/\sqrt n$)が打ち消し合い、一定値 $\pi/2$ に保たれる。この関係は、$n\to\infty$ の極限を取ることでウォリスの積公式
\frac{\pi}{2}=\prod_{k=1}^{\infty}\frac{2k}{2k-1}\cdot\frac{2k}{2k+1}
を導くための鍵となる関係式でもある。
まとめ
- ウォリス積分 $I_n=\int_0^{\pi/2}\sin^n x,dx$ は、部分積分による漸化式 $I_n=\frac{n-1}{n}I_{n-2}$ に従い、$n$ の増加とともに単調に減少して0に収束する。
- 収束の速さは指数関数的ではなく、漸近的に $I_n\sim\sqrt{\pi/(2n)}$ という $O(1/\sqrt n)$ のオーダーで緩やかに0へ近づく。
- シンプソン法による数値積分は、漸化式から求めた厳密値と機械精度レベル($10^{-14}$〜$10^{-17}$)で一致し、$\sin^n x$ のような滑らかな関数に対するシンプソン法の高い精度を裏付けた。
- $I_n$ 単体は0に収束する一方、$n,I_n,I_{n-1}$ という組み合わせは $n$ によらず $\pi/2$ で一定に保たれることを理論・数値の両面で確認した。この不変量は、$n\to\infty$ の極限でウォリスの積公式($\pi/2$ の無限積表示)を導く基盤にもなっている。
参考文献
付録: Pythonコード全文
"""ウォリス積分 I_n = ∫[0,π/2] sin^n(x) dx を
シンプソン法で数値積分し、漸化式による厳密値と比較する。
n を増やしたときに I_n がどのように振る舞うか
(0への収束・漸近挙動・不変量 n*I_n*I_{n-1}=π/2)も検証する。
"""
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
def simpson(f, a, b, n_intervals):
"""合成シンプソン公式による数値積分(n_intervalsは偶数に丸める)"""
if n_intervals % 2 == 1:
n_intervals += 1
x = np.linspace(a, b, n_intervals + 1)
y = f(x)
h = (b - a) / n_intervals
return h / 3 * (y[0] + y[-1] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-1:2]))
def wallis_exact(n_max):
"""漸化式 I_n=(n-1)/n * I_{n-2}, I_0=π/2, I_1=1 による厳密値"""
I = np.zeros(n_max + 1)
I[0] = np.pi / 2
if n_max >= 1:
I[1] = 1.0
for n in range(2, n_max + 1):
I[n] = (n - 1) / n * I[n - 2]
return I
n_max = 60
n_intervals = 2000 # シンプソン法の分割数
I_exact = wallis_exact(n_max)
I_simpson = np.array([
simpson(lambda x, n=n: np.sin(x) ** n, 0, np.pi / 2, n_intervals)
for n in range(n_max + 1)
])
error = np.abs(I_simpson - I_exact)
print("n, シンプソン法, 厳密値, |誤差|")
for n in [0, 1, 2, 3, 5, 10, 20, 30, 50, 60]:
print(f"{n:3d} {I_simpson[n]:.10f} {I_exact[n]:.10f} {error[n]:.3e}")
# 不変量 c_n = n*I_n*I_{n-1} の検証(理論上すべてπ/2で一定)
ns = np.arange(1, n_max + 1)
c = ns * I_exact[1:] * I_exact[:-1]
print("\nn*I_n*I_(n-1) の推移(理論上は全てπ/2で一定):")
for n in [1, 2, 5, 10, 30, 60]:
print(f"n={n:3d}: {c[n-1]:.10f} (π/2={np.pi/2:.10f})")
# --- 図1: I_n の推移と漸近曲線 ---
n_arr = np.arange(0, n_max + 1)
asymptotic = np.sqrt(np.pi / (2 * np.maximum(n_arr, 1)))
fig, ax = plt.subplots(figsize=(8, 5.5))
ax.plot(n_arr, I_exact, "o-", ms=3, color="tab:blue", label="厳密値 $I_n$(漸化式)")
ax.plot(n_arr, I_simpson, "x", ms=5, color="tab:red", label="シンプソン法(数値積分)")
ax.plot(n_arr[1:], asymptotic[1:], "--", color="gray", label=r"漸近近似 $\sqrt{\pi/(2n)}$")
ax.set_xlabel("n")
ax.set_ylabel(r"$I_n=\int_0^{\pi/2}\sin^n x\,dx$")
ax.set_title("ウォリス積分 $I_n$ の n による変化")
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
fig.savefig("ウォリス積分の収束.png", dpi=150)
# --- 図2: シンプソン法の誤差 ---
fig2, ax2 = plt.subplots(figsize=(8, 5))
ax2.semilogy(n_arr, error + 1e-20, color="tab:red")
ax2.set_xlabel("n")
ax2.set_ylabel("|シンプソン法 − 厳密値|")
ax2.set_title("シンプソン法の誤差(分割数=2000)")
ax2.grid(True, alpha=0.3)
fig2.tight_layout()
fig2.savefig("ウォリス積分の誤差.png", dpi=150)
# --- 図3: 不変量 n*I_n*I_(n-1) ---
fig3, ax3 = plt.subplots(figsize=(8, 5))
ax3.plot(ns, c, "o-", ms=3, color="tab:green")
ax3.axhline(np.pi / 2, color="tab:red", linestyle="--", label=r"$\pi/2$")
ax3.set_xlabel("n")
ax3.set_ylabel(r"$n\,I_n\,I_{n-1}$")
ax3.set_title(r"不変量 $n\,I_n\,I_{n-1}=\pi/2$ の数値検証")
ax3.legend()
ax3.grid(True, alpha=0.3)
fig3.tight_layout()
fig3.savefig("ウォリス積分の不変量.png", dpi=150)
plt.show()