19

円周率を積分で近似計算する

2338
1

はじめに

どうもこんにちは、🐟️🍊みかん🍊🐟️です。 こんなこと を言ってしまったので、今回は円周率などを高校範囲で有理数評価する方法を書いていきます。この方法は中学生の時に思い付いた方法ですが、いくつかの計算において一般的な計算方法がわかったために簡略化が加わったため、計算自体は非常に簡単なものになっています。そして、僕が知る限り高校範囲内において、一般的に多少の桁数を(高校範囲のみで)求める際には 計算が軽い ことと、他の数学定数の計算にも応用できる点でかなり使い勝手のいい方法だと思います。
本稿を書くにあたって、先行研究を全く調べていないので恐らく既知だと思います。また、π50桁は手計算したことがありますが、 二度とやりたくありません

この記事を読む前に、 匿さんの某積分記事 を読んでおくとよいかもしれません。参考にはなりませんが。

円周率計算

概説

まず、次の等式は当然成り立ちます。

0111+x2dx=π4

こんなの当たり前ですね。そこで、区間(0,1)においてかなり小さい値をもつ関数f(x)を用意して、

01f(x)1+x2dx,01f(x)dx

の二つが高校範囲で計算できてかつ、表式にπと有理数のみが出てこればこちらのものです。具体的にはこんな不等式評価が成り立ちます。

1201f(x)dx<01f(x)1+x2dx<01f(x)dx

証明はほぼ当たり前なので少し天下り的になりますが、結果的に次のような取り方をしました。

f(x)=x4n(1x)4n

はい。ベータ関数です。計算自体は(部分積分を繰り返し用いるなどすれば)可能なのですが、面倒なのでズルをします。

01x4n(1x)4ndx=B(4n+1,4n+1)=(4n)!2(8n+1)!=1(8n+1)(8n4n)

みていただくとわかる通り、分母に中心二項係数が出てくるのでnが大きくなるとかなり急速に0に近づきます。細かいことが知りたければStirling近似とかを使えばいいですね。中の積分も面倒ではありますが、πと有理数しか出てこないことは技術的には簡単に示せるので、実質的問題は

01x4n(1x)4n1+x2dx

を求める作業に帰着されました。

n=1の場合

参考程度にnが小さい場合について具体的に計算します。まず、ベータ積分のほうは

B(5,5)=4!29!=1630

と計算できます。真ん中の積分は普通に割り算して計算してあげると

01x4(1x)41+x2dx=01(x64x5+5x44x2+441+x2)=(1746+143+4)π=227π

であるので、次の不等式を得ます。

11260<227π<1630

これを整理すると、

2271630<π<22711260

となり、この結果だけで余裕で小学校における定説(?)であるπ=3.14を算出できます。この計算はさほど大変ではありません。恐らく、手計算で実用的なレベルの評価としてはn=2くらいが限界な気がします。

n=2の有理数評価もしておくと、地道な計算により
B(9,9)=1218790,01x8(1x)81+x2dx=4π18868415015
なので、

11750320+18868460060<π<1875160+18868460060

であるから、頑張って割り算すると

3.14159231...<π<3.1415928...

によって有効数字7桁を決めることができます。これくらいの計算なら受験本番でもできていいと思います。ちなみに355133よりも少しだけいい近似なくらいなので、覚えたいならむしろこっちのほうがいいですね。

本題

とりあえず精度が欲しいので、n=25くらいにして積分を計算してみます。これは一度だけ手計算したことがあります。当時の僕はこれだと余裕で50桁に届くことがStirling評価からわかるので、50桁だけ計算しました。この結果、なんと(高校一年生の時に教室の後ろの黒板に書いて)写メ取ってたみたいです。見にくいので有理数に関するところだけ赤くします。また、計算過程なんてものはやりたくもないので結果だけ書いていきます。まずベータ積分の方は、

B(101,101)=118200251445876759514246239592574316938775422524758080705105320

となります。それに対して、不等式の内側の積分は

01x100(1x)1001+x2dx=2542258884594828791784019899031095526689322089125551749388128952434249608396205313579242909858566328749487633387736328132736580036532695803192488438380846224993685496448911426147450281474976710656π

にと(8時間くらい頑張って)計算できるので、結果的に不等式は次のような形になります。

136400502891753519028492479185148633877550845049516161410210640<2542258884594828791784019899031095526689322089125551749388128952434249608396205313579242909858566328749487633387736328132736580036532695803192488438380846224993685496448911426147450281474976710656π<118200251445876759514246239592574316938775422524758080705105320

これを両辺281474976710656で割って移項などの整理を行うと、(ここは赤くしてません)

254225888459482879178401989903109552668932208912555174938812895243424960839620531357924290985856638092261362051105765248285301685603166696861085875287087233705766880133617989065376859671442227200110245830703712488149502452713532562564917398605977053684060123252079292579840<π<25422588845948287917840198990310955266893220891255517493881289524342496083962053135792429098585663809226136205110576524828530168560316669686108587528708723370576688013361798906537685967144222720015122915351856244074751226356766281282458699302988526842030061626039646289920

となり、とりあえず凄そうな分数が出てくるみたいです。50桁を超える十進小数だけ何故か過去の研究ノートに結果が書いてないのでWolframに入れようかと思ったのですが...なんと、文字列長の制限を食らいました。式が長すぎるようです。そこで、こんな計算やってられないので こちらの計算機 を使うことにしました。ここからは手計算じゃないですよ。先の分数を十進小数に直した不等式にしてあげると、

3.141592653589793238462643383279502884197169399375105820974944592307816406286267527922991584361811747032194196784852017799295308862<π<3.141592653589793238462643383279502884197169399375105820974944592307816406286169927247271102405380017296915582316302269556264686866

となるらしいので、最右辺と最左辺がどこまで一致しているか色を付けてみると、

LHS=3.141592653589793238462643383279502884197169399375105820974944592307816406286267527922991584361811747032194196784852017799295308862RHS=3.141592653589793238462643383279502884197169399375105820974944592307816406286169927247271102405380017296915582316302269556264686866

ええと、これは何桁なんだ。数えるのも面倒だったので こちらのサイト に数えてもらったところ、76桁の精度であるようだ。精度高くない?

まとめ

円周率の近似値は、おおよそ3.141592653589793238462643383279502884197169399375105820974944592307816406286であることがわかりました!

ln2も計算してみる

2の自然対数も求めることができます。具体的に言うと

1201xn(1x)ndx<01xn(1x)n1+xdx<01xn(1x)ndx

で不等式を作ればよいです。対数の引数の値が大きくなると挟む積分が極限で消えてくれないというバグが発生するのでln5などを求めるのには使えません。(個人的には、210=1024213=8192を繰り返しかけて10のべき乗の近似値を作ってlc2を求めてから、間接的にln5を求めるほうが現実的だと思います)ここではn=1,2,3,4,8のときのln2の有理数近似を計算します。まず積分の結果については

01x(1x)1+xdx=322ln201x2(1x)21+xdx=4ln211401x3(1x)31+xdx=111208ln201x4(1x)41+xdx=16ln26215601x8(1x)81+xdx=256ln242629549240240

になります。ベータ関数のほうは

B(2,2)=16B(3,3)=130B(4,4)=1140B(5,5)=1630B(9,9)=1218790

であることがわかるので、n=1,2,3,4,8の順に

23<ln2<341241116+1240<ln2<1116+112011116011120<ln2<11116012240621896+120160<ln2<621896+1100804262954961501440+1112020480<ln2<4262954961501440+156010240

で、小数表示では大まかには(ノーテーションも雑になっているのは許してください)

0.67<ln2<0.700.6916<ln2<0.69580.6928<ln2<0.69330.69312<ln2<0.693170.6931471775<ln2<0.6931471864

なので、ln2の近似値が0.69314717であることが分かりました。n=4のときによく知られている0.6931まで計算できるのは初めて計算した時は驚きましたね。

おわりに

一番最初にこの記事を書こうと思ったのは この記事 までさかのぼりますが、研究ノートの発掘が面倒だったのでずっとさぼっていました。どうせツイートしてしまったので、作ってみるか、という動機で(対数のほうの)計算をしていたのですが、n=8に関してはさすがにやったことがなく、かなり高い精度で計算できたというのが少し面白かったです。

近似計算が数学において本質的に重要か、と言われれば必ずしもそうでないと思いますが、何処までの精度をもって計算できるかというのはときたま意外な一面を見せるので面白いですね。あと誰か、この方法(あるいはより一般的なもの)を知っている人がいたら教えてください。

最後まで読んでいただきありがとうございました。

投稿日:20231013
OptHub AI Competition

この記事を高評価した人

高評価したユーザはいません

この記事に送られたバッジ

バッジはありません。
バッチを贈って投稿者を応援しよう

バッチを贈ると投稿者に現金やAmazonのギフトカードが還元されます。

投稿者

級数

コメント

他の人のコメント

コメントはありません。
読み込み中...
読み込み中
  1. はじめに
  2. 円周率計算
  3. 概説
  4. n=1の場合
  5. 本題
  6. まとめ
  7. ln2も計算してみる
  8. おわりに