Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to interpolate by Newton's method.
!****************************************************
! ニュートン補間
!
! * 入力はテキストファイルをパイプ処理
! 1行目: 点の数
! 2行目以降: x, y
!
! date name version
! 2018.12.13 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
! SP: 単精度(4), DP: 倍精度(8)
integer, parameter :: SP = kind(1.0)
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const
module newton
use const
implicit none
private
public :: interpolate
contains
! ニュートン補間
!
! :param(in) real(8) pts(2, num_pt): あらかじめ与えられた点の配列
! :param(in) integer(4) num_pt: あらかじめ与えられた点の数
! :param(in) real(8) x: 補間する x 値
! :return real(8) y: 補間後の y 値
function interpolate(pts, num_pt, x) result(y)
implicit none
real(DP), intent(in) :: pts(2, num_pt), x
integer(SP), intent(in) :: num_pt
real(DP) :: y
integer(SP) :: i, j
real(DP) :: c(num_pt), w(num_pt)
! 配列初期化
c(:) = 0.0_DP ! 係数用
w(:) = 0.0_DP ! 作業用
! 差分商
do i = 1, num_pt
w(i) = pts(2, i)
do j = i - 1, 1, -1
w(j) = (w(j + 1) - w(j)) / (pts(1, i) - pts(1, j))
end do
c(i) = w(1)
end do
! 総和
y = c(num_pt)
do i = num_pt - 1, 1, -1
y = y * (x - pts(1, i)) + c(i)
end do
end function interpolate
end module newton
program newton_interpolation
use const
use newton
implicit none
real(DP),allocatable :: pts(:, :) ! あらかじめ与える点
integer(SP) :: num_pt, i
real(DP) :: x, y
! 点の数の取得
read (*, *) num_pt
! 配列メモリ確保
allocate(pts(2, num_pt))
! あらかじめ与えられた点の読み込み
do i = 1, num_pt
read (*, *) pts(:, i)
end do
! ニュートン補間
! (1行1点(x, y)で出力)
do i = int(pts(1, 1)), int(pts(1, num_pt)) * 2
x = real(i, DP) / 2
y = interpolate(pts, num_pt, x)
print '(F7.2, X, F7.2)', x, y
end do
! 配列メモリ解放
deallocate(pts)
end program newton_interpolation
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.