数学の考え方参照。
解析的に解く場合は自由に計算することができるが、 数値計算では総当たり、かつ誤差の少ないような計算方法を探しながら、 求めていくことになる。 掃き出し法を使った数値計算はいくつか知られている。 Jordan-Gauss の方法や、Gauss の消去法がよく知られている。
行基本変形を Ruby 言語で再現してみよう。
正方行列の場合、0 ベクトルの行を含む場合など、 解が縮退することがある。 縮退しない空間の次元をランクという。
何度も使用するプログラムを外へ出し、分かりやすくする書き方をする。 Ruby 言語では method (def -- end) であった。 もっとも最後に書かれている値を返す。
def 関数(仮引数) 定義 返す値 end 関数(変数)
Fortran では、例えば module subroutine を使用して繰り返しを外に出す。
module モジュール implicit none contains subroutine サブルーチン(仮引数a,仮引数b) 型, intent(in) :: 仮引数a 型, intent(out) :: 仮引数b ... : 実行文 : end subroutine サブルーチン end module モジュール program 主プログラム use モジュール implicit none 型 :: 変数x = 初期値設定 型 変数y ! 関数からの結果を得るために用意する変数 : : call サブルーチン(変数x,変数y) : end program 主プログラム 関数(変数)
となる。 取り込む仮引数に intent(in), 戻す仮引数に intent(out) を明示する。
第 2 講で扱った連立方程式を行列計算により求めることができる。 x=[x,y] 、b=[e,f] という縦ベクトル、 A = [[a,b],[c,d]] という行列で書くとき、 連立方程式は
A x = b a b x e ( ) ( ) = ( ) c d y f
と表せる。これを、Gauss-Jordan の方法によって、解けばよい。 すなわち、(A|b)=(I|b˜) なるb˜ が得られる解となる。
先頭の成分を 1 にするような行基本変形を行う作業を、pivot という。 Gauss-Jordan では、上から順に pivot を作り、 その他の行では pivot した行に、 その行で初めてゼロでない行の成分をかけたものを引く行基本変形を行う。 すなわち pivot 行ととった成分の存在する列がゼロとなるように行基本変形を行う。 pivot 行の列のみ 1 で、その他が 0 となったら、 次の行へ進み、その行の pivot を作る。 最後の行まで単位行列を作ると、(I|b˜) を得る。
Ruby 言語による Gauss-Jordan の解法例は、例えば次のようになる。
#!/usr/koeki/bin/ruby
require 'mathn'
def tr_mat(a)
tr = Array.new
for i in 0 .. a[0].length - 1
col = Array.new
for j in 0 .. a.length - 1
col << a[j][i]
end
tr << col
end
tr
end
def id_mat(a)
id = Array.new
for i in 0 .. a.length - 1
col = Array.new
for j in 0 .. a[0].length - 1
if j == i
col << 1
else
col << 0
end
end
id << col
end
id
end
def mul(p,q)
multi = Array.new
mul_p = tr_mat(p)
for i in 0 .. mul_p.length - 1
col = Array.new
for j in 0 .. q[0].length - 1
sum = 0
for k in 0 .. mul_p[0].length - 1
sum += mul_p[i][k] * q[k][j]
end
sum
col << sum
end
if col != Array.new(mul_p.length,0)
multi << col
end
end
multi
end
def p_mat(a,i,c)
p_m = id_mat(a)
p_m[i][i] = c
mul(p_m, a)
end
def q_mat(a,i,j,c)
q_m = id_mat(a)
q_m[i][j] = c
mul(q_m, a)
end
def r_mat(a,i,j)
r_m = id_mat(a)
k = r_m[i]
l = r_m[j]
r_m[j] = k
r_m[i] = l
r_m
mul(r_m, a)
end
def pivot(a,i)
c = Array.new
for j in 0 .. a[0].length - 1
c << a[i][j] / a[i][i]
end
a[i] = c
a
end
mat = [[1,2,3,1],[4,5,6,2],[7,8,0,3]]
for i in 0 .. mat.length - 1
mat = pivot(mat,i)
for j in 0 .. mat.length - 1
if j != i
c = -1 * mat[j][i]
mat = q_mat(mat,i,j,c)
end
end
end
ans = Array.new
for i in 0 .. mat.length - 1
ans << mat[i].pop
end
p ans
tr_mat は転置行列を作る。 id_mat は単位行列を作る。 mul は行列の積を計算する。 上記の行基本変形に対応する行列をそれぞれ、p_mat, q_mat, r_mat とした。
[A,I] を新たな行列(A|I)とし、 行基本変形の規則にそって、(I|A-1) に変形する。 Gauss-Jordan の方法で、先頭の n × n 行列を単位行列にすると、 逆行列が後半の行列に出来上がる。 答えを出す部分で、n × n+1 行列ではない場合に、行列そのものを返すことにする。
if mat[0].length - mat.length == 1
ans = Array.new
for i in 0 .. mat.length - 1
ans << mat[i].pop
end
p ans
else
p mat
end
今までの例では、0 、あるいは 0 に近い小さな数で割らないような例を扱う場合だけであった(正則行列)。 pivot 行の先頭のある列全てが 0 の場合に計算を止め、また、 pivot にする行 a(k,k) をそれより下の列からなる部分列 (l,k)(l = k, .. , i) で最大の値を取るものを使用することに変更したい。 この際、本来 pivot にした行と、探し出した行を交換することにする。
選択的に大きな数で割るので、除算による誤差が少なくなる。
pp = Array.new
for l in i .. mat.length - 1
pp << mat[l][i]
end
if pp.rindex(pp.sort.pop) != 0
mat = r_mat(mat, i, pp.rindex(pp.sort.pop) + i)
end
module 内の subroutine や function に配列を渡すときは、
Fortran は C 言語と同じように、配列名とその大きさをもって、
配列の受け渡しを行うことも可能である。これを
module subpart
implicit none
contains
subroutine dummy_print(e,k)
integer, intent(in) :: k
real(8), intent(in) :: e(0:k-1,0:k-1)
integer i, j
do i = 0, k - 1
write(*,'(100e12.4)') (e(i,j), j = 0, k-1)
end do
end subroutine dummy_print
end module subpart
program explicit
use subpart
implicit none
real(8), allocatable :: e(:,:)
integer :: n = 5
integer i, j
allocate (e(0:n-1,0:n-1))
do i = 0, n - 1
do j = 0, n - 1
if (j == i) then
e(i,j) = 1.0d0
else
e(i,j) = 0.0d0
end if
end do
end do
call dummy_print(e,size(e(1,:)))
end program explicit
大きさがいくつかどうかを調べたくない場合は、size(配列(次元,:)) も使える。
Ruby 言語で method を定義したときのように、配列をそのまま受け渡すこともできる。これを
外部 subroutine で、形状引き継ぎ配列を使用する場合には、interface module を用意する必要がある。 当面は、内部 subroutine のみを考える。
module subpart
implicit none
contains
subroutine dummy_print2(e, lower)
integer, intent(in) :: lower
real(8), intent(in) :: e(lower:,lower:)
integer i, j, k
k = size(e(1,:))
do i = 0, k
write(*,'(100e12.4)') (e(i,j), j = 0, k)
end do
end subroutine dummy_print2
end module subpart
program assume
use subpart
implicit none
real(8), allocatable :: e(:,:)
integer :: n = 5
integer i, j
allocate (e(0:n-1,0:n-1))
do i = 0, n - 1
do j = 0, n - 1
if (j == i) then
e(i,j) = 1.0d0
else
e(i,j) = 0.0d0
end if
end do
end do
call dummy_print2(e,lbound(e,1))
end program assume
module subroutine の場合いずれも intent(out) を明示し、 配列に値を渡して戻してもよいが、 渡す際にも答えを得る配列を用意しておかなくてはならない。 inner subroutine で答えを得たい場合には、module function の方が便利である。
以前の FORTRAN 77 ではスカラ関数のみであったが、 配列を返すことができるようになった。 また、部分配列の組み込み演算も可能である。
function 関数(仮引数) result(返す変数) 型, intent(in) :: 仮引数 型 内部変数 返す変数への代入文 end function 関数
module op
implicit none
contains
function unit(mat) result(e)
real(8), intent(in) :: mat(0:,0:)
real(8), allocatable :: e(:,:)
integer j, l, k
k = size(mat,1)
l = size(mat,2)
allocate(e(0:k-1,0:l-1))
e(0:k-1,0:l-1) = 0.0d0
do j = 0, k-1
e(j,j) = 1.0d0
end do
end function unit
function pci(c,i,mat) result(p)
real(8), intent(in) :: c, mat(0:,0:)
real(8), allocatable :: p(:,:)
integer, intent(in) :: i
integer s, t, k, l
k = size(mat,1)
l = size(mat,2)
allocate(p(0:k-1,0:l-1))
p = unit(p)
p(i,i) = c * p(i,i)
p = matmul(transpose(p),mat)
end function pci
function qcij(c,i,j,mat) result(q)
real(8), intent(in) :: c, mat(0:,0:)
real(8), allocatable :: q(:,:)
integer, intent(in) :: i,j
integer s, t, k, l
k = size(mat,1)
l = size(mat,2)
allocate(q(0:k-1,0:l-1))
q = unit(q)
q(i,j) = c * 1.0d0
q = matmul(transpose(q),mat)
end function qcij
function rij(i,j,mat) result(r)
real(8), intent(in) :: mat(0:,0:)
real(8), allocatable :: r(:,:), vec(:)
integer, intent(in) :: i,j
integer s, t, k, l
k = size(mat,1)
l = size(mat,2)
allocate(r(0:k-1,0:l-1))
allocate(vec(0:l-1))
r = unit(r)
vec(0:l-1) = r(i,:)
r(i,0:l-1) = r(j,0:l-1)
r(j,0:l-1) = vec(0:l-1)
r = matmul(transpose(r),mat)
end function rij
end module op
program gauss
use op
implicit none
integer :: row = 4, col = 5
integer m
integer i, j, p, q
real(8) c
real(8), allocatable :: mat(:,:), cal(:,:)
real(8), allocatable :: b(:), x(:)
allocate(mat(0:row-1,0:col-1))
allocate(cal(0:row-1,0:col-1))
allocate(b(0:row-1))
allocate(x(0:row-1))
mat(0,0:col-1) = (/ 2.0d0, 3.0d0, 1.0d0, -3.0d0, 1.0d0/)
mat(1,0:col-1) = (/ -1.0d0, 2.0d0, 2.0d0, 4.0d0, 6.0d0/)
mat(2,0:col-1) = (/ 4.0d0, 1.0d0, -3.0d0, 5.0d0, 3.0d0/)
mat(3,0:col-1) = (/ 5.0d0, -4.0d0, -4.0d0, 1.0d0, 3.0d0/)
do m = 0, row - 1
write(*,'(100e12.4)') mat(m, 0:col-1)
end do
write(*,*)
cal = mat
do i = 0, row - 1
if (cal(i,i) == 1.0d0) stop 'coeff is nan'
do j = i + 1, row - 1
if (abs(cal(i,j)) > cal(i,i)) then
cal = rij(i,j,cal)
end if
end do
c = 1.0d0 / cal(i,i)
cal = pci(c, i, cal)
do j = i+1, row -1
cal = qcij(-1.0d0 * cal(j,i),i,j,cal)
end do
end do
write(*,*)
do m = 0, row - 1
write(*,'(100e12.4)') cal(m, 0:col-1)
end do
write(*,*)
if (row + 1 == col) then
do m = 0, row - 1
p = col - 2 - m
b(p) = cal(p,col - 1)
x(p) = b(p)
do q = p + 1, col - 2
x(p) = x(p) - 1.0d0 * cal(p,q) * x(q)
end do
end do
write(*,'(100e12.4)') x
end if
end program gauss
pivot の代入が本文内で短く済むため、わざわざ呼び出さなかったりしているが、 組み込み関数がいくつか定義されていたり、 Ruby で先に作った Gauss の消去法によるプログラムを焼き直したものである。