情報科学 第７回　数値解析(2)

数値解析(2) トピック 連立1次方程式の解 消去法 反復法 乱数とモンテカルロ法 乱数とは？ 疑似乱数 モンテカルロ法

連立1次方程式の数値解法 小規模な連立1次方程式の解法 大規模な連立1次方程式の解法 消去法 反復法 Gauss消去法 Gauss-Jordan法 大規模な連立1次方程式の解法 反復法 Jacobi法

消去法による解法 小規模な連立1次方程式は消去法によって解くことができる 消去法とは以下の演算を何回か行うことによって未知数を消去し、方程式をより簡単な方程式に変形して解を求める方法（筆算と同様） 1つの方程式に(0でない)ある数を掛ける 1つの方程式にある数を掛けて他の方程式に加える 多くの消去法では、未知数を1回に1個ずつ消去することによって係数行列を三角行列に変形し、後退置換によって未知数の値を逐次求める

Gauss消去法 与えられた連立1次方程式を行列演算で表現し、係数行列を定義する 連立1次方程式 係数行列 定数ベクトル 消去法の説明を簡略化するため 右のような表記を導入する 筆算で式1に3を掛けて式2から 減算する操作を　　　　　　で表す

Gauss消去法（続き） 係数行列を三角行列に変形できるように各行に演算操作を行う 変形前 例えば で示された係数を0にするため 行1に3を掛けて行2から減算する 同様に で示された係数を消去 するために という操作を行う 変形後 この操作を左のように表記する 一番右の列は、変形における 行の演算操作を表す

はその都度あらわす値が変わっていることに注意 Gauss消去法（続き） 今度は で示された係数を調整する 対角成分は1に、非対角成分は0 になるように行の演算操作を行う 　　　　　　　　はその都度あらわす値が変わっていることに注意 最終的に係数行列が上三角行列 に変形されたら、後退置換により、 各変数の値が求まる

Gauss-Jordan法 Gauss消去法において、係数行列を単位行列に変形する方法 後退置換が不必要になる 例)　前出の連立方程式の解法

Pivoting 消去法では、係数行列の要素での割算における問題がある 割算の分母が0になる場合（係数行列の対角要素が0） 行または変数を入れ替える必要がある 割算の分母が0に近い値になる場合 割算の結果非常に大きな数字が結果として表れると数値計算の精度が失われる（前回のスライドを参照） この場合も行または変数を入れ替える必要がある この問題を解決するため以下のPivotingを行う 番目の行の処理で　　が0に近い場合、　　　 の絶対値 が大きくなり、　　行に 　　をかけて　　行から引いた 時に情報落ち誤差が生じるので良くない そこで、　　の中で絶対値が最大の　　　を選び　　　　　　　　 行と　　行を入れ替える 元の連立方程式では式の順序を入れ替えることに相当する　　 　行 　行

反復法 Gauss消去法は、単純明快という利点があるが、未知数の個数が多い場合演算回数が非常に多くなるという欠点がある 個の連立方程式から、１つの未知数を消去するために に比例する演算が必要となる これを 個の未知数に対して繰り返すためには、合計 　　　に比例する演算が必要となる また、物理現象の解析などの応用においては、疎な(sparse)係数行列を扱う場合や　　が大きい場合が多く　　　の演算回数は膨大になることがある 多くの未知変数を含む連立方程式の場合は、消去法の他に反復法が用いられることもある Jacobi法などは、行列のかけ算（　　の演算）を繰り返し、解を近似する 　回のかけ算で解が近似できる場合、計算量は となり、　　　　 の場合には反復法の方が有利となる

Jacobi法 Jacobi法は、 適当な初期値ベクトルに、行列の演算を繰り返し、解ベクトルに収束させる方法 連立方程式を次の形で表現する 係数行列 未知数ベクトル 定数ベクトル

Jacobi法（続き） 係数行列 を対角行列 と非対角行列 に分解する を定義する さらに （注意） ここで対角行列 の逆行列 係数行列 を対角行列 と非対角行列 に分解する さらに を定義する （注意） ここで対角行列 の逆行列 も対角行列であり簡単に求まる 一般に逆行列の計算は連立方程式の解法と同様の計算が必要 ※行列Aのある対角成分が0の場合は、 行を入れ替える（式の順番を入れ替える）か 列を入れ替える（変数を入れ替える）などの 操作を行い対角成分が0でないようにし、 　行列Dの逆行列 D-1が存在するようにする

Jacobi法（続き） を用いてベクトル列 を右の漸化式で定義する このとき の極限で は方程式 の解ベクトル に収束する   (証明のスケッチ) の両辺に をかけて  ：単位行列  ( と の定義 )

ノイマン級数 Jacobi法では、 という方程式を解く際に の逆行列として、無限級数 を考えた この無限級数 をノイマン級数と呼ぶ ノイマン級数　　　が収束するための十分条件は 行列 の最大ノルム が1より小さくなること (最大ノルム：行の要素の絶対値の総和をとりその中で最大となるもの)

Jacobi法の例 注意： Bの最大ノルム＜１

Jacobi法の例（続き） 一方、消去法による厳密解は、 となる は の周りを振動しながら に収束する

Jacobi法の計算例 スライド15と同じ例を用いる (1)行列A(ma)とベクトルb(b)の初期値設定 (2) 各変数の初期化 Ruby スライド15と同じ例を用いる (1)行列A(ma)とベクトルb(b)の初期値設定 (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] 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] (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])

Jacobi法の計算例（続き） (5)プログラムの続き:逐次的に解を計算する これを実行し計算ステップの数と誤差をプロットすると図の様になる Ruby (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] 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" } これを実行し計算ステップの数と誤差をプロットすると図の様になる 誤差は指数関数的に減少していることがわかる

プログラム別解 def jacob(a,b) Ruby 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 nx[i] = (b[i]-r)/a[i][i] d = dist(x,nx) print "nx=", nx, " d=", d, "\n" x[i] = nx[i] x プログラム別解 def dist(x,y) d = 0 for i in 0..x.size-1 d += (x[i]-y[i])**2 end Math.sqrt(d) jacob([[5.0,0,1.0], [3.0,4.0,0.0], [1.0,1.0,3.0]], [3,2,1])

2．乱数とモンテカルロ法 乱数とは？　 一様乱数 正規乱数 疑似乱数 モンテカルロ法

乱数とは？ 乱数とは、ある一つの数が現れたとき、それに続く数が前の数と全く関係なく現れるような列（乱数列）の要素 乱数列とは、ある分布に従う（互いに独立な事象を表す）確率変数の実現値の列をいう 乱数とは、乱数列の各要素である もう少し丁寧に定義すると… ある確率変数 の実現値の列である乱数列 において が の値をとることは からは予測できない 確率分布によって、一様乱数、正規乱数などがある。

一様乱数 出現確率が値によらず一定である乱数を一様乱数という（確率分布が一様である場合） 例えば、 サイコロの振って出た数字を並べたものは、続けて出現 する二つの数の間には因果関係がないため乱数列である また、（サイコロに細工がしてなければ）1から6までの数字が 1/6の（一定の）確率で出現するため一様乱数である 年末ジャンボ宝くじの当選番号（数字部分）を毎年記録した ものは一様乱数だろうか？ 矢を射ることによって全く公平に行われるとすると これも一様乱数となる。

自然現象における乱数列 放射性原子は放射線を出して崩壊する （原子核が崩壊して別の核種になるか、励起状態の原子核 がより低いエネルギー状態へ遷移する） 放射性物質の単位時間あたりの崩壊数を計測する 個々の原子核が他と無関係に崩壊するならば乱数列 同一の核種は平均として同じ割合で崩壊するが、 　崩壊数は同じ確率で出現するわけではない 　（平均からのずれがある） このように、一般には、一様ではない乱数列を扱う必要もある

正規乱数 出現確率が正規分布に従うような乱数を正規乱数という 確率変数 が よりも小さい値をとる確率 が以下の正規分布に従う場合 確率変数　　が 　よりも小さい値をとる確率　　 　が以下の正規分布に従う場合 自然現象を扱うシミュレーション（計算機で模擬実験を行うこと） では、一様乱数と共に正規乱数をよく用いる

疑似乱数 (pseudorandom numbers) 乱数列のように見えるが、計算機である計算を行うことによって求めることができる数列をいう 区間[0,1)の一様乱数を近似する数列をさすことが多い 疑似乱数は、n個目までのm個の既知要素 を用いて以下の漸化式によりn+1番目の要素を計算するものが多い 疑似乱数列 は周期が大きいことが望ましい （実際には、一様分布に関する適合度、統計的独立性など を様々な乱数検定法で検定し疑似乱数としての資格を与える）

モンテカルロ法 数学的問題の解法に乱数（疑似乱数）を用いる方法の総称を（広義の）モンテカルロ法という 一般に以下の２種類の問題を扱う 決定論的問題 乱数を用いる多次元積分 例えば、最も簡単な例では円周率の近似計算などがある。 非決定論的（確率的変動を含む）問題 自然科学・社会科学におけるシミュレーション 　　（実際に実験が困難であったり多大な費用がかかる場合） 　　例えば、ランダムウォーク（ブラウン運動） モンテカルロとはモナコにある賭博で有名な都市 （賭博場が国営）であり、多くの賭博は乱数を利用 して行われるためこのように呼ばれる

精度を上げるためには、 非常に多くの 点を用いる必要がある 決定論的問題 円周率の近似計算 Ruby 平面上の正方領域　　　　　　　　 にランダムに　　個の点を配置し、その中で、 四半円領域 にある点の数　　を求める 1 一様乱数を用いる場合、配置される点の数は 領域の面積に比例するはずであるから、 つまり円周率は　　　　　　　　で近似できる n=1000000 m=0 n.times{ x=rand　#0以上１未満の一様疑似乱数を返す y=rand if (x*x+y*y)<1.0 m+=1 end } puts 4*Float(m)/n O 1 のうち 　を満たす点 　を満たす点 精度を上げるためには、 非常に多くの　　　　　　　　点を用いる必要がある

Random Walk 実行結果 Ruby m=30000 #30000個の粒子 n=100 #粒子の移動の回数 非決定論的問題 Random Walk Ruby 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"} この確率過程は正規分布に従う？

Random Walk(補足） Ruby 非決定論的問題 一次元のRandom Walkを解析的に解く １つ粒子のみについて考え、n回の移動のあとの粒子の位置をxとする n回の移動のうち正方向にn1回、 負方向にn2回移動すれば 粒子が位置xに到達する場合の数は 粒子がn回の移動のあと位置xにいる確率は、n+x, n-xが奇数のときは0, それ以外では n>>1におけるスターリングの公式 および、y<<1におけるTaylor展開 を用いる 前スライドの確率過程はx<<nの極限で正規分布に従うことがわかる

数値解析(2) まとめ 連立1次方程式の解 乱数とモンテカルロ法 消去法 反復法 乱数とは？ 疑似乱数 モンテカルロ法 ランダムウォーク 円周率の計算

ノイマン級数の収束（補足） 物理学で扱うラプラス方程式やポアソン方程式では、ノイマン級数にあらわれる行列が対角優越となる場合が多い 対角優越の行列では、行列要素の適当な操作によって最大ノルムを1より小さくできる