Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fortran95 source code to solve approximate equation with least squares method.
!***********************************************************
! 最小二乗法
!***********************************************************
!
!***********************************************************
! 定数モジュール
!***********************************************************
module constants
implicit none ! 暗黙の型指定を使用しない
! データ数
integer, parameter :: N = 7
! 予測曲線の次数
integer, parameter :: M = 5
! 測定データ
integer, parameter :: X(7) = (/ -3, -2, -1, 0, 1, 2, 3 /)
integer, parameter :: Y(7) = (/ 5, -2, -3, -1, 1, 4, 5 /)
end module constants
!***********************************************************
! 計算モジュール
!***********************************************************
module calc
use constants ! 定数モジュール
implicit none ! 暗黙の型指定を使用しない
private ! モジュール外非公開
public :: calc_lsm ! calc_lsm のみモジュール外公開
contains
!*********************************************************
! 最小二乗法
!*********************************************************
subroutine calc_lsm()
! 配列定義
real :: a(M + 2, M + 1)
real :: s(2 * M + 1) = 0, t(M + 1) = 0
! s[], t[] 計算
call calc_st(s, t)
! a[][] に s[], t[] 代入
call ins_st(a, s, t)
! 掃き出し
call sweap_out(a)
! y 値計算&結果出力
call display(a, s, t)
end subroutine calc_lsm
!*********************************************************
! 以下、 private subroutine
!*********************************************************
! s[], t[] 計算
subroutine calc_st(s, t)
real, intent(inout) :: s(:), t(:)
integer :: i, j
do i = 1, N
do j = 1, 2 * M + 1
s(j) = s(j) + x(i) ** (j - 1)
end do
do j = 1, M + 1
t(j) = t(j) + x(i) ** (j - 1) * y(i)
end do
end do
end subroutine calc_st
! a[][] に s[], t[] 代入
subroutine ins_st(a, s, t)
real, intent(inout) :: a(:,:)
real, intent(in) :: s(:), t(:)
integer :: i, j
do i = 1, M + 1
do j = 1, M + 1
a(j, i) = s(i + j - 1)
end do
a(M + 2, i) = t(i)
end do
end subroutine ins_st
! 掃き出し
subroutine sweap_out(a)
real, intent(inout) :: a(:,:)
integer :: i, j, k
real :: p, d
do k = 1, M + 1
p = a(k, k)
do j = k, M + 2
a(j, k) = a(j, k) / p
end do
do i = 1, M + 1
if (i /= k) then
d = a(k, i)
do j = k, M + 2
a(j, i) = a(j, i) - d * a(j, k)
end do
end if
end do
end do
end subroutine sweap_out
! 結果出力
subroutine display(a, s, t)
real, intent(in) :: a(:,:)
real, intent(in) :: s(:)
real, intent(in) :: t(:)
integer :: i, px
real :: p
do i = 1, M +1
write (*, '(a,i1,a,f10.6)'), "a", i - 1, " = ", a(M + 2, i)
end do
write (*, '(a)'), " x y"
do px = -3 * 2, 3 * 2
p = 0
do i = 1, M + 1
p = p + a(M + 2, i) * ((px / 2.0) ** (i - 1))
end do
write (*, '(f5.1,f5.1)'), (px / 2.0), p
end do
end subroutine display
end module calc
!***********************************************************
! 主処理
!***********************************************************
program least_squares_method
use constants ! 定数モジュール
use calc ! 計算モジュール
implicit none ! 暗黙の型指定を使用しない
! ==== 最小二乗法
call calc_lsm()
end program least_squares_method
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment