円周率の近似分数に 22/7 がありますが、この値が円周率に如何ほど近いのかは、積分によって評価できることが知られています。
これは分子を展開して地道に計算することで証明されます。
これを知れば、も同様にして求められるのではないかと思うのが人情でしょう。私も昔、計算機の記号積分に冪を高めた式を入れて、式が出てこないかと試してみたりしましたが、単純には求まらない感じでした。
先日ネットをぶらぶらブラウジングしていると、の積分に出会いました。積分式だけは紙切れにメモしておいたのですが、元ページは失念してしまいました。ただググるとwikipediaにより詳しい記事がありました。
円周率が22/7より小さいことの証明 - Wikipedia
式は結構ごてごてしていて、以下の通りです。
この積分は、被積分関数の端点での微分値が0なので台形公式でいい値が出そうな気がします。そこで台形積分で計算してみることにしました。
プログラム
面倒なので文関数などを使ってみました。解析解による近似の度合いの見積もりは引き算による桁落ちで精度が下がるので、倍々精度で計算しています。
program pi_approx implicit none integer, parameter :: kd = kind(0.0q0) real(kd), parameter :: pi = 4 * atan(1.0_kd) integer, parameter :: n = 64 + 1 real(kd) :: x(n), y(n), h, s, f integer :: i ! statement function f(x) = x**8 * (1 - x)**8 * (25 + 816 * x**2) / ( 3164 * (1 + x**2 ) ) ! h = (1.0_kd - 0.0_kd) / (n - 1) forall (i = 1:n) x(i) = 0.0_kd + (i - 1) * h y(i) = f(x(i)) end forall ! Trapezoidal rule s = sum(y) - (y(1) + y(n)) / 2 s = s * h print *, s, 355.0_kd / 113.0_kd - pi end program pi_approx
実行結果
分割が少ない割に解析値と良い一致が出ています。
2.667641890624304779248581088571583E-0007 2.667641890624223123689328865557821E-0007