Fortran 95 source code to interpolate by Newton's method.

newton_interpolation.f95

!****************************************************
! ニュートン補間
!
! * 入力はテキストファイルをパイプ処理
!     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