情報科学 - lecture.ecc.u

情報科学
第７回数値解析(2)
1
数値解析(2) トピック
• 連立1次方程式の解
– 消去法
– 反復法
• 乱数とモンテカルロ法
– 乱数とは？
– 疑似乱数
– モンテカルロ法
2
連立1次方程式の数値解法
• 小規模な連立1次方程式の解法
– 消去法
• Gauss消去法
• Gauss-Jordan法
• 大規模な連立1次方程式の解法
– 反復法
• Jacobi法
3
消去法による解法
• 小規模な連立1次方程式は消去法によって解くこ
とができる
• 消去法とは以下の演算を何回か行うことによって
未知数を消去し、方程式をより簡単な方程式に変
形して解を求める方法（筆算と同様）
– 1つの方程式に(0でない)ある数を掛ける
– 1つの方程式にある数を掛けて他の方程式に加える
• 多くの消去法では、未知数を1回に1個ずつ消去す
ることによって係数行列を三角行列に変形し、後
退置換によって未知数の値を逐次求める
4

Gauss消去法
• 与えられた連立1次方程式を行列演算で表現し、
係数行列を定義する
x yz2
3x  5y  7z  0
2x  3y  z  5
連立1次方程式

消去法の説明を簡略化するため
右のような表記を導入する
筆算で式1に3を掛けて式2から
減算する操作を r2  3r1 で表す
1 1 1x  2

   
3 5 7y  0

   
2 3 1 z  5
係数行列
1
1
1 2
定数ベクトル
r1
3 5 7 0 r2
2 3 1 5 r3
5
Gauss消去法（続き）
係数行列を三角行列に変形できるように各行に演算操作を行う
変形前
1
1
1 2
r1
3 5 7 0 r2
2 3 1 5 r3
例えば
行1に3を掛けて行2から減算する r2  3r1
同様に
1
0 2
0 5
1

2
r1
この操作を左のように表記する
r3  2r1
一番右の列は、変形における
行の演算操作を表す

4 6 r2  3r1
3
で示された係数を消去
するために r3  2r1 という操作を行う
変形後
1
で示された係数を0にするため
1
6
Gauss消去法（続き）
1
1
1
2
r1
0 2 4 6 r2
0 5 3 1 r3
1 1 1
2
今度は
対角成分は1に、非対角成分は0
になるように行の演算操作を行う
r1 , r2 , r3 はその都
r1
0 1 2 3
r2 /2
0 0 7 14 r3  5 /2r2
1 1 1
2
r1
0 1 2 3
r2
0 0 1 2 r3 / 7
で示された係数を調整する
度あらわす値が変わ
っていることに注意
最終的に係数行列が上三角行列
に変形されたら、後退置換により、
各変数の値が求まる
x  3 y 1 z  2
7
Gauss-Jordan法
• Gauss消去法において、係数行列を単位行列に
変形する方法
• 後退置換が不必要になる
例) 前出の連立方程式の解法
1
1
1 2
r1
3 5 7 0 r2
2 3 1 5 r3
1 0
1
5
r1  r2 /2
0 1 2 3
r2 /2


0 0 7 14 r3  5 /2r2
1
1
1
2
r1
0 2 4 6 r2  3r1
0 5 3 1 r3  2r1
1 0 0 3
r1  r2 / 7
0 1 0 1 r2  2 / 7r3
0 0 1 2
r3 /7
x  3 y 1 z  2
8
Pivoting
• 消去法では、係数行列の要素での割算における問題がある
– 割算の分母が0になる場合（係数行列の対角要素が0）
行または変数を入れ替える必要がある
– 割算の分母が0に近い値になる場合
割算の結果非常に大きな数字が結果として表れると数値計算の
精度が失われる（前回のスライドを参照）
この場合も行または変数を入れ替える必要がある
この問題を解決するため以下のPivotingを行う
1
k行
0
akk
0 ak 1,k
i行
0
ai,k
ak,k 1
ak 1,k 1

ai,k 1
k 番目の行の処理で akk が0に近い場合、 aik / akk の絶対値
が大きくなり、 k 行に aik / akk をかけて i 行から引いた
時に情報落ち誤差が生じるので良くない
そこで、aik の中で絶対値が最大の

 a pk を選び ( p  k )
k 行と p 行を入れ替える



9
元の連立方程式では式の順序を入れ替えることに相当する
反復法
• Gauss消去法は、単純明快という利点があるが、未知数の個
数が多い場合演算回数が非常に多くなるという欠点がある
– n 個の連立方程式から、１つの未知数を消去するために
n 2 に比例する演算が必要となる

– これをn 個の未知数に対して繰り返すためには、合計
n 3 に比例する演算が必要となる
 • また、物理現象の解析などの応用においては、疎な(sparse)

係数行列を扱う場合や n が大きい場合が多く n 3 の演算回
数は膨大になることがある
• 多くの未知変数を含む連立方程式の場合は、消去法の他に
反復法が用いられることもある


2
– Jacobi法などは、行列のかけ算（ n の演算）を繰り返し、解を近似する
2
– k 回のかけ算で解が近似できる場合、計算量は k  n となり、 k  n
10
の場合には反復法の方が有利となる

Jacobi法
• Jacobi法は、適当な初期値ベクトルに、行列の
演算を繰り返し、解ベクトルに収束させる方法
連立方程式を次の形で表現する
a11 a12

a21

A


an1
係数行列
a1n 
x1 

 
 x  x 2 

 


 
ann 
x n 
A x  b
b1 
 
b2 

b
 
 
bn 
未知数ベクトル定数ベクトル
11
Jacobi法（続き）
係数行列 A を対角行列 Dと非対角行列 C に分解する
A  D C
A  (aij ) D  (dij ) C  (c ij )
aij  dij  c ij dii  aii dij  0 c ij  aij (i  j)


 さらに
B  D1  C c D1  b を定義する



（注意）
1
ここで対角行列 D の逆行列 D も対角行列であり簡単に求まる
一般に逆行列の計算は連立方程式の解法と同様の計算が必要

d 1
0
0 
11
※行列Aのある対角成分が0の場合は、


 0 d 1
行を入れ替える（式の順番を入れ替える）か
1


22
D 
列を入れ替える（変数を入れ替える）などの


0
操作を行い対角成分が0でないようにし、
12


1
行列Dの逆行列 D-1が存在するようにする
0 dnn 
 0
Jacobi法（続き）
B  D  C c D  b を用いてベクトル列 x
1
1
を右の漸化式で定義する
x
(0)
0 x
(n 1)
(n )
 B x
(n )
c
n  の極限で x は方程式 A x  b

( )
の解ベクトル x に収束する
(n )
このとき

(証明のスケッチ)
A  x  (D  C)  x  
b の両辺に D 1 をかけて


1
D
(D  C)  x  D1b  (I  D1  C)  x  D1b
I ：単位行列
c
 I  B  B 2  ... B n  ...c
IB
 B  ... B  B
 0  c   c  c ...13  c
(I  B) x  c  x 
( B と c の定義)


ノイマン級数
• Jacobi法では、 (I  B) x  c という方程式を解く際
に (I  B) の逆行列として、無限級数 Q を考えた
c
x
 I  B  B2  ... Bn  ... c  Q c
IB



n
Q  I  B  B 2  ... B
 ...
この無限級数 Q をノイマン級数と呼ぶ
Q が収束するための十分条件は
ノイマン級数

n
行列 B の最大ノルム

B   max bij が1より小さくなること
i
j
(最大ノルム：行の要素の絶対値の総和をとりその中で最大となるもの)

14
Jacobi法の例
5x  z  3
3x  4 y  2
x  y  3z  1

5 0 1x  3

   
3
4
0

y  2

   
1
1
3

z  1
0 0 1
5 0 0
1/5 0
0 






C  3 0 0 D  0 4 0 D1   0 1/4 0 


 



1
1
0
0
0
3
0
0
1/3






1/5 0
0 0 0 1  0
0
1/5


 

B  D1  C   0 1/4 0 3 0 0 3/4
0
0 



 

0
0
1/3
1
1
0
1/3
1/3
0


 

1/5 0
0 3 3/5
注意：






c  D1  b   0 1/4 0 2 1/2 
Bの最大ノルム＜１

   
15
0
0
1/3

1 1/3
Jacobi法の例（続き）
 0
0
1/53/5 3/5  0.53333

    

(2)
(1)
x  B  x  c  3/4
0
0 1/2  1/2   0.05 

    

0 1/3 1/3 0.03333
1/3 1/3
 0
0
1/5 0.53333 3/5 0.606666


   

(3)
(2)
x  B  x  c  3/4
0
0  0.05  1/2   0.1 


   

0 0.03333 1/3 0.138888
1/3 1/3
0.580444
0.572222



 (5)
(4 )
(4 )
(3)
x  B  x  c   0.045  x  B  x  c  0.070833




0.127593
0.097777
0.5763


一方、消去法による厳密解は、 x*  0.0678 となる


0.1186

16


(n )
x は x* の周りを振動しながら x* に収束する
Ruby
Jacobi法の計算例
スライド15と同じ例を用いる
(1)行列A(ma)とベクトルb(b)の初期値設定
ma=[[5.0,0,1.0],[3.0,4.0,0.0],[1.0,1.0,3.0]]
b=[3,2,1]
(3) 行列D(md)，D-1(mdi)， C(mc) の初期化
for i in 0..2
for j in 0..2
if (i==j)
md[i][i]=ma[i][i]
mdi[i][i]=1.0/ma[i][i]
else
md[i][j]=mdi[i][j]=0
end
mc[i][j]=ma[i][j]-md[i][j]
end
end
(2) 各変数の初期化
mb=[[0,0,0],[0,0,0],[0,0,0]]
mc=[[0,0,0],[0,0,0],[0,0,0]]
md=[[0,0,0],[0,0,0],[0,0,0]]
mdi=[[0,0,0],[0,0,0],[0,0,0]]
c=[0,0,0]
(4) 行列B (mb)，ベクトルc(c) の初期化
for i in 0..2
for j in 0..2
for k in 0..2
mb[i][j]+=(-mdi[i][k]*mc[k][j])
end
c[i]+=(mdi[i][j]*b[j])
end
end
17
Ruby
Jacobi法の計算例（続き）
(5)プログラムの続き:逐次的に解を計算する
iter=15 #計算ステップの回数
xs=[0.5763,0.0678,0.1186] #解析解
x=[0,0,0]
iter.times{ |k|
nx=[0,0,0]
for i in 0..2
for j in 0..2
nx[i]+=(mb[i][j]*x[j])
end
nx[i]+=c[i]
end
x=nx
err=Math.sqrt((x[0]-xs[0])**2+(x[1]-xs[1])**2+(x[2]-xs[2])**2)
print "k= ",k,"\t","x=[ ",x[0],",",x[1],",",x[2],"] err= ", err, "\n"
}
#解析解との誤差
これを実行し計算ステップの数と誤差をプロットすると図の様になる
誤差は指数関数的に減少していることがわかる
18
Ruby
プログラム別解
def dist(x,y)
d = 0
for i in 0..x.size-1
d += (x[i]-y[i])**2
end
Math.sqrt(d)
end
jacob([[5.0,0,1.0],
[3.0,4.0,0.0],
[1.0,1.0,3.0]],
[3,2,1])
def jacob(a,b)
n = b.size
x = Array.new(n,0)
nx = Array.new(n,0)
d = 1.0
while d>0.000001
for i in 0..n-1
r = 0
for j in 0..n-1
if j!=i
r += a[i][j]*x[j]
end
end
nx[i] = (b[i]-r)/a[i][i]
end
d = dist(x,nx)
print "nx=", nx, " d=", d, "\n"
for i in 0..n-1
x[i] = nx[i]
end
end
x
19
end
2．乱数とモンテカルロ法
1.
乱数とは？
•
•
2.
3.
一様乱数
正規乱数
疑似乱数
モンテカルロ法
20
乱数とは？
• 乱数とは、ある一つの数が現れたとき、それに続く数が前の数と
全く関係なく現れるような列（乱数列）の要素
もう少し丁寧に定義すると…
• 乱数列とは、ある分布に従う（互いに独立な事象を表す）確率変
数の実現値の列をいう
• 乱数とは、乱数列の各要素である
ある確率変数 x の実現値の列である乱数列
{x1, x 2,...,x n , x n 1,...} において
x が x n 1 の値をとることは x1, x 2 ,...,x n からは予測できない



確率分布によって、一様乱数、正規乱数などがある。
21
一様乱数
• 出現確率が値によらず一定である乱数を一様乱数
という（確率分布が一様である場合）
例えば、
サイコロの振って出た数字を並べたものは、続けて出現
する二つの数の間には因果関係がないため乱数列である
また、（サイコロに細工がしてなければ）1から6までの数字が
1/6の（一定の）確率で出現するため一様乱数である
年末ジャンボ宝くじの当選番号（数字部分）を毎年記録した
ものは一様乱数だろうか？
矢を射ることによって全く公平に行われるとすると
22
これも一様乱数となる。
自然現象における乱数列
放射性原子は放射線を出して崩壊する
（原子核が崩壊して別の核種になるか、励起状態の原子核
がより低いエネルギー状態へ遷移する）
放射性物質の単位時間あたりの崩壊数を計測する
• 個々の原子核が他と無関係に崩壊するならば乱数列
• 同一の核種は平均として同じ割合で崩壊するが、
崩壊数は同じ確率で出現するわけではない
（平均からのずれがある）
このように、一般には、一様ではない乱数列を扱う必要もある
23
正規乱数
• 出現確率が正規分布に従うような乱数を正規乱数
という
• 確率変数 x が  よりも小さい値をとる確率 P(x)
が以下の正規分布に従う場合

1
 P(x   )  2
 x 2 
 exp 2 dx


自然現象を扱うシミュレーション（計算機で模擬実験を行うこと）
では、一様乱数と共に正規乱数をよく用いる

24
疑似乱数
(pseudorandom numbers)
• 乱数列のように見えるが、計算機である計算を
行うことによって求めることができる数列をいう
• 区間[0,1)の一様乱数を近似する数列をさすこと
が多い
疑似乱数は、n個目までのm個の既知要素 x n , x n1,...,x nm 1
を用いて以下の漸化式によりn+1番目の要素を計算するものが多い
x n 1  f (x n , x n1,...,x nm 1)

疑似乱数列 {x n } は周期が大きいことが望ましい

（実際には、一様分布に関する適合度、統計的独立性など
を様々な乱数検定法で検定し疑似乱数としての資格を与える）
25
モンテカルロ法
• 数学的問題の解法に乱数（疑似乱数）を用いる方
法の総称を（広義の）モンテカルロ法という
モンテカルロとはモナコにある賭博で有名な都市
（賭博場が国営）であり、多くの賭博は乱数を利用
して行われるためこのように呼ばれる
• 一般に以下の２種類の問題を扱う
– 決定論的問題
• 乱数を用いる多次元積分
例えば、最も簡単な例では円周率の近似計算などがある。
– 非決定論的（確率的変動を含む）問題
• 自然科学・社会科学におけるシミュレーション
（実際に実験が困難であったり多大な費用がかかる場合）
例えば、ランダムウォーク（ブラウン運動）
26
決定論的問題
円周率の近似計算
Ruby
平面上の正方領域 P  {0  x  1, 0  y  1}
y
にランダムに n 個の点を配置し、その中で、
2
2
四半円領域 Q  {x  y 1, 0  x 1, 0  y 1}
1
にある点の数 m を求める
一様乱数を用いる場合、配置される点の数は


領域の面積に比例するはずであるから、
P  {0  x  1, 0  y  1}


m
/ n  Q / P   /4
つまり円周率は
  4m/n で近似できる
n=1000000
O
1
m=0

n.times{
Q  {x 2  y 2 1, 0  x 1, 0  y 1}
x=rand #0以上１未満の一様疑似乱数を返す
y=rand
x i  [0,1) y i  [0,1) i 1...n のうち
if (x*x+y*y)<1.0
x 2  y 2 1 を満たす点
i
i
m+=1
2
2
x i  y i 1 を満たす点
end
}



27
puts 4*Float(m)/n
x


精度を上げるためには、
非常に多くの
(n 1)
点を用いる必要がある
非決定論的問題
Ruby
Random Walk
1次元格子上をm個の粒子が原点を起点とし、互いに独立に確率1/2で+1あるいは-1移動する。
各粒子がn回移動した後、格子点上のどの点にどれだけいるか統計をとり、
粒子の位置の確率分布をしらべよう
実行結果
m=30000 #30000個の粒子
n=100
#粒子の移動の回数
hist=Array.new(n*2,0)
m.times{
x=0
n.times{ |i|
p=rand #一様乱数を確率過程に用いる
if (p>0.5)
x+=1
# 確率1/2で+1移動
else
x-=1
# 確率1/2で-1移動
end
}
hist[n+x]+=1
}
hist.size.times{ |i| print i-n,"\t",hist[i],"\n"}
28
この確率過程は正規分布に従う？
非決定論的問題
Ruby
Random Walk(補足）
一次元のRandom Walkを解析的に解く
１つ粒子のみについて考え、n回の移動のあとの粒子の位置をxとする
n回の移動のうち正方向にn1回、負方向にn2回移動すれば
x  n1  n 2
n  n1  n 2
(n1  n 2 )!
n1!n 2
!
粒子がn回の移動のあと位置xにいる確率は、n+x, n-xが奇数のときは0,
それ以外では

1
1
1
n!
P(n, x)  ( ) n1 ( ) n 2 W (n1,n2 )  ( ) n
2
2
2 ( n  x )!( n  x )!
2
2
n x
nx
ln P(n, x)  n ln2  ln n!ln(
)!ln(
)!
2
2
n>>1におけるスターリングの公式
2
n x
x
1
1
1
 ln
 { ( ) 2  o( ) 3 }  o(1/n)

lnn! (n  )lnn  n  ln2  o( )
2
2
n
n
2 n
n
粒子が位置xに到達する場合の数は W (n1,n 2 ) 

P(n, x) 
および、y<<1におけるTaylor展開
2
2
x
exp( )
n
2n

1
ln(1 y)  y  y 2  o(y 3 ) を用いる
2
29
前スライドの確率過程はx<<nの極限で正規分布に従うことがわかる
数値解析(2) まとめ
• 連立1次方程式の解
– 消去法
– 反復法
• 乱数とモンテカルロ法
– 乱数とは？
– 疑似乱数
– モンテカルロ法
• ランダムウォーク
• 円周率の計算
30

Download Report