fortran66のブログ

fortran について書きます。

355/113 - π

円周率の近似分数に 22/7 がありますが、この値が円周率{\pi}に如何ほど近いのかは、積分によって評価できることが知られています。

{
\int_0^1{x^4(1-x)^4\over1+x^2}dx=22/7-\pi
}

これは分子を展開して地道に計算することで証明されます。

これを知れば、{355/113-\pi}も同様にして求められるのではないかと思うのが人情でしょう。私も昔、計算機の記号積分に冪を高めた式を入れて、式が出てこないかと試してみたりしましたが、単純には求まらない感じでした。

先日ネットをぶらぶらブラウジングしていると、{355/113-\pi}積分に出会いました。積分式だけは紙切れにメモしておいたのですが、元ページは失念してしまいました。ただググるwikipediaにより詳しい記事がありました。

円周率が22/7より小さいことの証明 - Wikipedia

式は結構ごてごてしていて、以下の通りです。
{\int_0^1{x^8(1-x)^8(25+816x^2)\over3164(1+x^2)}dx}

この積分は、被積分関数の端点での微分値が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