Newton法

少し改良。

実行結果

ソース・プログラム

module m_sub
  implicit none
  
  abstract interface
    pure elemental real function t_f(a, x)
      real, intent(in) :: a, x
    end function t_f
  end interface

contains

  recursive function xiter(f, a, xlist) result(res)
    procedure(t_f) :: f
    real, intent(in) :: a, xlist(:)
    real, allocatable :: res(:)
    real :: f0, f1
    f0 = xlist(size(xlist))
    f1 = f(a, f0)
    if ( abs(f1 - f0) < 0.5 * spacing(f0) ) then
      res = xlist
    else  
      res = xiter(f, a, [xlist, f1])
    end if
    return
  end function xiter
  
end module m_sub

program iteration
  use m_sub
  implicit none
  
  call pr(3.0, xiter(f, 3.0, [0.1]), 'inverse')
  call pr(7.0, xiter(f, 7.0, [0.1]), 'inverse')
  call pr(3.0, xiter(g, 3.0, [1.0]), 'square root')
  call pr(7.0, xiter(g, 7.0, [1.0]), 'square root')
  call pr(3.0, xiter(h, 3.0, [1.0]), 'qubic  root')
  call pr(7.0, xiter(h, 7.0, [1.0]), 'qubic  root')
  stop

contains

  subroutine pr(a, y, text)
    real, intent(in) :: a, y(:)
    character(len = *), intent(in) :: text
    print '(*(g0))', "Newton's method : ", text, ' of ', a, ' = ', y(size(y)) 
    print '((7f11.6))', y
    print *
    return
  end subroutine pr

  pure elemental real function f(a, x)  ! 1 / a
    real, intent(in) :: a, x
    f = x * (2 - a * x)
    return
  end function f

  pure elemental real function g(a, x)  ! sqrt(a)
    real, intent(in) :: a, x
    g = 0.5 * (x + a / x)
    return
  end function g

  pure elemental real function h(a, x)  ! a^1/3
    real, intent(in) :: a, x
    h = (2 * x + a / (x * x)) / 3
    return
  end function h

end program iteration