Skip to content

Instantly share code, notes, and snippets.

@komasaru
Last active December 13, 2018 01:40
Show Gist options
  • Save komasaru/30af36adf434d782d5dca4e099096d27 to your computer and use it in GitHub Desktop.
Save komasaru/30af36adf434d782d5dca4e099096d27 to your computer and use it in GitHub Desktop.
Fortran 95 source code to solve nonlinear-equation with newton method.
!****************************************************
! 非線形方程式の解法(ニュートン法)
! * 方程式: y = x**3 - x + 1
!
! date name version
! 2018.10.12 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
program nonlinear_equation_newton
implicit none
! SP: 単精度(4), DP: 倍精度(8)
integer, parameter :: SP = kind(1.0)
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
integer, parameter :: NMAX = 20
real(DP), parameter :: EPS = 1.0e-6
logical :: stat
integer :: n
real(DP) :: x, y, dx, dy
write (*, "(a)", advance="no") "x : "
read (*,*) x
stat = .false.
do n = 1, NMAX
! 次の値の推定
y = f(x)
dy = df(x)
dx = -y / dy
x = x + dx
! 収束判定
if (abs(dx) < EPS) then
stat = .true.
write (*, '("収束 [", i4, "]")') n
exit
else
write(*,fmt='("誤差 [", i4, "] = ", e20.8)') n, abs(y)
end if
end do
! 結果出力
if (.not. stat) then
write (*, *) "近似不可!"
end if
write (*, '("近似値 = ", e20.8)') x
write (*, '("誤差 = ", e20.8)') abs(y)
stop
contains
! 方程式
! * f = x**3 - x + 1
!
! :param real(8) x
! :return real(8) f
real(DP) function f(x)
implicit none
real(DP), intent(in) :: x
f = x**3 - x + 1
end function f
! 方程式の導関数
! * f' = 3* x**2 - 1
!
! :param real(8) x
! :return real(8) df
real(DP) function df(x)
implicit none
real(DP), intent(in) :: x
df = 3 * x**2 - 1
end function df
end program nonlinear_equation_newton
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment