4 連立方程式の解法

２．４連立方程式の解法
大きく分類して，以下の方法がある。
①直接法
②間接法（別名：反復法）
［直接法］
①ガウス（Gauss)の消去法
②ガウス・ジョルダン（Gauss-Jordan）法
（別名：掃出し法）
③ＬＵ分解法
Johann Carl Friedrich Gaus(1777-1855)
［間接法］（別名反復法）
①ヤコビの反復法（連立変位法）（ヤコブ法ともいう）
②ガウス・ザイデル（Gauss-Seidel）法
②ＳＯＲ法
（Successive Over-Relaxation method:逐次過緩和法）
Carl Gustav Jacob Jacobi（1804年-1851）
ちょっと横道にそれて…
同一人物と間違えそうな二人のジョルダン
①ドイツの測地学者 Wilhelm Jordan(1819-1904)
ガウス・ジョルダン法の発案者
②フランスの数学者 Camille Jordan (1771-1821)
ジョルダン標準形の発案者
ちょっと横道にそれて…
ルートヴィヒ・ザイデル（Philipp Ludwig von Seidel）ってどんな人？
（1821年8月24日 - 1896年8月13日）
ツバイブリュッケン生まれ。ドイツの数学者、光学者、天文学者。
光学の研究を行いザイデル収差で知られる。
1846年：ミュンヘン大学で望遠鏡の光学系の研究で博士号。
1847年：ミュンヘン大学教授。
1856年：収差理論に関する著書を出版。
1874年：連立一次方程式の解法である
ガウス＝ザイデル法を発表。
1879年～1882年までボーゲンハウゼン天文台の所長。
（１）ガウス（Gauss）の消去法
前進消去（forward elimination ）と
後退代入（backward substitution）に分かれる。
［前進消去］
6 5 4 8 


12
13
10
16


18 21 17 27
２行目－１行目×２→２行目
３行目－１行目×３→３行目
6 5 4 8 


0
3
2
0


0 6 5 3
３行目－２行目×２→３行目
［後退代入］
x3  3
b2  2  x3 0  2  3
x2 

 2
3
3
b1  5  x2  4  x3 8  5  (2)  4  3
x1 

1
6
6
6 5 4 8 


0
3
2
0


0 0 1 3
ガウスの消去法の考え方
一般化
［前進消去］
処理効率をあげるため共通の係数を先に計算(対角項を１にしてもよい）

aik
akk
(i  k  1, , n)
aij  aij    akj
bi  bi    bk
［後退代入］
( j  k  1, , n)
(i  j  1, , n)
対角項を１にすれば，除算は必要ない。
bn
xn 

ann
bk 
xk 
n
 a x
j  k 1

akk
kj
j
(k  n  1,  ,1)
 a12

a11
 0 a
22




0
0
 a1n   x1   b1 
 a2 n   x2  b2 

       
   

 ann   xn  bn 
ＶＢでプログラムを作る
①シートの定義（なおＡ＊ｘを計算する部分を用意しておく）
配列は，Ｎ×（Ｎ＋１）とし，最後にＢの値を入れる。
②データの設定部分
Public A() As Double
Public B() As Double
Public C() As Double
Public ESP As Double
Private Function データ設定()
N = 3
N1 = N + 1
ReDim A(N, N1)
EPS = 0.000001
N1 = N + 1
With Worksheets("Sheet1")
For i = 1 To N
For J = 1 To N1
A(i, J) = Val(.Cells(i, J))
Next J
Next i
End With
データ設定 = N
End Function
③結果設定部分とボタンイベントハンドラ
Private Sub 結果表示(N)
N1 = N + 1
With Worksheets("Sheet1")
For i = 1 To N
‘ 以下3行は１行に書くこと
.Cells(i, 6) = Right(Space(12) &
Format(A(i, N1), "0.00000E+00"),
12)
Next i
End With
End Sub
Sub ボタン1_Click()
N = データ設定()
If N >= 2 Then
R = 消去法(A, N, ESP)
結果表示 N
End If
End Sub
④処理本体（１）
Public Function 消去法(A, N, ESP) As Boolean
On Error GoTo エラー
消去法 = False
N1 = N + 1
For K = 1 To N
' 前進消去
If Abs(A(K, K)) < ESP Then
消去法 = True
MsgBox “解を求めることができません"
Exit Function
End If
For I = K + 1 To N
ALFA = A(I, K) / A(K, K)
For j = K + 1 To N1
A(I, j) = A(I, j) - ALFA * A(K, j)
Next
Next
Next K
（続きあり）
④処理本体（２：後退代入）
For I = N To 1 Step -1 ' 後退代入
T = A(I, N + 1)
For j = I + 1 To N
T = T - A(I, j) * A(j, N + 1)
Next
A(I, N + 1) = T / A(I, I)
Next
復帰: Exit Function
エラー:
消去法 = True
MsgBox Error(Err)
Resume 復帰
End Function
（２）ガウス－ジョルダン（Gauss-Jordan）法
（掃出し法）
消去方法
1
2

1

3
1 1 1 10
1 3 2 21
3 2 1 17 

2 1 1 14
１行目/1
２行目－１行目×２
３行目－１行目
４行目－１行目×３
1
0

0

0
11 
 1 
9 

0  3  2  17 
0 2
1 1
0 3
1
0
0
１行目－３行目×（２／３）
２行目－３行目×（－１／３）
３行目／３
４行目－１行目×（－３／３）
1
1
10 
1 1
0  1 1

0
1


0 2
1
0
7 


0

1

2

2

16


１行目－２行目×（－１）
２行目／（－１）
３行目－２行目×２×（－
１）
４行目－１行目×（－１）
1
0

0

0





0 0  2  8
0 0
1 0
0 1
1
0
0
5
2
3
１行目－４行目×（－１／２）
２行目－４行目×０
３行目－４行目×０
４行目／（－２）
1
0

0

0
11 
 1 
9 

0  3  2  17 
0 2
1 1
0 3
1
0

0

0
1
0
0
0 0 0 1
1 0 0 2
0 1 0 3

0 0 1 4
④処理本体（１）
Public Function 掃出法(A, N, ESP) As Boolean
On Error GoTo エラー
掃出法 = False
N1 = N + 1
For K = 1 To N
If Abs(A(K, K)) < ESP Then
掃出法 = True
MsgBox "行列式の解を求めることができません"
Exit Function
End If
K1 = K + 1
For J = K1 To N1
A(K, J) = A(K, J) / A(K, K)
Next
（続きあり）
④処理本体（２）
For i = 1 To N
If i <> K Then
For J = K1 To N1
A(i, J) = A(i, J) - A(i, K) * A(K, J)
Next
End If
Next
Next K
復帰: Exit Function
エラー:
掃出法 = True
MsgBox Error(Err)
Resume 復帰
End Function
（３）ピボット選択付きのガウスの消去法
最大値を対角項にもってくるよう行を入替える。
（行を入替えても解は変わらないことに注意）
 0 5 4   x1   8 
12 0 10  x   16 

 2   
18 21 0   x3  27
5 x2  4 x3  8
12 x1
18 x1  21x2
両者とも
同じ連立方程式で
あることに注意
18 21 0   x1  27
 0 5 4 x    8 

 2   
12 0 10  x3  16 
18 x1  21x2
 10 x3  16
 27
 27
5 x2  4 x3  8
12 x1
 10 x3  16
④処理本体（１）
Private Function 消去法(A, N, ESP) As Boolean
On Error GoTo エラー
消去法 = False
N1 = N + 1
For K = 1 To N
K1 = K + 1
Max = Abs(A(K, K)): IR = K
' ピボット入替え
If K <> N Then
For i = K1 To N
If Abs(A(i, K)) >= Max Then
Max = Abs(A(i, K))
IR = i
End If
Next i
End If
（続きあり）
ピボット入替え部分
④処理本体（２）
If Max < ESP Then
消去法 = True
MsgBox "行列式の解を求めることができません"
Exit Function
End If
If IR <> K Then
For j = K To N1
AA = A(K, j): A(K, j) = A(IR, j): A(IR, j) = AA
Next j
End If
W = A(K, K)
' 前進消去
For j = K To N1
A(K, j) = A(K, j) / W
Next j
（続きあり）
④処理本体（３）
If K = N Then Exit For
For i = K1 To N
For j = K1 To N1
A(i, j) = A(i, j) - A(i, K) * A(K, j)
Next j
Next i
Next K
For K = N - 1 To 1 Step -1
' 後退代入
For j = K + 1 To N
A(K, N1) = A(K, N1) - A(K, j) * A(j, N1)
Next j
Next K
復帰: Exit Function
エラー:
消去法 = True
MsgBox Error(Err)
Resume 復帰
End Function
（４）逆行列
定義から
18 21 0  b11 b12 b13  1 0 0
 0 5 4  b b
  0 1 0 
b

  21 22 23  

12 0 10 b31 b32 b33  0 0 1
次のような方程式を解き，bij を求めることに帰着する。
18 21 0  b11  1
 0 5 4  b   0

  21   
12 0 10 b31  0
18 21 0  b12  0
 0 5 4  b   1

  22   
12 0 10 b32  0
18 21 0  b13  0
 0 5 4  b   0

  23   
12 0 10 b33  1
すなわち，Ｎ×（２N）の行列を用意し，連立方程式を解けばよい。
ＶＢでプログラムを作る
①シートの定義（なおＡ＊Ａ－１を計算する部分を用意しておく）
配列サイズは，Ｎ×２Ｎとし，Ｎ＋１から２Ｎの範囲に
対角項のみ１，それ以外に０を入れてから掃出し法を
適用する。
②データの設定部分
Public A() As Double
Public B() As Double
Public C() As Double
Function データ設定()
N=3
ReDim A(N, N * 2)
With Worksheets("Sheet1")
For i = 1 To N
For j = 1 To N
A(i, j) = Val(.Cells(i, j))
Next j
Next i
End With
データ設定 = N
End Function
③結果設定部分とボタンイベントハンドラ
Sub 結果表示(N)
With Worksheets("Sheet1")
For k = 1 To N
For i = 1 To N
‘ 以下２行は１行に記述する
.Cells(i, k + 4) = Right(Space(12) & Format(A(i, k
+ N), "0.00000E+00"), 12)
Next i
Next
.Select
End With
End Sub
Sub ボタン1_Click()
N = データ設定()
ESP = 0.00001
R = 掃出法(A, N, ESP)
結果表示 N
End Sub
④処理本体（１）
Public Function 掃出法(A, N, ESP) As Boolean
'On Error GoTo エラー
For i = 1 To N
For j = 1 To N
A(i, j + N) = 0
If i = j Then A(i, j + N) = 1
Next
Next
掃出法 = False
N2 = N * 2
For k = 1 To N
If Abs(A(k, k)) < ESP Then
掃出法 = True
MsgBox "行列式の解を求めることができません"
Exit Function
End If
（続きあり）
④処理本体（２）
S = 1 / A(k, k)
For j = k + 1 To N2
A(k, j) = A(k, j) * S
Next
For i = 1 To N
If i <> k Then
S = A(i, k)
For j = k + 1 To N2
A(i, j) = A(i, j) - S * A(k, j)
Next
End If
Next
Next
復帰: Exit Function
エラー:
掃出法 = True
MsgBox Error(Err)
Resume 復帰
End Function
（５）ＬＵ分解による方法
［例］
 6 5 4   x1  1 0 0 6 5 4  x1   8 
12 13 10  x   2 1 0 0 3 2  x   16 

 2  

 2   
18 21 17  x3  3 2 1 0 0 1  x3  27
6 5 4  x1   y1 
1 0 0  y1   8 
0 3 2  x    y  と置くと 
  y   16 
2
1
0
2
2

   

 2   
0 0 1  x3   y3 
3 2 1  y3  27
y1  8, 2 y1  y2  16, 3 y1  2 y2  y3  27
y2  16  2 y1  16  16  0,
y3  27  3 y1  2 y2  27  24  0  3
6 5 4  x1  8
0 3 2   x   0 

 2   
0 0 1  x3  3
この演算の意味
元の行列に施した処理を
列ベクトルにも施す。
下三角行列は，
元の行列を
上三角行列に変換する行列の
逆行列である。
x3  3, 3 x2  2 x3  0, 6 x1  5 x2  4 x3  8
0  2 x3  2  3

 2
3
3
8  5 x2  4 x3 8  5  (2)  4  3 6
x1 

 1
6
6
6
x2 
確認
（例）下三角行列の逆行列を元の行列に乗じると上三角行列になる
0 0 1 0 0  1
0
0  1 0 0
1
  2 1 0   2 1 0    2  2
  0 1 0 
0

1
0

0


 
 

 1  2 1 3 2 1 1  4  3 0  2  2 0  0  1 0 0 1
すなわち
1 0 0 
 2 1 0


3 2 1
の逆行列
0 0
1
  2 1 0


 1  2 1
を元の行列に
乗じると
0 0  6 5 4  
6
5
4
1
 6 5 4
 2 1 0 12 13 10    12  12
  0 3 2

10

13

8

10


 
 

 1  2 1 18 21 17 6  24  18 5  26  21 4  20  17 0 0 1
ＶＢでプログラムを作る
①シートの定義（行列のＬＵ分解で作ったシートに追加する）
（なおＡ＊Ａ－１を計算する部分を用意しておく）
この部分を追加
②データの設定部分
Private A(3, 3) As Double
Private L(3, 3) As Double
Private U(3, 3) As Double
Private Y(3) As Double
Private X(3) As Double
Private B(3) As Double
Private Sub データ設定()
With Worksheets("Sheet1")
For i = 1 To 3
For j = 1 To 3
A(i, j) = .Cells(i + 1, j + 1)
Next
Next
For i = 1 To 3
B(i) = .Cells(i + 1, 16)
Next
End With
End Sub
③結果設定部分とボタンイベントハンドラ
Private Sub 結果設定(L, U, Y, X)
With Worksheets("Sheet1")
For i = 1 To 3
For j = 1 To 3
.Cells(i + 1, j + 6) = L(i, j)
.Cells(i + 1, j + 11) = U(i, j)
Next
Next
For i = 1 To 3
.Cells(i + 1, 17) = Y(i)
.Cells(i + 1, 18) = X(i)
Next
End With
End Sub
Sub ボタン3_Click()
データ設定
LU分解による連立方程式 A, L, U, B, Y, X, 3
結果設定 L, U, Y, X
End Sub
④LU分解処理（2.3に同じ。再掲）
Private Sub LU分解(A, L, U, N)
For i = 1 To N
For j = 1 To N
L(i, j) = 0
U(i, j) = 0
Next
Next
L(1, 1) = 1
For j = 1 To N
U(1, j) = A(1, j)
Next
For i = 2 To N
U(i, 1) = 0
L(1, i) = 0
L(i, 1) = A(i, 1) / U(1, 1)
Next
For i = 2 To N
L(i, i) = 1
T = A(i, i)
For k = 1 To i
T = T - L(i, k) * U(k, i)
Next
U(i, i) = T
For j = i + 1 To N
U(j, i) = 0: L(i, j) = 0
T = A(j, i)
For k = 1 To i
T = T - L(j, k) * U(k, i)
Next
L(j, i) = T / U(i, i)
T = A(i, j)
For k = 1 To i
T = T - L(i, k) * U(k, j)
Next
U(i, j) = T
Next
Next
End Sub
⑤LU分解した結果を使って連立方程式を解く処理
Sub LU分解による連立方程式(A, L, U, B, Y, X, N)
LU分解 A, L, U, N
For i = 1 To N
T = B(i)
For k = 1 To i - 1
T = T - L(i, k) * Y(k)
Next
Y(i) = T
Next
For i = N To 1 Step -1
T = Y(i)
For k = i + 1 To N
T = T - U(i, k) * X(k)
Next
X(i) = T / U(i, i)
Next
End Sub
（６）ピボット選択によるＬＵ分解
①クラウト法によるＬＵ分解は分かりやすいが，ピボット選択（対角
項の絶対値を最大値にとる）が適用できない。
②ガウスの消去法における前進消去は，上三角行列を作成する
過程と同じ。したがって，ガウスの消去法における前進消去を用
いてＬＵ分解することが可能。
考え方の基本
①ガウスの消去法の前進消去は，注目する i 列に対して
a
g ki   ki (k  i ) それ以外は，対角項=1, 非対角項=0
aii
という変換行列で表すことができる。
［例］
0 0 18 21 17 18 21 17 
 1
 12 18 1 0 12 13 10   0  1  4 3


 

  6 18 0 1  6 5 4   0  2  5 3
②ピボット選択は置換行列で表現できる
0 0 1  a11 a12
0 1 0 a

  21 a22
1 0 0 a31 a32
a13  a31 a32
a23   a21 a22
a33   a11 a12
a33 
a23 
a13 
ガウスの消去法によるＬＵ分解の考え方（１）
（以下，式の展開を覚える必要はない。どのような手順で計算されるのかを理解するこ
と）
ピボット選択を行うガウスの消去法の前進消去は，置換行列 Pi (i  1,, n  1)
および変換行列 Gi (i  1,, n  1) を用いると
Gn 1 Pn 1 G1 P1 Ax  Gn 1 Pn1 G1 P1 x
右辺の Gn 1 Pn 1 G1 P1 を計算すると

 Gn 1 Pn 1Gn  2 Pn  2 G1 P1 P1
1
P
1
 Gn 1  Pn 1Gn  2 Pn  2  P2G1 I P1

P
P
P

1
 Gn 1 Pn 1Gn  2  P2G1 IP2 P2 P1
 Gn 1
 Gn 1
 Gn 1
同様にして

n 1G n  2  P2G1 IP2
1
P P
2

P P P P
 P P G IP P P P P
n 1G n  2 Pn  2  P2G1 IP2
n 1G n  2 Pn  2
 
1
 Gn1 Pn1  P3 P2G1 IP2
1
1
3
3
1
3
1
2
1
2
P P
1
3
3
1
2
1
1
3
3

2
1
 Pn11 Pn1  P3 P2 P1
ガウスの消去法によるＬＵ分解の考え方（１）
（以下，式の展開を覚える必要はない。どのような手順で計算されるのかを理解すること）

 
Gn1 Pn1  P3 P2G1 IP2
1
P P
1
3
3
1

 Pn11 Pn1  P3 P2 P1
この式を計算するために漸化式の形で表現する。
H0  I,
H i 1  PiGi 1 H i  2 Pi
1
(i  2,  , n  2),
n 1
P   Pk
k 1
ガウスの消去法によるＬＵ分解の考え方（２）
さらに A0  A,
Ai  Gi Pi Ai 1
(i  1,, n  1) とすると
An 1 は，前進消去の最後であるため，上三角行列となる。したがって，
Gn 1 Pn 1 G1 P1 Ax  Gn 1 Pn1 G1 P1 x
は，次のように表現できる。
Ux  Gn 1 H n  2 Px
すなわち，
下三角行列は
Gn1 H n2 1Ux  Px
L  Gn1 H n2 
1
となる。
ガウスの消去法によるＬＵ分解のアルゴリズム
ピボット選択を行う前進消去を，置換行列 Pi (i  1,, n  1)
および変換行列Gi (i  1,, n  1) を用いて表すと，
① P1 , G1 を求め， A1  P1 G1 A ， H 0  I とする。
（ピボット選択を行う前進消去を1桁だけ実行したことと同じ）
P1 の替わりに置換した元のピボット位置を記憶してもよい。
ガウスの消去法によるＬＵ分解のアルゴリズム
② 2  i について Pi , Gi を求め
Ai  Pi Gi Ai 1 , Hi 1  Pi Gi 1 Hi 2 Pi
（置換行列の性質から Pi  Pi
1
1
であることに注意）
③ステップ②を i  n 1 まで繰返し， i  n 1 のときの
An1  Pn1 Gn1 An2 が上三角行列，
L  (Gn 1 H n  2 ) 1 が下三角行列となる。
ＶＢでプログラムを作る
①シートの定義
置換行列があることに注意
②行列演算等，共通的な処理が多いので，できる限り共通ルーチン化する。
③連立方程式の解法は，クラウト法による解法を参考にして作成してみよう。
なお，列ベクトルに置換行列を乗じてから前進代入，後退代入を行うこと。
③データの設定部分
Private A(3, 3) As Double
Private U(3, 3) As Double
Private L(3, 3) As Double
Private multP(3, 3) As Double
Private Sub データ設定()
With Worksheets("Sheet1")
For i = 1 To 3
For J = 1 To 3
A(i, J) = .Cells(i
Next
Next
End With
End Sub
' 対象とする配列
' 上三角行列
' 下三角行列
'置換行列の積
+ 1, J + 1)
④結果設定とボタンイベントハンドラ
Sub 結果設定()
With Worksheets("Sheet1")
For i = 1 To 3
For J = 1 To 3
.Cells(i + 1, J + 6) = L(i, J)
.Cells(i + 1, J + 11) = U(i, J)
.Cells(i + 1, J + 16) = multP(i, J)
Next
Next
End With
End Sub
Sub ボタン2_Click()
TraceMode = True
データ設定
ピボット選択LU分解 A, L, U, multP, 3
結果設定
End Sub
⑤共通ルーチン（１）
Private Function 最大対角項(A, K, N) As Integer
ID = K
For J = K + 1 To N
If Abs(A(ID, K)) < Abs(A(J, K)) Then ID = J
Next
最大対角項 = ID
End Function
Private Sub 行入替え(K1, K2, P, N) '行入替え
For i = 1 To N
T = P(K1, i): P(K1, i) = P(K2, i): P(K2, i) = T
Next
End Sub
⑤共通ルーチン（２）
Private Sub 単位行列設定(H, N) ' 単位行列の設定
For i = 1 To N
For J = 1 To N
H(i, J) = 0
If i = J Then H(i, J) = 1
Next
Next
End Sub
Private Sub 置換行列(K1, K2, P, N) '置換行列Pの設定
単位行列設定 P, N
If K1 = K2 Then Exit Sub
行入替え K1, K2, P, N
End Sub
⑤共通ルーチン（３）
Private Sub 行列乗算(A, B, C, N1, N2) '行列乗算
For i = 1 To N1
' C = A * B
For J = 1 To N2
C(i, J) = 0
For K = 1 To N1
C(i, J) = C(i, J) + A(i, K) * B(K, J)
Next
Next
Next
End Sub
⑤共通ルーチン（４）
Private Sub 変換行列(A, K, G, N) '変換行列Gの設定
単位行列設定 G, N
AKK = A(K, K)
For i = K + 1 To N
G(i, K) = -A(i, K) / AKK
Next
End Sub
⑤共通ルーチン（５）
Public Function 掃出法(A, N, ESP) As Boolean
'On Error GoTo 掃出エラー
For i = 1 To N
For J = 1 To N
For i = 1 To N
A(i, J + N) = 0
If i <> K Then
If i = J Then A(i, J + N) = 1
S = A(i, K)
Next
For J = K + 1 To N2
Next
A(i, J) = A(i, J) - S * A(K, J)
掃出法 = False
Next
N2 = N * 2
End If
For K = 1 To N
Next
If Abs(A(K, K)) < ESP Then
Next
掃出法 = True
掃出復帰: Exit Function
MsgBox "解を求めることができません"
掃出エラー:
Exit Function
掃出法 = True
End If
MsgBox Error(Err)
S = 1 / A(K, K)
Resume 掃出復帰
For J = K + 1 To N2
End Function
A(K, J) = A(K, J) * S
Next
⑤共通ルーチン（６）
Private Sub 掃出結果設定(A, R, N)
For i = 1 To N
For J = 1 To N
R(i, J) = A(i, J + N)
Next
Next
End Sub
Private Sub 配列複写(A, R, N)
For i = 1 To N
For J = 1 To N
R(i, J) = A(i, J)
Next
Next
End Sub
⑥処理本体（初期設定）
Private Sub ピボット選択LU分解(A, L, U, multP, N)
Dim G() As Double ' 置換行列
Dim GB() As Double ' Gn-1を保存するための作業領域
Dim P() As Double ' 置換行列
Dim H() As Double ' 下三角行列用行列 Hi
Dim C() As Double ' 乗算用作業領域1
Dim CC() As Double ' 乗算および逆行列計算用作業領域
ReDim G(N, N), P(N, N), H(N, N), C(N, N), CC(N, 2 * N), GB(N, N)
単位行列設定 H, N
' H0を単位行列とする
単位行列設定 multP, N
' multPを単位行列とする
（実行部分の続きあり）
⑥処理本体（実行部分）
For i = 1 To N - 1
ID = 最大対角項(A, i, N) ' 最大対角項の選択
置換行列 i, ID, P, N
' 置換行列Piの設定
行列乗算 P, A, C, N, N
' C=Pi*Ai-1
変換行列 C, i, G, N
' 前進代入の変換行列Giを求める
行列乗算 G, C, A, N, N
' Ai = Gi*C = Gi*Pi*Ai-1
If i > 1 Then
' (i=1のときGi-1がない)
行列乗算 GB, H, C, N, N '
Hi = Pi*Gi-1*Hi-2*Pi^(-1)
行列乗算 P, C, CC, N, N '
(Pi^(-1) =Piに注意）
行列乗算 CC, P, H, N, N
End If
行列乗算 P, multP, C, N, N ' Pmult=Pi*Pmult
配列複写 C, multP, N
配列複写 G, GB, N
' Gi-1を保存
Next
行列乗算 G, H, CC, N, N
' CC=Gn-1*Hn-2
ESP = 0.0001
R = 掃出法(CC, N, ESP)
' CCの逆行列を計算して
掃出結果設定 CC, L, N
' 下三角行列とする。
配列複写 A, U, N
' An-1を上三角行列とする
End Sub
（７）ピボット選択によるＬＵ分解の改善
①置換行列は結局のところ，行の入替えである。ピボットの位置だ
けを持っていれば可能。
②下三角行列の対角項は１に決まっている。
③上三角行列の左下非対角項は０，下三角行列の右上非対角項
は０である。
⑤したがって，左下非対角項には下三角行列，対角項と右上非対
角項には上三角行列の各項を入れることにすれば，配列サイズ
を削減できる。
ＶＢでプログラムを作る
①シートの定義
②データの設定部分
Public A() As Double
Public B() As Double
Public C() As Double
Public X() As Double
Public LU() As Double
Public IPivot() As Integer
Function データ設定()
N = 4
ReDim A(N, N), LU(N, N), IPivot(N)
With Worksheets("Sheet1")
For I = 1 To N
For J = 1 To N
A(I, J) = Val(.Cells(I + 1, J + 1))
Next
Next
End With
データ設定 = N
End Function
③結果の設定とイベントハンドラ
Sub 結果表示(LU, IPivot, N)
With Worksheets("Sheet1")
For I = 1 To N
For J = 1 To N
.Cells(I + 1, J + 7) = LU(I, J)
Next
Next
For I = 1 To N
.Cells(I + 1, 13) = IPivot(I)
Next
End With
End Sub
Sub ボタン1_Click()
N = データ設定()
EPS = 0.001
LU分解 A, LU, IPivot, N, EPS
結果表示 LU, IPivot, N
End Sub
④共通ルーチン
Sub 配列移動(A, B, N)
For I = 1 To N
For J = 1 To N
B(I, J) = A(I, J)
Next
Next
End Sub
⑤処理本体
Sub LU分解(A, LU, IPivot, N, EPS)
配列移動 A, LU, N
For k = 1 To N
IPivot(k) = k
Next
For k = 1 To N
I = k
For J = k + 1 To N
If Abs(LU(IPivot(J), k)) > Abs(LU(IPivot(I), k)) Then I = J
Next
J = IPivot(k): IPivot(k) = IPivot(I): IPivot(I) = J
pivot = LU(IPivot(k), k)
If Abs(pivot) < EPS Then MsgBox (k & "番目の対角項が" & pivot & "です")
For I = k + 1 To N
T = LU(IPivot(I), k) / pivot: LU(IPivot(I), k) = T
For J = k + 1 To N
LU(IPivot(I), J) = LU(IPivot(I), J) - T * LU(IPivot(k), J)
Next
Next
Next
End Sub
⑥確認用乗算と設定（１）
Sub 確認用乗算と表示()
N = 4
ReDim A(N, N), B(N, N), C(N, N), IPivot(N)
With Worksheets("Sheet1")
For I = 1 To N
IPivot(I) = .Cells(I + 1, 13)
For J = 1 To N
If J = I Then
A(I, J) = Val(.Cells(I + 1, J + 7))
B(I, J) = 1
ElseIf J > I Then
A(I, J) = Val(.Cells(I + 1, J + 7))
B(I, J) = 0
Else
A(I, J) = 0
B(I, J) = Val(.Cells(I + 1, J + 7))
End If
Next
Next
（続きあり）
⑥確認用乗算と設定（２）
For I = 1 To N
For J = 1 To N
C(I, J) = 0
For k = 1 To N
C(I, J) = C(I, J) + B(I, k) * A(k, J)
Next
.Cells(I + 1, J + 18) = C(I, J)
Next
Next
End With
End Sub
Sub ボタン2_Click()
確認用乗算と表示
End Sub
ＬＵ分解による連立方程式の解法，逆行列の計算
参考プログラム１（連立方程式の解法）を開く
参考プログラム２（逆行列の計算）を開く
（８）ヤコビ（Jacobi）の反復法(連立変位法）
連立方程式
a11x1  a12 x2    a1n xn  b1
a21x1  a22 x2    a2 n xn  b2

an1 x1  an 2 x2    ann xn  bn
を，次のように並び替えて考える。

 b

x1  b1  a12 x2  a13 x3    a1n xn a11
x2
2



 a21x1  a23 x3    a2 n xn a22

xn  bn  an1 x1  an 2 x2    an ,n 1 xn ann
ヤコビ反復法の手順
①適当な初期値 x1( 0 ) , x2( 0 ) ,  , xn( 0 ) を与え，以下のように計算する。

 b

a
x1(i )  b1  a12 x2(i 1)  a13 x3(i 1)    a1n xn(i 1) a11
x2(i )
( i 1)
( i 1)
( i 1)

a
x

a
x



a
x
2
21 1
23 3
2n n


22

xn(i )  bn  an1 x1(i 1)  an 2 x2(i 1)    an ,n 1 xn(i11) ann
② xk( i )  xk( i 1)
(k  1, , n) のとき収束したものとする。
この判定には色々考えられるが，たとえば以下のようにする。
xn(i )  xn(i 1)
x1(i )  x1(i 1) x2(i )  x2(i 1)



(i )
(i )
(i )
x1
x2
xn
ＶＢでプログラムを作る
①シートの定義
②データの設定部分
Private A(3, 3) As Double
Private B(3) As Double
Private X(3) As Double
Private Sub データ設定()
With Worksheets("Sheet1")
For i = 1 To 3
B(i) = Val(.Cells(i + 1, 6))
For j = 1 To 3
A(i, j) = Val(.Cells(i + 1, j + 1))
Next
Next
End With
End Sub
③結果の設定とイベントハンドラ
Private Sub 結果設定()
With Worksheets("Sheet1")
For i = 1 To 3
.Cells(i + 1, 8) = X(i)
Next
End With
End Sub
Sub ボタン1_Click()
データ設定
EPS = 0.001
numLoop = Jacobi反復(A, B, X, EPS, 1000, 3)
MsgBox "収束回数=" & numLoop
結果設定
End Sub
④処理本体（１）
Private Function Jacobi反復(A, B, X, EPS, iterMax, N) As Integer
Dim Y() As Double: ReDim Y(N)
For i = 1 To N
X(i) = B(i) / A(i, i)
Next
E = EPS * 100: R = 0
Do While E > EPS And R < iterMax
R = R + 1
S = ""
For i = 1 To N
Y(i) = B(i)
For j = 1 To N
If i <> j Then Y(i) = Y(i) - A(i, j) * X(j)
Next
Y(i) = Y(i) / A(i, i)
Next
（続きあり）
④処理本体（２）
dev = 0
For i = 1 To N
EI = (Y(i) - X(i)) / Y(i): dev = dev + Sqr(EI * EI)
X(i) = Y(i)
Next
E = dev
Loop
Jacobi反復 = R
End Function
ヤコビの反復法の発散について
非対角項の絶対値に比べ，対角項の絶対値がかなり大きいとき
は収束するが，同じくらいか小さい場合はなかなか収束しない。
逆に発散することもある。
ただし，工学的には，非対角項に比べ対角項の絶対値が大きい
問題が圧倒的に多い。
［演習］
以下の連立方程式がヤコブの反復法で発散することを確認せよ。
 6 5 4   x1   8 
12 13 10  x   16 

 2   
18 21 17  x3  27
（９）ガウス・ザイデル（Gauss-Seidel）法
ヤコビの反復法における j 行目の計算で xk( i 1) ( k  j ) の替わりに，
既に計算されている xk( i ) (k  j ) を用いて収束を早める。
すなわち，
j 1
x 
(i )
j
bi   a jk xk(i ) 
k 1
n
( i 1)
a
x
 jk k
k  j 1
aii
収束条件は，
x (i )  x (i 1)  
ＶＢでプログラムを作る
①シートの定義（ヤコビの反復法で収束しなかった例を使う）
②データの設定部分（ヤコビの反復法と同じ）
Private A(3, 3) As Double
Private B(3) As Double
Private X(3) As Double
Private Sub データ設定()
With Worksheets("Sheet1")
For i = 1 To 3
B(i) = Val(.Cells(i + 1, 6))
For j = 1 To 3
A(i, j) = Val(.Cells(i + 1, j + 1))
Next
Next
End With
End Sub
③結果の設定とイベントハンドラ（関数呼出し以外はヤコビの反復法と同じ）
Private Sub 結果設定()
With Worksheets("Sheet1")
For i = 1 To 3
.Cells(i + 1, 8) = X(i)
Next
End With
End Sub
Sub ボタン1_Click()
データ設定
EPS = 0.00001
numLoop = GaussSeidel(A, B, X, EPS, 500, 3)
MsgBox "収束回数=" & numLoop
結果設定
End Sub
③処理本体
Private Function GaussSeidel(A, B, X, EPS, iterMax, N) As Integer
Dim Y() As Double: ReDim Y(N)
For i = 1 To N
X(i) = B(i) / A(i, i)
Next
E = EPS * 100: R = 0
Do While E > EPS And R < iterMax
R = R + 1: dev = 0
For i = 1 To N
Y(i) = B(i)
For j = 1 To N
If i <> j Then Y(i) = Y(i) - A(i, j) * X(j)
Next
Y(i) = Y(i) / A(i, i): EI = (X(i) - Y(i))
dev = dev + EI * EI:
X(i) = Y(i)
Next
E = Sqr(dev)
Loop
GaussSeidel = R
End Function
ヤコビの反復法とガウス・ザイデル法のプログラム上の違い
［演習１］両ソースプログラムを比較し，関数宣言以外で異なっている部分を指摘せよ。
Private Function Jacobi反復(A, B, X, EPS, iterMax, N) As Integer
Dim Y() As Double: ReDim Y(N)
For i = 1 To N
Private Function GaussSeidel(A, B, X, EPS, iterMax, N) As Integer
X(i) = B(i) / A(i, i)
Dim Y() As Double: ReDim Y(N)
Next
For i = 1 To N
E = EPS * 100: R = 0
X(i) = B(i) / A(i, i)
Do While E > EPS And R < iterMax
Next
R=R+1
E = EPS * 100: R = 0
For i = 1 To N
Do While E > EPS And R < iterMax
Y(i) = B(i)
R = R + 1: dev = 0
For j = 1 To N
For i = 1 To N
If i <> j Then Y(i) = Y(i) - A(i, j) * X(j)
Y(i) = B(i)
Next
For j = 1 To N
Y(i) = Y(i) / A(i, i)
If i <> j Then Y(i) = Y(i) - A(i, j) * X(j)
Next
Next
dev = 0
Y(i) = Y(i) / A(i, i)
For i = 1 To N
EI = (Y(i) - Y(i))
EI = (Y(i) - X(i)) / Y(i)
dev = dev + EI * EI: X(i) = Y(i)
dev = dev + Sqr(EI * EI)
Next
X(i) = Y(i)
E = Sqr(dev)
Next
Loop
E = dev
GaussSeidel = R
Loop
End Function
Jacobi反復 = R
End Function
ヤコビの反復法とガウス・ザイデル法の収束回数の違い
［演習２］
以下の連立方程式を両手法で解き，収束回数を確認し，比較せよ。
3 1 1  x1  10 
1 5 2  x    21

 2   
1 2 5  x3  30
ヤコビの反復法の収束回数
＝（
）
ガウスザイデル法の収束回数＝（
）
したがって（
短い回数で収束する。
）のほうが
（１０）ＳＯＲ法（Successive Over-Relaxation method）
逐次過緩和法ともいい，以下の反復式を用いる。
他は，ガウス・ザイデル法と同じ。
j 1

(i )
j
x
(i )
j

bi   a jk xk(i ) 
x
k 1
( i 1)
j

 
n
( i 1)
a
x
 jk k
k  j 1
aii
(i )
j
x
( i 1)
j

where 0    2
［演習］
  1 のとき，ガウス・ザイデル法と同じであることを示せ。
ＶＢでプログラムを作る
①シートの定義（ヤコビの反復法で収束しなかった例を使う）
②データの設定部分（加速度パラメータの部分が異なる）
Private A(3, 3) As Double
Private B(3) As Double
Private X(3) As Double
Private Omega As Double
’ 加速度パラメータ
Private Sub データ設定()
With Worksheets("Sheet1")
For i = 1 To 3
B(i) = Val(.Cells(i + 1, 6))
For j = 1 To 3
A(i, j) = Val(.Cells(i + 1, j + 1))
Next
Next
Omega = Val(.Cells(6, 5)) ’ 加速度パラメータ設定
End With
End Sub
③結果の設定とイベントハンドラ
Private Sub 結果設定()
With Worksheets("Sheet1")
For i = 1 To 3
.Cells(i + 1, 8) = X(i)
Next
End With
End Sub
Sub ボタン1_Click()
データ設定
EPS = 0.00001
numLoop = SOR(A, B, X, EPS, 500, Omega, 3)
MsgBox "収束回数=" & numLoop
結果設定
End Sub
③処理本体
Private Function SOR(A, B, X, EPS, iterMax, Omega, N) As Integer
Dim Y() As Double: ReDim Y(N)
For i = 1 To N
X(i) = B(i) / A(i, i)
Next
E = EPS * 100: R = 0
Do While E > EPS And R < iterMax
R = R + 1: dev = 0
For i = 1 To N
Y(i) = B(i)
For j = 1 To N
If i <> j Then Y(i) = Y(i) - A(i, j) * X(j)
Next
Y(i) = Y(i) / A(i, i): EI = (Y(i) - X(i))
dev = dev + EI * EI:
X(i) = X(i) + Omega * EI
Next
E = Sqr(dev)
Loop
SOR = R
End Function
ガウス・ザイデル法とＳＯＲ法のプログラム上の違い
［演習１］両ソースプログラムを比較し，関数宣言以外で異なっている部分を指摘せよ。
Private Function GaussSeidel(A, B, X, EPS, iterMax, N) As Integer
Dim Y() As Double: ReDim Y(N)
For i = 1 To N
X(i) = B(i) / A(i, i)
Next
E = EPS * 100: R = 0
Do While E > EPS And R < iterMax
R = R + 1: dev = 0
For i = 1 To N
Y(i) = B(i)
For j = 1 To N
If i <> j Then Y(i) = Y(i) - A(i, j) * X(j)
Next
Y(i) = Y(i) / A(i, i): EI = (Y(i) - Y(i))
dev = dev + EI * EI: X(i) = Y(i)
Next
E = Sqr(dev)
Loop
GaussSeidel = R
End Function
Private Function SOR(A, B, X, EPS, iterMax, Omega, N) As Integer
Dim Y() As Double: ReDim Y(N)
For i = 1 To N
X(i) = B(i) / A(i, i)
Next
E = EPS * 100: R = 0
Do While E > EPS And R < iterMax
R = R + 1: dev = 0
For i = 1 To N
Y(i) = B(i)
For j = 1 To N
If i <> j Then Y(i) = Y(i) - A(i, j) * X(j)
Next
Y(i) = Y(i) / A(i, i): EI = (Y(i) - X(i))
dev = dev + EI * EI: X(i) = X(i) + Omega * EI
Next
E = Sqr(dev)
Loop
SOR = R
End Function
加速度パラメータによる収束回数の違い
［演習２］
以下の連立方程式で，加速度パラメータを変化させＳＯＲ手法で解いた場合の収束回
数を確認せよ。なお，加速度パラメータ＝０または２のとき収束しないことを確認せよ。
3 1 1  x1  10 
1 5 2  x    21

 2   
1 2 5  x3  30
加速度パラメータ＝１．０のときの収束回数＝（
）
（ガウス・ザイデル法と同じ）
加速度パラメータ＝１．２のときの収束回数＝（
）
加速度パラメータ＝１．５のときの収束回数＝（
）
加速度パラメータ＝１．７のときの収束回数＝（
）
（１１）３重対角な連立方程式の場合
問題によっては，はじめから３重対角な行列式になることが分かっている場合が
ある。このような場合，Ｎ×Ｎの行列を用意する必要はない。
 b1
a
 2
0

0


 0
c1
0
0

b2
c2
0

a 3 b3
c3

0




 a n 1 b n 1
0

0
an
0 
0 
0 

 
c n 1 

b n 
b1
c1 
 0
a

b
c
2
2
2


 a3
b3
c3 







a n 1 b n 1 c n 1 


bn
0 
 a n
ガウスの消去法と同じ方法で，次のような行列に変換。ただし，隣あう行・列要素
だけに処理すればよい。
0  0 
1 c1 0
0 1 c 

0

0
2


0 0 1 c  3  0 


0
0






    0 1 cn 1 


1 
0 0  0 0
ＶＢでプログラムを作る
①シートの定義
②データの設定部分
Private A(6, 3) As Double
Private D(6) As Double
Private X(6) As Double
Function データ設定() As Integer
With Worksheets("Sheet1")
N = 6: AA = ""
For i = 1 To N
A(i, 1) = .Cells(i + 1, i)
A(i, 2) = .Cells(i + 1, i + 1)
A(i, 3) = .Cells(i + 1, i + 2)
D(i) = .Cells(i + 1, 8)
If i = 1 Then A(i, 1) = 0
If i = N Then A(i, 3) = 0
Next
End With
データ設定 = N
End Function
③結果設定部分
Sub 結果設定(N)
With Worksheets("Sheet1")
For i = 1 To N
.Cells(i + 1, 9) = X(i)
Next
End With
End Sub
Sub ボタン1_Click()
N = データ設定()
TriDiag A, D, X, N
結果設定 N
End Sub
④計算部分
Sub TriDiag(A, D, X, N)
Dim CC() As Double
ReDim CC(N)
btemp = A(1, 2): X(1) = D(1) / btemp
For j = 2 To N
CC(j) = A(j - 1, 3) / btemp
btemp = A(j, 2) - A(j, 1) * CC(j)
X(j) = (D(j) - A(j, 1) * X(j - 1)) / btemp
Next
For j = N - 1 To 1 Step -1
X(j) = X(j) - CC(j + 1) * X(j + 1)
Next
End Sub
（Ｘ）演習
（問題１）次の連立方程式を，サンプルプログラムを利用して解け
（場合によっては，配列サイズ等を変更する必要があるので注意すること）
①
2 x1  2 x2  4 x3  2 x4  10
x1  3 x2  2 x3  x4  17
3 x1  x2  3 x3  x4  18
x1  3 x2  4 x3  2 x4  27
x1  3 x2  x3
①
 x5  2 x6  4 x7  1
 x4
2 x2
 2 x6
3x3  x4  x5
 2 x4
2 x2
2 x1  6 x2  2 x3
3
 x7  5
 4 x6  3 x7  1
 6 x5  5 x6  9 x7  3
x1  5 x2  x3  x4  x5  5 x6  6 x7  7
x1
 4 x3  x4  2 x5  2 x6  8 x7  8
（Ｘ）演習
（問題２）回路内に流れる電流を求めよ。なお，R1 = R2 = R3 = R4 = 1Ωとする。
Ａ
Ｃ
Ｂ
I4
I2
Ｅ
I5
I3
I1
I1
Ｄ
Ｆ
I 4+I 5
Ｇ
I5
キルヒホッフの法則およびＢ，Ｃ点に流れる電流の関係から，
R1 I1  R2 I 2
 R2 I 2  R3 I 3  R4 I 4
 E1
（Ａ－Ｂ－Ｆ－Ｅ－Ａ）
0
（Ｂ－Ｃ－Ｇ－Ｆ－Ｂ）
 R4 I 4  R5 I 5   E2
I1
 I2
 I3
I3
 I4
（Ｃ－Ｄ－Ｈ－Ｇ－Ｃ）
0
（Ｂ点）
 I5  0
（Ｃ点）
Ｈ
（Ｘ）演習
R1 = R2 = R3 = R4 = 1Ωだから，次の連立方程式を解くことに帰着する。
1 1 0 0 0   I1   2 
0  1 1 1 0   I   0 

 2   
0 0 0  1 1   I 3    1

   
1

1

1
0
0

I 4   0 
0 0 1  1  1  I 5   0 
あるいは，対角項に非ゼロ項がくるように整理した連立方程式を解く。
1 1 0 0 0   I1   2 
1  1  1 0 0   I   0 

 2   
0  1 1 1 0   I 3    0 

   
0
0
1

1

1

I 4   0 
0 0 0  1 1   I 5   1
1行目
４行目
２行目
５行目
３行目
（Ｘ）演習
Pspiceによる結果は次のとおりである。結果と比較せよ。
（Ｘ）演習
（問題３）下図のように，構造物に分布荷重，集中荷重がかかったときの，節点の偏角と
水平変位は次のような連立方程式で表現される。
x
3
wL3
41   2  x 
L
24 EI
3
 wL3
1  4 2  x 
L
24 EI
4
PL3
1   2  x 
L
6 EI
1
w
P
L
L
L  2m, I  0.02m4 , w  2kg / m, P  20kg, E  2 106 kg / m2
のときの
1 , 2 , x
2
を求めよ。
（なお，できる限り，Ｅｘｃｅｌにおける式定義等を使い，手計算を行わないこと）

Download Report