使用環境など
この記事ではWindows10にインストールしたmingw32-gfortranを使っています。
言語は、fortran95を使っています。
LU分解とは
LU分というのは、行列Aを2つの行列LとUに分解することを言います。
これを用いることで、連立一次方程式の計算や行列式計算ができるようになります。
LU分解を使って連立一次方程式Ax=bを解くことのメリットとしては、一度AをLU分解をしたらそのAについての解をbの値を変えるだけで、解を求めることができます。
LUを用いて連立一次方程式を解く方法
先ず、
$$ A\vec{x} = \vec{b} $$
をLU分解します。つまり、以下の式が求まります。
$$ LU\vec{x} = \vec{b} $$
この式を以下のように考えます。
$$ L\vec{y} = \vec{b} $$
$$ U\vec{x} = \vec{y} $$
下の式を上の式に代入することで、解を簡単に求めることができます。
ソースコード
LU_decomposition.f95
program main
implicit none
real(8), allocatable, dimension(:,:) :: a,L,U
real(8), allocatable, dimension(:) :: b,x,y
real(8) temp, alpha, reserve, sum1, sum2
integer i, j, k, m, n, dim, maxline
print *, 'input the number of dimension'
read(*,*) dim
allocate (a(dim, dim), b(dim), L(dim,dim), U(dim,dim),x(dim),y(dim))
! Initialize ------------------------------------------------------
a = 0.0
b = 0.0
L = 0.0
U = 0.0
! Input -----------------------------------------------------------
! A
do i = 1, dim
do j = 1, dim
print *,"x(",i,',',j,')'
read *,a(i,j)
end do
end do
! b
do i = 1, dim
print *,"b(",i,')'
read *,b(i)
end do
! print the input -----------------------------------------------
print *, 'A'
do i = 1, dim
print *, a(i, 1:dim)
end do
print *, 'B'
print *, b(1:dim)
! Solve the problem------------------------------------------------
! LU --------------------------------------------------------------
do k = 1, dim
L(k,k) = 1.0
end do
do i = 1, dim
do j = 1, dim
if ( i .le. j ) then
sum1 = 0.0
do m = 1, i -1
sum1 = sum1 + L(i,m) * U(m,j)
end do
U(i,j) = a(i,j) - sum1
end if
if (i .gt. j) then
sum2 = 0.0
do n = 1, j-1
sum2 = sum2 + L(i,n) * U(n,j)
end do
L(i,j) = (a(i,j) - sum2) / U(j,j)
end if
end do
end do
! ----------------------------------------------------------------
! print for check-------------------------------------------------
print *, 'L'
do i = 1,dim
print *, L(i,1:dim)
end do
print *, 'U'
do i = 1, dim
print *, U(i,1:dim)
end do
! solve the problem-----------------------------------------------
! Ly=b------------------------------------------------------------
do i = 1, dim
y(i) = b(i)
end do
do j = 1, dim - 1
do i = j + 1, dim
y(i) = y(i) - y(j) * L(i,j)
end do
end do
! Ux=y------------------------------------------------------------
do i = 1, dim
x(i) = y(i)
end do
do j = dim, 1, -1
x(j) = x(j) / U(j,j)
do i = 1, j - 1
x(i) = x(i) - U(i,j) * x(j)
end do
end do
! ----------------------------------------------------------------
! print-----------------------------------------------------------
print *, 'x'
do i = 1, dim
print *, x(i)
end do
! ----------------------------------------------------------------
deallocate (a,b,L,U,x,y)
end
実行結果
LUdecomposition.exe
input the number of dimension
2
x( 1 , 1 )
1
x( 1 , 2 )
2
x( 2 , 1 )
3
x( 2 , 2 )
1
b( 1 )
2
b( 2 )
1
A
1.0000000000000000 2.0000000000000000
3.0000000000000000 1.0000000000000000
B
2.0000000000000000 1.0000000000000000
L
1.0000000000000000 0.0000000000000000
3.0000000000000000 1.0000000000000000
U
1.0000000000000000 2.0000000000000000
0.0000000000000000 -5.0000000000000000
x
0.0000000000000000
1.0000000000000000
参考文献：
「物理のかぎしっぽ」http://hooktail.org/computer/index.php?LU%CA%AC%B2%F2
「numerical algorithms group」http://www.nag-j.co.jp/fortran/