消去法は、直接解法であり、元の数により演算回数が決まる。 ただし除算を含むため、相対誤差が積み重なることがある。 部分 pivot を用いて防ぐことが可能であった。
この一方、
近似解のデータが固定されている場合は定常反復解法、と呼ばれ、
Jacobi 法、Gauss -- Seidel 法、
方程式を数値的に解く場合、数値解がある値へと収束しない場合は、 計算を続けても求まらない場合がある。 ある方程式において、k 回目に得られる数値解を vk とする。 十分大きな数 N 回目までに収束する場合は、 十分小さな数 ε (>0) を用いて、
|vk+1 - vk| < ε
が成立することを意味する。
n 元連立 1 次方程式から作られる n × n 行列 A の各要素について、
|ai,i| > &Sigmanj = 1, j ≠ i |ai,j|
(1 ≤ i ≤ n)
が成り立つとき、A は
0 の要素が多い
行列 A x = b として表したものを、x について解き、成分で表す。
&Sigmaj = 0n-1 ai, j xj = bi &Sigmaj = 0i - 1 ai, j xj + ai,i xi + &Sigmaj = i + 1n-1 ai, j xj = bi
数値解を求める際には、仮置きの数値を代入する。k - 1 回目の演算の際、 得られた解の組を {xi(k-1)} とする。 ai,jが対角成分と、対角成分を持たない下三角行列、 および対角成分を持たない上三角行列に分けられていることに注意しよう。 これらをそれぞれ D, L, U と書く。A = D + L + U となるから、
(D + L + U)x = b (D + L) x= b - U x
これを漸化式と見なし、k 回目の演算は {xi(k-1)} で書けることを用いると、
(D + L) x(k)= b - U x(k-1) D x(k) = b - U x(k-1) - L x(k)
成分で書いて、xi(k) について解くと、
xi(k) = ai,i-1(bi - &Sigmaj = 0i - 1 ai, j xj(k-1) - &Sigmaj = i + 1n-1 ai, j xj(k))
Ruby によるデータの読み込みは、open -- end を使用し、
open('ファイル名','r') do |ファイル変数|
while 内部変数 = ファイル変数.gets
if /(1列目の正規表現)\s+(2列目の正規表現).../ =˜ 内部変数
格納する変数 = $1
格納する変数 = $2
:
end
end
end
であった。形式は学部では、TAB 区切りによるものを、正規表現を使用して分解した。 数値列は \d+ や \d+\.\d+ を用い、to_i または to_f 変換を必要とした。
Fortran では、open -- read -- close を使用する。データを読み込む例を見てみよう。
program read
implicit none
integer :: fr = 15, fw = 26
integer a(0:2,0:3), i, j
real(8) r(0:2,0:3)
open(fr, file = 'test.dat')
read(fr, *) a(0,0:3)
read(fr, *) a(1,0:3)
read(fr, *) a(2,0:3)
close(fr)
do i = 0, 2
write(*,*) a(i,:)
end do
do i = 0, 2
do j = 0, 3
r(i,j) = dble(a(i,j))
end do
end do
open(fw, file = 'result.dat')
do i = 0, 2
write(fw,'(100e12.4)') r(i,:)
end do
close(fw)
end program read
データは整数値を読み込み、dble 関数で整数値を倍精度変換している。 read / write の () の左側は、ファイルのラベルで、 5, 6 以外の番号を自由につけてよい。5 は標準入力 (キーボード入力)、 6 は標準出力 (画面出力) である。read(5,*) は read(*,*), write(6,*) は write(*,*) である。 読み込み専用のファイルには、action='read' などのオプションや、 ファイルの読み込み確認のための iostat なども付け加えることができる。
LDU 分解するプログラムをつくろう。
program ldu
implicit none
integer :: fr = 15
integer a(0:2,0:2), b(0:2), i, j
real(8) r(0:2,0:2), br(0:2)
real(8) :: d(0:2,0:2) = 0.0d0 , l(0:2,0:2) = 0.0d0 ,
& u(0:2,0:2) = 0.0d0
open(fr, file = 'test.dat')
do i = 0, 2
read(fr, *) a(i,0:2), b(i)
end do
close(fr)
do i = 0, 2
br(i) = dble(b(i))
do j = 0, 2
r(i,j) = dble(a(i,j))
end do
end do
do i = 0, 2
do j = 0, 2
if (j == i) then
d(i,j) = r(i,j)
else if (i > j) then
l(i,j) = r(i,j)
else
u(i,j) = r(i,j)
end if
end do
end do
do i = 0, 2
write(*,*) d(i,:)
end do
do i = 0, 2
write(*,*) l(i,:)
end do
do i = 0, 2
write(*,*) u(i,:)
end do
end program ldu
長い宣言文になったため、6 列目以前に & をつけ、 前の行からの続きであることを示している。
b≠ 0 かつ aii≠ 0、ならば、 T = (D - L)-1 U から作られる行列が収束行列、すなわち k が無限大で T は限りなく零行列に近づく。 このとき、行列 T の 固有値の最大値を max|λi| とすると、
ρ(T) = max|λi| > 1 (i > i > n)
が成り立つ。A と A からなる対角行列が正定値行列の場合は、収束解を得ることができる。
module matmod
implicit none
contains
function sddm(ra) result(sd)
real(8), intent(in) :: ra(0:,0:)
real(8) sd(0:size(ra,1)-1,0:size(ra,2)-1), v(0:size(ra,1)-1)
integer i, j, k, l
k = size(ra(1,:))
v(0:k-1) = 0.0d0
sd(:,:) = ra(:,:)
do i = 0, k - 1
do j = 0, i - 1
v(i) = v(i) + abs(sd(i, j))
end do
do j = i + 1, k
v(i) = v(i) + sd(i,j)
end do
if (abs(ra(i,i)) > v(i)) then
sd(0:k - 1,0:k - 1) = ra(0:k - 1,0:k - 1)
else
sd(0:k - 1,0:k - 1) = 0.0d0
end if
end do
end function sddm
function trace(ra) result(tr)
real(8), intent(in) :: ra(0:,0:)
real(8) tr
integer i, k
k = size(ra,1)
tr = 0.0d0
do i = 0, k
if (ra(i,i) == 0.0d0) then
stop 'diag i is 0'
end if
tr = tr + ra(i,i)
end do
end function trace
function prev(ra) result(pr)
real(8), intent(in) :: ra(0: ,0:)
real(8) pr(0 : size(ra, 1) - 1, 0:size(ra, 2) - 1)
real(8) pv(0 : size(ra, 2) - 1), pv2(0 : size(ra, 2) - 1)
integer i, j, k, l
l = size(ra, 2)
k = size(ra, 1)
pr = ra
do i = 0, k - 2
j = i + 1
if (ra(i, i) == 0.0d0) then
pv(:) = ra(i, :)
pv2(:) = ra(j, :)
pr(i, :) = pv2(:)
pr(i + 1, :) = pv(:)
end if
end do
end function prev
subroutine gauss_seidel(ra, rb, rx, maxit, eps)
real(8), intent(in) :: ra(0: ,0:), rb(0:), eps
integer, intent(in) :: maxit
real(8), intent(out) :: rx(0: size(ra, 1) - 1)
real(8) ry(0:size(ra,1) - 1)
integer i, j, k, ct
k = size(ra,2)
rx(:) = 0.0d0
do ct = 1, maxit
do i = 0, k - 1
rx(i) = -1.0d0 * dot_product(ra(i, 0: i - 1), rx(0 : i - 1))
rx(i) = rx(i)
& -1.0d0 *
& dot_product(ra(i, i + 1 : k -1), rx(i + 1 : k -1))
rx(i) = rx(i) + rb(i)
rx(i) = rx(i) / ra(i, i)
ry(:) = rb(:) - matmul(ra,rx)
end do
if (dot_product(ry,ry) <= eps) then
write(*,*) ct
exit
end if
end do
end subroutine gauss_seidel
end module matmod
program gs
use matmod
implicit none
integer :: n = 0, fr = 15, maxit = 10000
real(8) :: eps = 1.0d-6
integer l
integer io
real(8), allocatable :: a(:,:), b(:), x(:), sum(:)
integer i, j
open(fr, file = 'test.dat', action = 'read')
do
read(fr, *, iostat = io) l
if (io < 0) then
exit
else
n = n + 1
end if
end do
close(fr)
allocate (a(0:n-1,0:n-1))
allocate (b(0:n-1))
allocate (x(0:n-1))
allocate (sum(0:n-1))
open(fr, file = 'test.dat', action = 'read')
do i = 0, n-1
read(fr, *, iostat = io) a(i, 0 : n - 1), b(i)
a(i,0:n-1) = dble(a(i, 0: n - 1))
b(i) = dble(b(i))
end do
close(fr)
a = sddm(a)
do i = 0, n-1
write(*,'(100e12.4)') a(i, 0 : n - 1), b(i)
end do
a = prev(a)
write(*,*) trace(a)
do i = 0, n-1
write(*,'(100e12.4)') a(i, 0 : n - 1), b(i)
end do
call gauss_seidel(a, b, x, maxit, eps)
do i = 0, n-1
write(*,'(100e12.4)') x(i)
end do
end program gs
Gauss--Seidel 法のプログラムのおよその形はこのようになる。