Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to interpolate by Lagrange'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 lagrange
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) :: p
! 総和初期化
y = 0.0_DP
! 総和
do i = 1, num_pt
p = pts(2, i)
! 総積
do j = 1, num_pt
if (i /= j) p = p * (x - pts(1, j)) / (pts(1, i) - pts(1, j))
end do
y = y + p
end do
end function interpolate
end module lagrange
program lagrange_interpolation
use const
use lagrange
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 lagrange_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.