計算物理学２第６回：連立方程式の解法
ver. 2018/5/21
n 元 1 次連立方程式
a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn =b1 ,a21 x1 + a22 x2 + a23 x3 + · · · + a2n xn =b2 ,a31 x1 + a32 x2 + a33 x3 + · · · + a3n xn =b3 ,..
.
an1 x1 + an2 x2 + an3 x3 + · · · + ann xn =bn ,の解法について説明します。大きく分けて、消去法 (elimination) と反復法 (iteration) による解法があります。消去法は係数 a や b の値を変更して解を導出する方法、反復法は近似解を更新してゆくことで真の解に近づけてゆく方法です。

1 消去法
1.1 Gauss の消去法
先程の連立一次方程式は行列の形で

Ax = b(1)と書けます。ここで A は n × n の行列

 (1)
a11
 (1)
a21
 (1)
(1)
A =
a31
 ..
 .
(1)

an1(1)(1)a12(1)
a22(1)a32..
.

a13(1)
a23(1)a33..
.

·········..
.

(1)an3(1)···an2(1) a1n(1) 
a2n (1) a3n .. . (1)ann(2)

で、b と x は列ベクトル

 (1) b1
 (1) b2  (1) (1)b =b3  , ..  . x1
 x2   x =  x3  , ..  . (3)xn

(1)bnです。これから A と b に操作を加えて値を変えていきますのでその初期値ということで (1) をつけています。
まずはこの連立一次方程式を変換し、A′ を上三角行列とする形に持っていきます。

A′ x = b ′

1(4)

ここで

 ′a11 0
A′ =  0 ..
 .

a′12a′220..
.

a′13a′23a′33..
.

·········..
.

000···

a′1na′2n 
a′3n ,.. . ′ann

 ′b1
 b′2  ′ b′ =  b3  , .. .
b′n(5)

式 (2) から (5) への変換は以下のように行います。a11 がゼロでないとして (ゼロの場合は他の行と順番を入れ替えます)、i 行目の式 (2 ≤ j ≤ n)

ai1 x1 + ai2 x2 + ai3 x3 + · · · + ain xn = bi

(6)

から 1 行目の式に ai1 /a11 掛けたもの

ai1 x1 +a13 ai1ain ai1
ai1
a12 ai1x2 +x3 · · · +xn =
bi
a11
a11
a11
a11(7)を引くことで、j 行目の式から x1 の項を消すことができます。この操作で行列 A と b は

 (1)
a11
 0(2)A = 0 ..
 .

a12(2)
a22(2)a32..
.

a13(2)
a23(2)a33..
.

·········..
.

(1) a1n(2) a2n (2) a3n ,... 0an2(2)an3(2)

···

ann

(1)

(1)

 (1) b1
 (2) b2 
 (2) (2)b =b3 
 .. 
 . (2)

(8)

(2)

bn

行列要素は 1 行目は不変で、2 行目以降は A(1) のものと値が変わります。
(1) (1)ai1 a1j(2)(1),
aij =aij −(1)
a11(1) (1)a b(2)(1)bi =bi − i1 (1)1a11(2 ≤ i, j ≤ n)(9)(10)(3)(3)(3)となる。続いて、この行列 A(2) と列ベクトル b(2) に対して、同様に a32 , a42 , · · · an2 がゼロとなるように変換し、行列 A(3) と列ベクトル b(3) を作ります。

 (1)
a11
 0

(3)
A =
 0
 ..
 .

a12(2)
a220..
.

a13(2)
a23(3)a33..
.

·········..
.

(1) 
a1n(2) a2n (3) a3n ,.. . 00

an3(3)···ann(1)(1)(3)

 (1) b1
 (2) b2 
 (3) (3)b =b3  ..  . (11)

(3)

bn

ここで同様に各成分は (3) が付いているもののみ(2) (2)(3)(2)

aij =aij −ai2 a2j(2)(12)

a22(2) (2)

(3)

bi(2)

=bi−

2

ai2 b2(2)a22(13)

となります。さらに (A(4) , b(4) ), (A(5) , b(5) )、と繰り返してゆくと、最終的に A(n) として上三角行列の形になったものが得られます。プログラムを書くときは A(k) と b(k) を用いて A(k+1) と b(k) を計算するときにその行列要素は(k) (k)aik akj(k+1)(k)
aij=aij −,(k)akk(k) (k)a b(k+1)(k)bi=bi − ik (k)kakk(k ≤ i, j ≤, n)(14)(15)で与えられる、という式でまとめて書くことができます。A(n) までこれらの操作を行うと最終的に

 (1)
a11

 0

(n)
A =
 0
 ..
 .

a12(2)
a220..
.

a13(2)
a23(3)a33..
.

·········..
.

(1) a1n(2) a2n (3) 
a3n ,.. . 

 (1) b1
 (2)  b2 
 (3) b= b3  , .. 
 . 000···ann(n)

bn(1)

(1)

(16)

(n)となり、A(n) は上三角行列となります。この操作を前進消去 (forward elimination) と言います。

n 行目の方程式は(n)a(n)
nn xn = bn(17)となり解は(n)xn =bn(18)(n)

annとして求まります。その一つ手前の n − 1 行目の式は
(n−1)

(n−1)

(n−1)

an−1n−1 xn−1 + an−1n xn = bn−1

(19)ですが、式 (18) を使うと xn−1 が求められます。

xn−1 =1

(n−1)

(n−1)
an−1n−1(n−1)(bn−1 − an−1n xn )(20)

このように、n 行目から始めて順番に xk を求めることができます。xk+1 , · · · xn の解が求まっているとき、xkは

xk =1(k)akk[(k)bk −n∑](k)aki xi(21)

i=k+1で与えられます。これは後退代入 (backward substitution) と呼ばれ、n 行目から 1 行目まで順番に行うことで x のすべての成分を求めることができます。

1.2 Gauss-Jordan の消去法
同じく連立一次方程式を解く方法ですが、行列 A を変換するときに上三角行列にせずに単位行列に変換する方法です。

3

(1)

はじめに A(1) からスタートし、次のステップでは一行目の方程式の x1 の係数を 1 とするため方程式を a11で割り、それを用いて他の方程式から x1 の項をすべて消去します。

1

0
A(2) = 0 ..
.

a12(2)
a22(2)a32..
.

a13(2)
a23(2)a33..
.

·········..
.

(2) a1n(2) a2n (2) a3n ,.. . 

0an2(2)an3(2)···

ann

(2)

(2)

 (2) b1
 (2) b2 
 (2) b(2) = b3  ..  . 

(22)

(2)

(2)

bn

ここで(1)a1j(2)
a1j = (1) ,a11(1 ≤ j ≤ n)(1) (1)ai1 a1j(2)(1)
aij =aij −(1)
a11(1)b(2)b1 = 1(1)a11(1) (1)a b(2)(1)bi =bi − i1 (1)ia11(23)(i ̸= 1)

(24)
(25)
(26)

(2)

次のステップでは 2 番目の方程式の x2 の係数を 1 とするため a22 で方程式を割り、それを用いて 1 番目を含むほかのすべての式から x2 の項を消去します。

10
(3)
A =
0
 ..
.

010..
.

a13(3)
a23(3)a33..
.

·········..
.

00

an3(3)···(3)(3) a1n(3) a2n (3) a3n ,... (3)ann

 (3) b1
 (3) b2 
 (3) (2)b =b3 
 ..  . (3)bn(27)

ここで(2)
a2j(3)
a2j = (2) ,a22(28)(2) (2)(3)(2)aij =aij −ai2 a2j(2)a22(i ̸= 2)

(29)

(2)

b(3)
b2 = 2(2)a22(30)(2) (2)

(3)

bi(2)

=bi−

ai2 bi(2)a22

4(i ̸= 2)(31)

同様に、A(k) から A(k+1) を作るときは k 行目の対角要素が 1 にできて、他の行列要素は
(k)
akj
(k+1)
akj
= (k) ,akk

(32)(k) (k)

(k+1)

aij

aik akj(k)=aij −(k)akk(i ̸= k)

(33)

(k)

(k+1)

bkb= k(k)akk(34)(k) (k)

(k+1)

bi(k)

=bi−

aik bi(i ̸= k)(k)akk(35)となります。n + 1 までこの操作を繰り返すと A(n+1) は単位行列となっているため、解は(n+1)

xi = bi(36)で与えられます。

2 反復法
反復法は連立方程式の問題に限らず数値計算で広く使われる方法です。求めたい解を一発で計算するのではなく、近似解として適当な値からスタートし、その値を更新してゆくことで真の解に収束させる方法です。ここでは Jacobi 法、Gauss-Seidel 法と SOR 法を紹介します。

2.1 Jacobi 法
連立方程式 Ax = b で、i 行目の式は

ai1 x1 + ai2 x2 + · · · ain xn = bi

(37)

ですが、xi のみを左辺に残すとbi −∑aij xj

j(̸=i)

xi =(38)

aiiとなります。これは、x1 , x2 , · · · xn (xi を除く) から xi を計算するという式になっており、反復 k 回目の近似(k)(k)(k)

解 x(k) = (x1 , x2 , · · · xn ) を右辺に代入し、

bi −(k+1)

xi

=

∑(k)

aij xjj(̸=i)(1 ≤ i ≤ n)aii(39)のように、k + 1 回目の近似解を計算することができます。これが Jacobi 法で、反復の結果収束すれば、真の解に到達します。収束したかどうかは、|x(k) − x(k+1) | の絶対値の大きさで判定します。収束は次のGauss-Seidel 法より遅いですが、線形性を保っているため、並列計算に向いたアルゴリズムです。

52.2 Gauss-Seidel 法
Jacobi 法を改良したもので、反復で解を更新するのは同じですが、Jacobi 法では x(k) を用いて x(k+1) を一斉に求めていました。Gauss-Seidel 法では利用できる最新のものを用います。つまり、x1 については Jacobi法と同じ

b1 −(k+1)

x1(k+1)

を用いますが、x2n∑(k)

a1j xjj=2=(40)a11(k+1)を計算する際には x1 についてはすでに計算した x1がありますのでそれを使います。
(k+1)

b2 − a21 x1(k+1)

x2(k+1)x3(k+1)

では x1(k+1)と x2(k)を x1(k)(k+1)

同様に xi(k+1)

では x1(k+1), x2(k+1)− a32 x2(k)a3j xj(42)

a33
(k+1)

, · · · xi−1を用いて

bi −(k+1)n∑−j=4=xi(41)の代わりに用います。
(k+1)

(k+1)

(k)

a2j xja22b3 − a31 x1x3n∑j=3=と x2−i−1∑(k+1)

aij xj(k)aij xjj=i+1

j=1=n∑−(43)

aiiとなります。Jacobi 法より早く解に収束しますが、並列計算には向きません。

2.3 SOR 法
逐次加速緩和法 (Successive Over-Relaxation 法) は Gauss-Seidel 法にさらに加速パラメターを導入して高速化を図ったものです。SOR 法ではn∑b1 −(k+1)x1(k)a1j xjj=2(k)=(1 − ω)x1 + ω,a11(k+1)

b2 − a21 x1(k+1)

x2−(44)n∑(k)

a2j xjj=3(k)=(1 − ω)x2 + ω(45)

a22..
.

(46)
bi −(k+1)xi(k)=(1 − ω)xi+ωi−1∑(k+1)

aij xjj=1

j=i+1aii6−n∑(k)

aij xj(47)

(k+1)とします。右辺の第二項は Gauss-Seidel 法での表式と全く同じですが、含まれている xiは SOR 法で決めたものを用います。加速パラメター ω は 0 < ω < 2 の範囲のものである必要があります。収束を加速させるには 1 < ω < 2 の範囲で最適な値を選ぶ必要があります。ω = 1 のときは Gauss-Seidel 法と一致します
し、ω < 1 の場合は Gauss-Seidel 法より収束は遅くなりますが、挙動が安定するため、Gauss-Seidel 法で解けない問題も解ける可能性があります。

3 演習問題
連立一次方程式

51−21−141−1

   
1 2x118x2  −92 1  =  −6 2 x3   14 
3 4x44を考える。

(22) Gauss 消去法でこの連立一次方程式を解け。
(23) Gauss-Jordan 消去法でこの連立一次方程式を解け。
(24) Jacobi 法 (反復法) でこの連立一次方程式を解け。
(25) Gauss-Seidel 法/SOR 法 (反復法) でこの連立一次方程式を解け。

7