アルキメデス法による円周率の算出について

−円に接する正多角形の辺数と近似円周率の正しい桁数との関係 -

田中一夫 2002.3.25.

tanaka.kazuo205@canon.co.jp, User494704@aol.com

要旨

円に内接ない外接する正 [Maple Math] 角形の周囲長より求められる近似円周率はそれぞれ [Maple Math] ,

[Maple Math] 桁正しい。

また、内接周囲長と外接周囲長とを1:2に内分した値は [Maple Math] 桁正しい近似円周率を与える。

解析及び数値データを与える

1 . はじめに

円の直径と周囲長との比で定義される円周率の算出は、Archimedesが円に接する正多角形の周長を

求めることにより始まった。円周率の歴史的な経緯及び逸話は参考文献を参照されたい。

本論は、円に接する正多角形の周囲長より近似円周率を求める方法にて、多角形の辺数と得られる

近似円周率の正しい桁数との関係を議論するものである。

2.解析

> restart:

> readlib(hypergeom):
with(linalg):

Warning, new definition for norm

Warning, new definition for trace

2.1内接円による円周率の算出

直径1の円に内接する辺数 [Maple Math]  なる正多角形を考える。その周囲長は円周率の近似値であり   

> pi[In,N] := N*sin(Pi/N);

[Maple Math]

で与えられる。このときの誤差は

> Error[In,N]:=pi-pi[In,N];

[Maple Math]

であり、これをTaylor展開をすると

> f(x):=sin(x):
taylor(f(x),x,12):
poly:=convert(%,polynom):
Pi-N*subs(x=Pi/N,poly):
Error_Series[In,N]:=expand(%);

[Maple Math]

であり、その第k項は

> Term[In,N]:=(-1)^(k+1)*Pi^(2*k+1)*N^(-2*k)/(2*k+1)!;

[Maple Math]

となり、これを用いて誤差は超幾何関数で表記できる。

> Error[In,N]:=value(Sum(Term[In,N],k=1..infinity));

[Maple Math]

Taylor展開の各項は次第に小さくなるが、第1項と第2項との比を求めると

> value(subs(k=2,Term[In,N])/subs(k=1,Term[In,N]));

[Maple Math]

であるた [Maple Math] であれば誤差のTaylor級数にて、支配的なのは第1項目のみであり、第2項以降は
近似円周率の正答桁数には影響しない。その第1項は

> First_Term[In,N]:=subs(k=1,Term[In,N]);
First_Term[In,n]:=subs(N=10^n,%);

[Maple Math]

[Maple Math]

であり、この常用対数を求めると

> log[10](First_Term[In,N]):
assume(N,positive):
a:=expand(%):
Digit[In,N]:=evalf(op(1,a)+op(2,a))+op(3,a);
log[10](First_Term[In,n]):
Digit[In,n]:=evalf(simplify(expand(%)));

[Maple Math]

[Maple Math]

である。従って、内接正 [Maple Math] 角形の周囲長より求められる円周率は [Maple Math] 桁正しいことを得る。

2.2.外接円による円周率の算出

> restart:

次に、直径1の円に外接する辺数 [Maple Math]  なる正多角形を考える。その周囲長も円周率の近似値であり

> pi[Out,N]:=N*tan(Pi/N);

[Maple Math]

で与えられる。このときの誤差は

> Error[Out,N]:=Pi-pi[Out,N];

[Maple Math]

であり、このTaylor展開は一般的に

> f(x):=tan(x):
taylor(f(x),x,12):
poly:=convert(%,polynom):
N*subs(x=Pi/N,poly)-Pi:
Error_Series[Out,N]:=expand(%);

[Maple Math]

これはBernoulli数を用いて

> Term[Out,N]:=2^(2*k)*(2^(2*k)-1)/(2*k)!*abs(bernoulli(2*k))*Pi^(2*k-1)*N^(-2*k+2):
Error[Series,Out,N]:=value(Sum(Term[Out,N],k=2..infinity));

[Maple Math]

となる。内接多角形の場合と同様に、正しい円周率の桁を決定するのは第1項

> First_Term[Out,N]:=value(Sum(Term[Out,N],k=2..2));
First_Term[Out,n]:=subs(N=10^n,%);

[Maple Math]

[Maple Math]

のみである。これの常用対数を求めると

> log[10](First_Term[Out,N]):
assume(N,positive):
a:=expand(%):
Digit[Out,N]:=evalf(op(1,a)+op(2,a))+op(3,a);
Digit[Out,n]:=simplify(subs(N=10^n,%));

[Maple Math]

[Maple Math]

となり、内接多角形の場合と同様に、外接正 [Maple Math] 角形の周囲長より求められる円周率も   [Maple Math]

桁正しいことを得る。

2.3.内接円、外接円の内分による円周率の算出

> restart:

内接ないし外接正 [Maple Math] 角形の周囲長より得られる近似円周率の誤差はそれぞれ [Maple Math] 及び [Maple Math] であるため、

両者の1:2 なる内分値

> pi[In]:=N*sin(Pi/N):
pi[Out]:=N*tan(Pi/N):
pi[Interpolation]:=(2/3)*pi[In]+(1/3)*pi[Out];

[Maple Math]

はより良い近似円周率が得られることが期待できる。ここで、この内分値の近似円周率の正しい桁数を検討する。

この近似値の誤差は

> Error[Interpolation,N]:=%-Pi;
Error[Interpolation,n]:=subs(N=10^n,%):

[Maple Math]

であり、Taylor展開すると

> f(x):=(2/3)*sin(x)+(1/3)*tan(x):
taylor(f(x),x,15):
poly:=convert(%,polynom):
Error[Series,Interpolation,N]:=expand(N*subs(x=Pi/N,poly)-Pi);

[Maple Math]

となり、内接ないし外接多角形の場合と同じく、近似円周率の正しい桁数はこの級数の第1項

> First_Term[Interpolation,N]:=op(1,Error[Series,Interpolation,N]);
First_Term[Interpolation,n]:=subs(N=10^n,%);

[Maple Math]

[Maple Math]

を検討すれば良い。これの常用対数は

> log[10](First_Term[Interpolation,N]):
assume(N,positive):
a:=expand(%):
Digit[Interpolation,N]:=evalf(op(1,a)+op(2,a))+op(3,a);
log[10](First_Term[Interpolation,n]):
b:=(simplify(expand(%))):
Digit[Interpolation,n]:=
evalf((op(1,numer(b))+op(2,numer(b))+op(3,numer(b)))/denom(b))
+simplify((op(4,numer(b))+op(5,numer(b)))/denom(b));

[Maple Math]

[Maple Math]

であり、内接、外接正 [Maple Math] 角形の周囲長を 1:2 に内分する近似円周率は   [Maple Math]

正しいく内接、外接正 [Maple Math] 角形から求めた値の倍の桁が得られる。

3.数値データ

ここで、数値データを示す。

> restart:

3.1.Archimedesの場合

Archimedesは内接、外接する正96角形の周囲長を求め 221/71<Pi<22/7なる不等式を得た

そして、円周率を 3.14とした。この場合を評価する。

> N:=96;
n:=(evalf(log[10](N)));
Digit[Theoretical]:=2*%;
Digits:=4*round(%):
pi[Exact_Value]:=evalf(Pi);
pi[In]:=evalf(N*sin(Pi/N));
Digit[In]:=1-log[10](evalf(pi[Exact_Value]-%));
pi[Out]:=evalf(N*tan(Pi/N));
Digit[Out]:=1-log[10](evalf(%-pi[Exact_Value]));
pi[Interpolation]:=evalf((2/3)*pi[In]+(1/3)*pi[Out]);
Digit[Interpolation]:=1-log[10](evalf(%-pi[Exact_Value]));

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

[Maple Math]

3.2.内接円による算出

nを1から125,000までの適当な値を設定して近似円周率を計算し、正しい桁数を評価した。

2n+1桁目を四捨五入すると、全ての場合にて、2n桁正しい円周率が得られた。

> restart:with(plots):with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> xdata:=[1,2,3,4,5,6,7,8,9,10,15,20,25,30,35,40,45,50,100,130,200,241,300,400,
500,1000,1301,2000,2401,5000,10000,13001,20000,24999,25000,27500,30000,
35000,50000,125000]:
ydata:=[2,4,6,8,10,12,14,16,18,20,30,40,50,60,70,80,90,100,200,260,400,482,
600,800,1000,2000,2602,4000,4802,10000,20000,26002,40000,49998,50000,
55000,60000,70000,100000,250000]:
A:=augment([seq(1,i=1..40)],xdata):
At:=transpose(A):
p:=evalm(inverse(At &* A) &* At &* ydata):
Correct_Digits:=dotprod(p,vector([1,x]));

[Maple Math]

> lplot:=loglogplot(Correct_Digits,x=0.5..300000,y=0.5..400000,color=brown,
labels=[`n=log[10](N)`,`Correct digits`],thickness=1,
title=`Relation between n=log[10](N) and Correct Digits of Pi`,
titlefont=[TIMES,BOLD,13]):
dplot:=loglogplot([seq([xdata[i],ydata[i]],i=1..40)],
style=point,symbol=circle,color=blue):
plots[display]([lplot,dplot]);

[Maple Plot]

3.3.内接、外接円及び両者の内分による結果

内接多角形にてn=1,50,100の場合は2n+1桁目を四捨五入した。

外接多角形にてn=100の場合は正値の2n桁目を四捨五入すると [Maple Math] となる。

また、内分値にて、n=50,100,1000の場合は正値の 4n桁目を四捨五入すると [Maple Math] となる。

> restart:

> matrix(4,11,[[n , 1, 10, 50, 100, 500, 1000, 5000, 10000, 50000,125000],
[Digit[In], 2, 20, 100, 200, 1000, 2000, 10000, 20000, 100000,250000],
[Digit[Ont], 1, 19, 99, 198, 999, 1999, 9999, 19999, 99999,NotEvaluated],
[Digit[Intpr], 3, 39, 198, 398, 1999, 3998, 19999, 39999, 199999,NotEvaluated]]);

[Maple Math]

3.4.過去のデータの検討

ここで、ArchimedesからLudolfまでのデータを評価する。

> restart:with(plots):

> xdata:=[log[10](96),log[10](96),log[10](1536),log[10](6*2^16),log[10](2^30),
log[10](60*2^29),log[10](2^44),log[10](2^62)]:
ydata:=[3,4,5,10,16,20,26,36]:
pline:=2*x:

> lplot:=plot(pline,x=-2.5..21,y=-5..44,color=aquamarine,
labels=[`n=log[10](N)`,`Correct digits`],
title=`Relation between n=log[10](N) and Correct Digits of Pi`,
thickness=2,titlefont=[TIMES,BOLD,9],axes=FRAME):
dplot:=plot([seq([xdata[i],ydata[i]],i=1..8)],
style=point,symbol=circle,color=red):
archimedes:=textplot([3,0,`Archimedes`],color=blue,font=[TIMES,BOLD,10]):
pisano:=textplot([1,8,`Pisano`],color=blue,font=[TIMES,BOLD,10]):
liu:=textplot([5,4.5,`LiuHui`],color=blue,font=[TIMES,BOLD,10]):
viete:=textplot([7,10,`Viete`],color=blue,font=[TIMES,BOLD,10]):
adriaan:=textplot([10,14,`Adriaan`],color=blue,font=[TIMES,BOLD,10]):
ludolph1:=textplot([12,20,`Ludolph`],color=blue,font=[TIMES,BOLD,10]):
kamata:=textplot([15,25,`Kamata`],color=blue,font=[TIMES,BOLD,10]):
ludolph2:=textplot([17,39,`Ludolph`],color=blue,font=[TIMES,BOLD,10]):
text:=textplot([7,24,`Digits=2n`],color=violet,font=[TIMES,BOLD,15]):
plots[display]([lplot,dplot,archimedes,pisano,liu,viete,
adriaan,ludolph1,kamata,ludolph2,text]);

[Maple Plot]

>

4.参考文献

1. "Encyclopedia of Mathematics" Ed.JpnSoc.Math. (Iwanami Shoten,Tokyo,1972) p.136 (in Japanese) .

2 . E.W.Weisstein, "Concise Encyclopedia of Mathematics" (Chapman & Hal/CRCnetBASE,1999)

3. Y.Kaneda, "Story of [Maple Math] " (Tokyo Tosho,Tokyo, 1997) (in Japanese) .

4. L.BerggrenCJ.Borwein and P.Borwein "Pi; A Source Book" (Springer, Germany, 1997)

5. S.Kobayashi, "Mathematics of a Circle" (Shokabo,Tokyo,1999) (in Japanese).

Part of this article has been read at the Annual Mtg of the Jpn.Soc.Indst&Appl.Math. (Ehime Univ., Matsuyama, Japan, 4th Oct., 1999)

>