LU分解で連立一次方程式を解いてみた。

使用環境など

この記事では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/