fortran66のブログ

fortran について書きます。

archimedes 式 Pi 計算

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

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

定式化

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

内接正方形から始める場合(より一般的には2角形というか、直径に等しい線分から出発したと考えられますが)ピタゴラスの定理を用いて少し計算すると、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角形の場合。

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