+ n(n-1) - Keio.ac.jp

数値・記号処理（４）
方程式を解く（一変数）
慶應義塾大学理工学部
櫻井彰人
線型方程式を解く
線型方程式を解くことは簡単である。
中学の数学に出てくる
 小学生高学年だって理解できる
 Cramer の公式で、一瞬！

本当に簡単か？

例えば、100変数ならどうか。 100 x 100 の行列式？
一般的に解く。解けることは分かる

係数を記号にして、四則のみで、解が記述可能
Mathematica で解いてみよう
http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear01.nb
クラメル(Cramer)の公式
連立一次方程式の基本的な解法であった
 2 2  1  x   5 
 1 1 3   y    7 

   
 5 1 2   z  13
x
5
2 1
2
5
1
2
2
5
7
1
1
7
3
1 1
7
3
13 1 2
2 2 1
y
5 13 2
2 2 1
z
5
2
1 13
2 1
1 1
3
1 1
3
1 1
3
5
2
5
2
5
2
1
1
1
Cramer の公式（一般の場合）
Cramer の公式 (Cramer’s rule) を使うと
a11 x1  a2 x2   an xn  d1
a21 x1  a22 x2   a2 n xn  d 2


  
am1 x1  am 2 x2   amn xn  d m
a11
a12
 a1n
a21 a22  a2 n
A


am1 am 2
am
xi 
A1 
d1
d2
a12  a1n
a22  a2 n

dn

an 2
ann
A2 
a11 d1  a1n
a21 d 2  a2 n

an1

dn
ann
Ai
A
http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear02-Cramer.nb
Cramer の公式の問題点
行列式の計算が大変

n行 n列の行列式の計算に (n-1) n! 回の乗算が必要
丸め誤差の影響がありうる

n! 個の項の加減
http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear03-det.nb
数値解析で用いる連立１次方程式の特徴
未知数が非常に多い
対称行列になることが多い
帯行列になることが多い
*
*

*

0
0

0
0

0
0

0
*
*
*
*
0
0
0
0
0
0
*
*
*
*
*
0
0
0
0
0
0
*
*
*
*
*
0
0
0
0
0
0
*
*
*
*
*
0
0
0
0
0
0
*
*
*
*
*
0
0
0
0
0
0
*
*
*
*
*
0
0
0
0
0
0
*
*
*
*
*
0
0
0
0
0
0
*
*
*
*
0
0

0

0
0

0
0

*
*

*
解法の分類
直接法

何らかの形で直接的に解く。理論的には有限回演算。丸め誤差
の可能性
Cramer 公式
消去法
直交化法
せいぜい６元
Gauss法
一般に使える
三角化
Crout法
Choleski法正定値対称の場合は標準的
掃き出し法(Gauss-Jorda法) 分かりやすい方法
安定性良. 複雑
反復法

トライアル・アンド・エラ－で解く。収束すれば丸め誤差の影響な
し。計算の打ち切りが必要
Gauss-Seidel法
線形反復法逐次過剰緩和法
Chebyshev反復
最急降下
勾配法
共役勾配
収束遅い
普通
速いが複雑
ナイーブ故
有限回
Gaussの消去法の例
Step 1

第一行を何倍かし、それを第二行以下に加減算し、
第二行以下の第一列成分を0とする
1
1
0   u1   1 4 
4
 1
2
0
 1 2 u2   1 8 

   

 1
0
2
 1 2 u3   1 8 

  
0

1
2

1
2
1

 u4  1 12 
1
0   u1   1 4 
4  1
0 7 4  1 4  1 2 u  3 16 


  2   

0  1 4 7 4  1 2 u3  3 16 

  
1  u4  1 12 
0  1 2  1 2
Gaussの消去法の例
Step 2

第二行を何倍かし、それを第三行以下に加減算し、
第三行以下の第一、二列成分を0とする
1
0   u1   1 4 
4  1
0 7 4  1 4  1 2 u  3 16 


  2   

0  1 4 7 4  1 2 u3  3 16 

  
0

1
2

1
2
1

 u4  1 12 
1
0   u1   1 4 
4  1
0 7 4  1 4  1 2  u   3 16 


  2   

0 0 12 7  4 7 u3   3 14 

  
0
0

4
7
6
7

 u4  23 168 
Gaussの消去法の例
Step 3

第三行を何倍かし、それを第四行以下に加減算し、第四
行以下の第一、二、三列成分を0とする
1
0   u1   1 4 
4  1
0 7 4  1 4  1 2  u   3 16 


  2   

0 0 12 7  4 7 u3   3 14 

  
0
0

4
7
6
7

 u4  23 168 
0   u1   1 4 
4  1  1
0 7 4  1 4  1 2  u   3 16 


  2   

0 0 12 7  4 7 u3   3 14 

  
0
0
0
2
3

 u4  5 24 
Gaussの消去法の例
Step 4

u4,u3,u2,u1の順に解を求めることができる
0   u1   1 4 
4  1  1
0 7 4  1 4  1 2  u   3 16 


  2   

0 0 12 7  4 7 u3   3 14 

  
0
0
0
2
3

 u4  5 24 
5 3 5
u4 
 
24 2 16
3 14   4 7 u4  3  4  5  7 11
u3 
      
12 7
14  7  16  12 48
 3  1  11  1  5  4 11
u2           
16  4  48  2  16  7 48
11
11  1 17
1
u1     1   1  
48
48  4 96
4
Gaussの消去法のプログラミング
通常のように
 a11 a12
a
a22
21

 a31 a32




am1 am 2
a13
a23
a33

am 3
 a1m   u1   f1 
 a 2 m   u2   f 2 
    
 a3 m   u 3    f 3 
   
   

   
 amm  um   f m 
前進消去
左下部分を0に
i  1　 m  1　 a11 a12
a
a22
 21
 a31 a32


 
am1 am 2
j  i  1　 m　 a ji 
f j  f j    f i
 aii 
k  i  1　 m　 a ji 
a jk  a jk   aik
 aii 
a13
a23
a33

am 3
 a1m   u1   f1 
 a 2 m   u2   f 2 
    
 a3 m   u 3    f 3 

       
   
 amm  um   f m 
後退代入
解を、下から、順番に
a11 a12
0 a
22

0
0



 0
0
fm
um 
amm
i  m  1　 1　( i が減少する順)
j  i  1　 m　f i  f i  aij f j
fi
ui 
aii
a13  a1m   u1   f1 
a23  a2 m   u2   f 2 
   
a33  a3m  u 3    f 3 
   
    

   
0  amm  um   f m 
プログラム例
http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear04-GE.nb
ガウスの消去法の計算量（１）
演算１ mi=aik/akk i=k+1, k+2, ..., n
演算２ aij = aij - mi×akj
i=k+1,
k+2, ..., n かつ j=k+1, k+2, ..., n
演算３ bi = bi - mi×bk
i=k+1, k+2, ...,
n
n
n
n
A
x = b
0
n
0
k
akk
n-k
ガウスの消去法の計算量（２）
前進部分での実行回数

乗除算


n(n-1) + n(n-1)(2n-1)/6
加減算

n(n-1)/2 + n(n-1)(2n-1)/6
後退部分での実行回数

乗除算


n(n-1)/2 + n
加減算

n(n-1)/2
行列式を求める
先ほど、行列式を求めるのは、計算量が大変といった。
しかし、ガウスの消去法を用いれば簡単。なぜなら
 a12

a11


a22

A  



  a1n 
a13
  a2 n 
a23

  a3 n 
a33

 
 
ann
 a22
 a33
  ann

det  A  det  A  a11
しかしなぜ？
行列式を求める
答え:
行列 A のある行の何倍かを他の行に加えて得ら
れる行列を A´ とすると
det(A) = det(A´)
それ（「行の何倍か、、、、」）はまた何故？
答え:
行列式の多重線形性
det(a1, a2,.., ai+ai´,.., an) =
det(a1, a2,.., ai,.., an) + det(a1, a2,.., ai´,.., an)
による
ガウス・ジョルダンの消去法
ガウスの消去法との違い：
ガウスの消去では


対角成分はもとのまま
（対角成分を除く）下三角成分のみを0とする
ガウス・ジョルダン消去法では

対角成分を1にし、
対角成分以外を0にする

後退代入は不要となる！

http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear05-GJ.nb
逆行列を求める
ガウス・ジョルダン消去法を繰り返せばよい
（LU分解と呼ばれる方法の方が効率がよい）
実際には、一気に計算する
 a11 a12

 a21 a22
 


a
 n1 an 2
 a11 a12
a
 21 a22
 


an1 an 2
 a1n  c11 c12

 a2 n  c21 c22
   


 ann  cn1 cn 2
 a1n  c11   1 
   

 a2 n  c21   0 

       
   
 ann  cn1   0 
 c1n   1 0
 
 c2 n   0 1


 
 
 
 cnn   0 0
 a11 a12
a
 21 a22
 


an1 an 2




0

0


1 
 0
 a1n  c1i   
  
 a2 n  c2i   
 1
      
    
 ann  cni   
 0
逆行列を計算する
実際には一気に計算する
 a11 a12

 a21 a22
 


a
 n1 an 2
1

0


0

0
1

0
 a1n  c11 c12

 a2 n  c21 c22
   


 ann  cn1 cn 2




0  c11 c12

0  c21 c22

  

1  cn1 cn 2
 c1n   1 0
 
 c2 n   0 1


 
 
 
 cnn   0 0
 a12

 c1n   a11
 
 a22

 c2 n   a21

    

 
 cnn   an1 an 2




0

0


1 
 a1n 


 a2 n 
  

 
 ann
http://www.sakurai.comp.ae.keio.ac.jp/classes/numsymbol-class/2004/Linear05-GJ.nb
ピボット選択
ピボット: ガウスの消去法で「対角要素」が 0 で
あったらどうしよう？
答え：別の行を持ってくる（絶対値最大が安全）
悪条件 ill-conditioned の方程式
2x1 – x2 = 3
2.1x1 – x2 = 3

Download Report