LoginSignup
0
0

More than 3 years have passed since last update.

@ochimusha01(yasunori matuki)

【数理考古学】とある円周率への挑戦？

Last updated at Posted at 2021-02-05
  • (2021年2月10日)初投稿
  • (2021年4月23日)Tex導入

古来数学者達は、それがx=1の時に$\frac{π}{2}$で収束して円周率に迫るのでアークタンジェント(ArcTangent)近似式の高速化に血道を上げてきたのです。
円周率の公式と計算法

ビエトの無限乗積式(Vieteian Infinite Product Formula)

そうした円周率近似演算の原風景として「(初めて既知数を記号化し記号代数の原理と方法を確立して当時の代数学を体系化した事から)代数学の父」と呼ばれるフランソワ・ビエト(François Viète、1540年~1603年）の無限乗積式(Infinite Product Formula)があり、それ自体は「最初のk項までの積=半径1の円に内接する正2k+1角形の面積πへの収束」あるいは「半円に内接する半多角形の周長への収束)」を求める内容だったのです。

#π=2×2/sqrt(2)×2/sqrt(2+sqrt(2))×2/sqrt(2+sqrt(2+sqrt(2)))…
> 2*2/sqrt(2)
[1] 2.828427
> 2*2/sqrt(2)*2/sqrt(2+sqrt(2))
[1] 3.061467
> 2*2/sqrt(2)*2/sqrt(2+sqrt(2))*2/sqrt(2+sqrt(2+sqrt(2)))
[1] 3.121445
>2*2/sqrt(2)*2/sqrt(2+sqrt(2))*2/sqrt(2+sqrt(2+sqrt(2)))*2/sqrt(2+sqrt(2+sqrt(2+sqrt(2))))
[1] 3.136548
>2*2/sqrt(2)*2/sqrt(2+sqrt(2))*2/sqrt(2+sqrt(2+sqrt(2)))*2/sqrt(2+sqrt(2+sqrt(2+sqrt(2))))*2/sqrt(2+sqrt(2+sqrt(2+sqrt(2+sqrt(2)))))
[1] 3.140331
#ところで√2/2(0.7071068)は√2/(√2^2)。
#なので分子分母共に√2乗すると1/√2。
#すなわち2×2/sqrt(2)(2.828427)=2×sqrt(2)となる。
[1] 2.828427
> sqrt(2)/2
[1] 0.7071068
> 1/sqrt(2)
[1] 0.7071068
> 2*sqrt(2)
[1] 2.828427
> 2*2/sqrt(2)
[1] 2.828427
#当然2/sqrt(2)=sqrt(2)である。
> 2/sqrt(2)
[1] 1.414214
> sqrt(2)
[1] 1.414214
#すなわち
> 2*sqrt(2)
[1] 2.828427
> 2*sqrt(2)*2/sqrt(2+sqrt(2))
[1] 3.061467
#つまりこの計算は半円に内接する半多角形の周長を求めているのに等しい。

半円に内接する半多角形の周長計算と考えた場合の「ビエトの円近似(Vieteian Circle Approximation)」
Rplot23.png

f0<-function(x) sqrt(1-(x-1)^2)
plot(f0,asp=1,xlim=c(0,2),ylim=c(0,1),main="Vieteian Circle Approximation",xlab="Real",ylab="Imaginal")
segments(0,0,2,0,col=rgb(1,0,0),lwd=2)
#ここでおもむろに「ガウスの巡回群による三角関数演算」概念を導入。
par(new=T)#上書き指定
c0=seq(0,pi,length=3)
cx0=cos(c0)+1
cy0=sin(c0)
plot(cx0,cy0,type="l",asp=1,xlim=c(0,2),ylim=c(0,1),main="",xlab="",ylab="",col=rgb(0,0,1),lwd=2)

par(new=T)#上書き指定
c1=seq(0,pi,length=5)
cx1=cos(c1)+1
cy1=sin(c1)
plot(cx1,cy1,type="l",asp=1,xlim=c(0,2),ylim=c(0,1),main="",xlab="",ylab="",col=rgb(0,1,0),lwd=2)

legend("bottomright", legend=c("2","2×2/√2=2.83","2×2/√2×2/√(2+√2))=3.06"),lty=c(1,1,1,1),col=c(rgb(1,0,0),rgb(0,0,1),rgb(0,1,0)),cex=3/4)
  • 「往復(1周)」でなく「片道(半周)」しか見ないので、私の投稿では1の冪乗根を求める演算z^n=1の結果集合たるガウスの巡回群を用いて計算を省力化した「(円を挟む内接多角形と外接多角形の辺数をひたすら増やし続けながら辺長計の差を求める)挟み撃ち定理（Squeeze Theorem）」の半分の計算範囲。
    【初心者向け】挟み撃ち定理による円周率πの近似
    image.gif
  • 一方ユークリッド計量(Euclidean metric)概念を用いた表現では「xy軸上の直線距離{0,0}-{2,0}の中心(0,1)に(ピタゴラスの定理$y=\sqrt{1-x^2}$に従って)等距離の経過点を置いて行く」イメージとなるが、これを一般化した計量関数(Metric Function)は通常xy軸上の直線距離{1,0}{0,1}(0→π/2=90度の範囲)しか見ないのでその倍の計算範囲。
    【初心者向け】方形描画関数②距離空間との関係。
    image.png

こうした視座(Perspective)の概念は、以下の各次元上の制約に対応します。

しばしば「収束が遅く実用的でない」なる一言で片付けられてしまうビエトの円近似法ですが挟み撃ち定理の特殊例に過ぎないとはいえ「(sqrt(2)を含む)素数族2^nのみ導入された環境下でもπ(3.141593)へは到達可能である事を示す意味合いにおいて重要です(sqrt(2)の逆数たる2^2=4すら登場しない)。

#ビエトの円近似には外接多角形の半周の計算に対応する逆算もありそう？
#現段階では未達…

#正方形
>f0(4)
[1] 4
> 2*(2*(2-1))
[1] 4

#正八角形
> f0(8)
[1] 3.313708
> 3.313708/4
[1] 0.828427
> 0.828427/2
[1] 0.4142135
> sqrt(2)-1
[1] 0.4142136
> 2*(2*(2-1))*(2*(sqrt(2)-1))
[1] 3.313708

#正16角形
> f0(16)
[1] 3.182598
> 3.182598/3.313708
[1] 0.9604341
> 0.9604341/2
[1] 0.480217
sqrt(x)-1
> 1.480217^2
[1] 2.191042
[1] 2.191042
> 2*(2*(2-1))*(2*(sqrt(2)-1))*(2*(sqrt(2.191042)-1))
[1] 3.182597

#現段階では、これが精一杯…
  • ちなみにビエトの無限乗積式(Vieteian Infinite Product Formula)に拠る演算結果集合は、0あるいは2を単位元と定め、πに収束する元に-1を掛け、2に収束する演算結果集合を逆元と設定すると結合条件も成立し群定義の少なくとも一部を満たす((興味深い事に単位元を0に定めた時は2,2に定めた時は0が極限値となるので両者が単位元として並び立つ事はない)
    【初心者向け】群論概念①基本定義
  • 問題は演算結果が閉じない事だが(そもそも近似過程で現れる全ての数字を含めた元をイメージする事自体がナンセンス?)、同じ問題なら式形(1+1/N)^Nの時にネイピア数exp(1)(2.718282)、式形(1-1/N)^Nの時にネイピア数の逆数exp(-1)(0.3678794)、式形(1±θi/N)^Nの時に単位元0/極限-1(±π)へと収束する自然指数関数(Natural Exponential Function)e^xも抱えており、こういう場合には実数全体を元に採用するしかない(ちなみに私の投稿では実数を「10進方法の演算結果集合」と規定し、無理数もとりあえず計算可能範囲まで計算した結果をその演算結果集合に「仮に」加えている。コンピューター時代ならではの対「ゴルディアスの結び目(Gordian Knot)」解決法?)。
    【初心者向け】指数・対数関数の発見とそれ以降の発展について。

ところで、ここまでビエトの無限乗積式への関心を掻き立てておいて何なのですが、残念なお知らせがあります。実はこの演算、オイラーの公式(Eulerian Formula,奇しくも著名なe^θi=Cos(θ)+Sin(θ)iと同名)Sin(x)/x(Sinc関数)=cos(x/2)×cos(x/4)×cos(x/6)×cos(x/8)×cos(x/16)…に((0±1i)^2x→e^πi変換の様に)x=π/2を代入した特殊ケースに過ぎなかったのです。
ヴィエトの無限積の公式
sinx/xの極限は？x→0とx→∞の場合を証明付きで東大医学部生が教えます！
image.gif

Macrolin_expansion_exp<-function(x){

#グラフのスケール決定
Graph_scale_x<-c(-40,40)
Graph_scale_y<-c(-1,1)
#関数と凡例の決定
switch(x,
"0"= f0<-function(x) {cos(x/2)},
"1"= f0<-function(x) {cos(x/2)*cos(x/4)},
"2"= f0<-function(x) {cos(x/2)*cos(x/4)*cos(x/8)},
"3"= f0<-function(x) {cos(x/2)*cos(x/4)*cos(x/8)*cos(x/16)},
"4"= f0<-function(x) {cos(x/2)*cos(x/4)*cos(x/8)*cos(x/16)*cos(x/32)},
"5"= f0<-function(x) {cos(x/2)*cos(x/4)*cos(x/8)*cos(x/16)*cos(x/32)*cos(x/64)}
)

switch(x,
"0"=text<-c("cos(x/2)"),
"1"=text<-c("cos(x/2)*cos(x/4)"),
"2"=text<-c("cos(x/2)*cos(x/4)*cos(x/8)"),
"3"=text<-c("cos(x/2)*cos(x/4)*cos(x/8)*cos(x/16)"),
"4"=text<-c("cos(x/2)*cos(x/4)*cos(x/8)*cos(x/16)*cos(x/32)"),
"5"=text<-c("cos(x/2)*cos(x/4)*cos(x/8)*cos(x/16)*cos(x/32)*cos(x/64)"),
)

f1=function(x){sin(x)/x}

plot(f1,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(0,0,0),main="Macrolin expansion", xlab="X", ylab="Y")
par(new=T)#上書き指定
plot(f0,xlim=Graph_scale_x,ylim=Graph_scale_y,type="l",col=rgb(1,0,0),main="", xlab="", ylab="")

legend("bottomleft", legend=c("e^x",text),lty=c(1,1),col=c(rgb(0,0,0),rgb(1,0,0)))
}

library("animation")
Time_Code=c("0","0","0","1","1","1","2","2","2","3","3","3","4","4","4","5","5","5")
saveGIF({
for (i in Time_Code){
　　Macrolin_expansion_exp(i)
}
}, interval = 0.1, movie.name = "sync01.gif") 

それってどういう事? また歴史的価値しかない式と思いきや、ロジスティック方程式/写像との絡みがあって大学受験数学の出題範囲にもなってる模様…
ロジスティック写像と漸化式
【初心者向け】ロジスティック方程式とその関連範囲
ロジスティック写像 - Wikipedia
image.png

オイラー変換(Eulerian Transformation)

オイラーはアークタンジェント(Arctangent)の求め方の高速化も図っています。具体的には1755年に1671年にスコットランドのジェームス・グレゴリーが発見したグレゴリー級数(Gregory Series)のアルゴリズムを大幅に改善したオイラー変換(Eulerian Transformation)を発表。
円周率の公式と計算法
グレゴリー級数
【Rで球面幾何学】オイラーの公式を導出したマクローリン級数の限界？
image.png

ライプニッツ級数(Leibniz series,1674年)とグレゴリー級数(Gregory series, 1671年)とテイラー級数(Taylor Series,1715年)

ライプニッツ級数(Leibniz series,1674年)

\sum_{n}^{\infty}\frac{-1^n}{2n+1}=\frac{π}{4}

グレゴリー級数(Gregory series, 1671年)

\sum_{n}^{\infty}\frac{-1^n}{2n+1}x^{2n+1}=\frac{π}{4}

x=1の時、ライプニッツ級数と一致する。

テイラー展開 - Wikipedia

数学において、テイラー級数(Taylor series)は関数のある一点での導関数たちの値から計算される項の無限和として関数を表したものである。そのような級数を得ることをテイラー展開という。

テイラー級数の概念はスコットランドの数学者ジェームズ・グレゴリーにより定式化され、フォーマルにはイギリスの数学者ブルック・テイラーによって1715年に導入された。0を中心としたテイラー級数は、マクローリン級数(Maclaurin Series)とも呼ばれる。これはスコットランドの数学者コリン・マクローリンにちなんでおり、彼は18世紀にテイラー級数のこの特別な場合を積極的に活用した。

関数はそのテイラー級数の有限個の項を用いて近似することができる。テイラーの定理はそのような近似による誤差の定量的な評価を与える。テイラー級数の最初のいくつかの項として得られる多項式はテイラー多項式と呼ばれる。関数のテイラー級数は、その関数のテイラー多項式で次数を増やした極限が存在すればその極限である。関数はそのテイラー級数がすべての点で収束するときでさえもテイラー級数に等しいとは限らない。開区間（あるいは複素平面の開円板）でテイラー級数に等しい関数はその区間上の解析関数と呼ばれる。

なるほど、以下で私が悩んだ「ベクトルの世界と関数列の世界をどう統合するか」問題の起源はまさにその大源流にあった訳です(ライプニッツの数列をグレゴリーが関数列化し、これにニュートンやライプニッツの手になる微積分概念が合流して「関数列」テイラー級数となった)。
【無限遠点を巡る数理】等差数列と等比数列③基底関数の線形結合?

そんな感じで以下続報…

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0