archimedes 式 Pi 計算 - fortran66のブログ

アルキメデスは、円の内接多角形と外接多角形を書いて円周の真値を挟み込む形で円周率を求める方法を後世に残しています。ここでは内接多角形から円周率を求めることにします。

ところで、小学生の知識で周長のわかる場合としては、円は内接六角形と外接正方形に挟まれる形しかないと思うので、小学生レベルで証明できるのは、円周率は３より大きく４より小さいということだけのような気がします。しかも作図した図形から、明らかに円周率は３に近い側にあるので円周率は約３と証明できるのみかと思います。（数年前の円周率３論争が、こういう論理だとは思いませんがｗ）

定式化

直径１の円を考えることにします。

内接正方形から始める場合（より一般的には２角形というか、直径に等しい線分から出発したと考えられますが）ピタゴラスの定理を用いて少し計算すると、 $2^{n}$ 角形の円周を $s_n$ とすると、 $s_{n+1}=2^n\sqrt{2(\,1-\sqrt{1-({s_n\over 2^n})^2}\,)}$ という漸化式が得られます。
この式は見ると分かるように、桁落ちしやすい形をしているので $S_n\equiv({s_n\over 2^n})^2$ と置き換えるとともに、定石に従って ${1+\sqrt(\...)\over1+\sqrt(\...)}$ を掛けて桁落ちしない形に直します。すると $S_{n+1}={S_n / 2 \over 1+\sqrt{1-S_n}}$ という漸化式になります。ここで初期値は $S_2={1\over2}\,\,(s_2={2\sqrt2})$ となります。

なお内接六角形から始めると（より一般的には内接正三角形から始まると考えるられますが）、 $3\cdot2^n$ 角形の円周に対して漸化式は $s_{n+1}=3\cdot2^n\sqrt{2(\,1-\sqrt{1-({s_n\over 3\cdot2^n})^2}\,)}$ となります。ここで $S_n\equiv({s_n\over3\cdot2^n})^2$ と置き換えをすると $S_n$ の漸化式は上の場合と同じになります。ここで初期値は $S_1={1\over4}\,\,(s_1=3)$ となります。

ソース・コード

PROGRAM archimedes
   IMPLICIT NONE
   INTEGER, PARAMETER :: kd = 8 
   INTEGER :: i
   REAL(kd) :: pi, ss, f
   f = 4.0_kd           ! 6.0_kd
   ss = 1.0_kd / 2.0_kd ! 1.0_kd / 4.0_kd
   DO i = 3, 30
    ss = ss * 0.5_kd / ( 1.0_kd + SQRT(1.0_kd - ss) )
    f  = 2.0_kd * f 
    pi = SQRT(ss) * f
    PRINT *, i, 2**i, pi, pi - 4.0_kd * ATAN(1.0_kd)
   END DO 
   STOP
END PROGRAM archimedes
! starting from square, 2r = 1
! square f = 2 ss = 0.5, hexagon f = 6 ss = 1 / 4
! s_n+1 = 2^n SQRT( 2 * (1 - SQRT(1 - (s_n / 2^n)^2) ) ), s_2 = 2 * SQRT(2)
! S_n+1 = S_n / 2 / (1 + SQRT(1 - S_n)), S_n = (s_n / 2^n)^2, S_2 = 1 / 2

実行結果

$2^n$ 角形の場合。

見ての通り、この方法では収束は遅いようです。