掃き出し法を用いて行列の逆行列を求めるプログラム

掃き出し法を用いて行列の逆行列を求めるプログラム

 

掃き出し法を用いて行列Aの逆行列A⁻¹を求めるプログラムを作ることにする。

行列Aを

　　

とし、その逆行列をB

　　

とすると、

　　

となる。

これより、

　　

という連立方程式を得ることができ、この（２）、（３）、（４）を掃き出し法のプログラムを用いて解かせればよい。

 

掃き出し法のプログラムは前回作ってあるので、それを利用することにし、たとえば、次のようなプログラムを作ればよいだろう。

 

! 掃き出し法を用いて逆行列を求める
real a(10,10), b(10,10), c(10,10)
data a,b,c/100*0, 100*0, 100*0/ ! 行列の初期化

! sample
    a(1,1)=1; a(1,2)=0; a(1,3)=2
    a(2,1)=4; a(2,2)=1; a(2,3)=3
    a(3,1)= 2; a(3,2)=-1; a(3,3)=0

n=3

call gyaku(a,b,n)

write(6,*) 'matrix'
do i=1,n
    write(6,*) (a(i,j), j=1,n)
end do
write(6,*)

write(6,*) 'inverse matrix'
do i=1,n
    write(6,*) (b(i,j), j=1,n)
end do

write(6,*)
write(6,*) 'check'
! 計算のチェック
do i=1, n
    do j=1,n
        sum = 0
        do k=1,n
            sum = sum + a(i,k)*b(k,j)
        end do
        c(i,j)=sum
    end do
end do

do i=1,n
    write(6,*) (c(i,j), j=1,n)
end do

end

! 逆行列を求める計算部
subroutine gyaku(p,q,n)
real p(10,10), q(10,10)
real a(50,51)

do k=1, n
    ! 連立方程式の係数をセット
    do i=1,n
        do j=1,n
            a(i,j)=p(i,j)
        end do
        if (i.ne.k) then
            a(i,n+1)=0
        else
            a(i,n+1) = 1
        end if
    end do
    ! 掃き出し法で解く
    call hakidashi(a,n)
    ! 逆行列をセット   
    do i=1, n
        q(i,k) = a(i,n+1)
    end do
end do
end


! 掃き出し法
subroutine hakidashi(a,n)
real a(50,51)

do k =1, n
    p= a(k,k)
    do j=k+1, n+1
        a(k,j)= a(k,j)/p
    end do
    do i=1, n
        if (i.ne.k) then
            do j=k+1, n+1
                a(i,j)=a(i,j)-a(i,k)*a(k,j)
            end do
        end if
    end do
end do

end

なお、上のプログラムでは、このプログラムで求めたAの逆行列Bが本当にAの逆行列なのか確かめるために、行列Aと行列Bの積ABを求め、それを最後に表示してある。

 

 

real p(10,10), q(10,10)

としているので、１０行１０列以下の逆行列までを求められるようにしてある。この部分を変更すれば、より大きな行列の逆行列を求めることもできる。