Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran 95 source code to compute a definite integral.
!****************************************************
! 定積分(台形則、シンプソン則)
! * f = sqrt(4.0 - x**2) (0.0 <= x <= 2.0)
!
! date name version
! 2018.10.13 mk-mode.com 1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module integrator
implicit none
! SP: 単精度(4), DP: 倍精度(8)
integer, parameter :: SP = kind(1.0)
integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
real(DP), parameter :: PI = atan(1.0_8) * 4.0_DP ! 円周率
contains
! 数値積分(台形則)
!
! :param real(8) f(:)
! :param real(8) dx
! :return real(8) ret
function trapezoid(f, dx) result(ret)
implicit none
real(DP), intent(in) :: f(:)
real(DP), intent(in) :: dx
real(DP) :: ret
integer(SP) :: i, n
n = size(f)
ret = 0.5_DP * (f(1) + f(n))
do i = 2, n - 1
ret = ret + f(i)
end do
ret = ret * dx
end function trapezoid
! 数値積分(シンプソン則)
!
! :param real(8) f(:)
! :param real(8) dx
! :return real(8) ret
function simpson(f, dx) result(ret)
implicit none
real(DP), intent(in) :: f(:)
real(DP), intent(in) :: dx
real(DP) :: ret
integer(SP) :: i, n
n = size(f)
! 端点を含めた配列サイズが奇数でなければ終了
if (mod(n, 2) == 0) then
write(*,*) '配列サイズが奇数でない!'
stop
end if
ret = f(1) + f(n)
! (偶数)
do i = 2, n - 1, 2
ret = ret + 4.0_8 * f(i)
end do
! (奇数)
do i = 3, n - 2, 2
ret = ret + 2.0_DP * f(i)
end do
ret = ret * dx / 3.0_DP
end function simpson
end module integrator
program definite_integral
use integrator
implicit none
integer(SP) :: n, nmax
real(DP) :: x_1, x_2, dx, integral
real(DP), allocatable :: f_1(:), f_2(:)
write(*, fmt='(A)', advance='no') '区間の数 : '
read(*,*) nmax
! 配列用メモリ確保
allocate(f_1(0:nmax))
allocate(f_2(0:nmax))
x_1 = 0.0_DP
x_2 = 2.0_DP
dx = (x_2 - x_1) / nmax
! 数値積分(台形則)
do n = 0, nmax
f_1(n) = f(x_1 + dx * real(n, 8))
end do
integral = trapezoid(f_1, dx)
write (*, fmt='(a15, ": ", f18.15, " (誤差= ", e18.12, ")")') &
& "台形則", integral, abs(integral / PI - 1.0_DP)
! 数値積分(シンプソン則)
do n = 0, nmax
f_2(n) = f(x_1 + dx * real(n, 8))
end do
integral = simpson(f_2, dx)
write(*, fmt='(a18, ": ", f18.15, " (誤差= ", e18.12, ")")') &
& "シンプソン則", integral, abs(integral / PI - 1.0_DP)
! 配列用メモリ解放
deallocate(f_1)
deallocate(f_2)
stop
contains
! 被積分関数
real(DP) function f(x)
implicit none
real(DP), intent(in) :: x
f = sqrt(4.0_DP - x * x)
end function f
end program definite_integral
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.