Fortranで二分法を使って円周率を求める

Fortranで二分法を実装する

入力条件

与えられたa,bに対して

a<b

であるとする

コード

誤差などは全然気にせず適当に書いている部分があるため厳密な用途には使えない

program nibun

implicit none

real(8):: a,b,c
real(8),parameter :: e =0.00000000001_8
integer(4)::i
real(8) , parameter :: zero = 0.00000000_8

a=0.0
b=5.0

  do i=1,1000
    c=(a+b)/2
    if(abs(pi_func(c)) < e) then
      exit
    else if( (pi_func(a)*pi_func(c)) > zero)then
      a=c
    else if( (pi_func(b)*pi_func(c)) > zero ) then
      b=c
    end if
  end do

write(*,*) c
  stop
contains



  real(8) function pi_func(cc)
    implicit none
    real(8):: cc
    pi_func= dsin(cc/4)-dcos(cc/4)
    return
  end function pi_func
end program nibun

出力

 3.1415926535919425 

(10/12追加) 指摘があったのでコードを修正して精度とかも悪かったのでちょっと手直しした

コンパイラはGNU Fortran (Ubuntu 5.5.0-12ubuntu1~16.04) 5.5.0 20171010

program nibun

implicit none

double precision:: a,b,c
double precision,parameter :: e =0.1E-19
integer(4)::i
double precision , parameter :: zero = 0.0_8

!read(*,*) a,b

a=0.0
b=5.0   
    do i=1,100000
        c=(a+b)/2
        if(abs(pi_func(c)) < e) then
            exit
        else if( (pi_func(a)*pi_func(c)) > zero)then
            a=c 
        else if( (pi_func(b)*pi_func(c)) > zero ) then
            b=c
        end if
    end do

write(*,*) c
write(*,*) '厳密解は 3.1415926535 8979323846 2643383279 5028841971'
contains
double precision function pi_func(cc)
    implicit none
    double precision:: cc
    pi_func= sin(cc/4)-cos(cc/4)
    end function pi_func 
end program nibun                           

出力結果

  3.1415926535897931
厳密解は 3.1415926535 8979323846 2643383279 5028841971

もうちょっと精度がよくなる実装があるのかもですがひとまずここまでとします

参考

https://qiita.com/Ren112358/items/f2c0bbc90810c040a2e5 実装は殆どそこの記事を参考にさせていただいた