● 講義資料

● 講義資料

▼ 講義予定

• 反復法による連立一次方程式の数値解法

● 前回の講義のまとめ

▼ 消去法

★ はき出し法

• 以下,A∈Mn(R),b,x∈Rⁿ として,連立一次方程式 Ax=b

を数値的に解くことを考える. （A が正則でない場合には, 「正則でない」ことが判定でき ればよいことにする.）

• 連立一次方程式を数値的に解く方法としては,大きく分けて,消去法と反復法がある. はじめ に,消去法の中でももっとも基本的な「掃出し法」（「Gauss-Jordanの消去法」）を考える.

• 掃出し法とは, 「行基本変形」を繰り返して, 行列 Aを単位行列に変形する方法である. こ こで,「行基本変形」は,以下の３種類の変形である.

1. i行目をλ倍する.

これは,方程式の両辺に左からDi(λ)をかけることに他ならない.

2. i行目とj 行目を入れ替える.

これは,方程式の両辺に左からTij をかけることに他ならない.

3. i行目にj 行目のλ倍を加える. （ただし,j < iを仮定する）

これは,方程式の両辺に左からEij(λ)をかけることに他ならない.

ただし,行列Di(λ),Tij,Eij(λ)は,

Di(λ) =

i

<

































1 0 . . . 0

0 . .. ...

... . .. ...

0 . . . λ . . . 0

... . .. ...

... . .. 0

0 . . . 0 1

































< i

Tij=

i

<

j

<































1 . . . 0 ... . .. . . ... 0 · · · 0 · · · 1 · · · 0

... . .. ...

0 · · · 1 · · · 0 · · · 0 ... . . .. ...

0 . . . 1































< i

< j

Eij(λ) =

j

<

i

<































1 . . . 0 ... . .. . . ... 0 · · · 1 · · · 0 · · · 0

... . .. ...

0 · · · λ · · · 1 · · · 0 ... . . .. ...

0 . . . 1































< j

< i

である.

• 掃出し法のアルゴリズムは以下の通りである. ただし,「枢軸選択」を行っていないので,正 則であっても,アルゴリズムが最後まで到達できないときがある.

– i= 1, . . . , nに対して,以下を繰り返す.

1. 対角成分aii が 1となるように,i 行目を定数倍する.

すなわち,「Di(a⁻¹_ii )を左からかける」または,「aij :=aij/aii (j = 1, . . . , n) と する」

2. j= 1, . . . , i−1, i+ 1, . . . , nに対して,非対角成分aji を0となるように,j 列目か らi列目の定数倍を引く.

すなわち,「Eji(−aji)を左からかける」または,「ajk:=ajk−ajiaik (k= 1, . . . , n) とする」

– もし, 1のステップでaii= 0 であるときには,プログラムを停止する.

★ Gaussの消去法

• 連立一次方程式を与える行列Aが上三角行列U である場合には, U y=b

は以下のようにしてO(n²)の計算量で解くことができる. いま,U y=bをexplicitに書くと, a11y1+· · ·+a1nyn=b1,

...

an−1n−1yn−1+an−1nyn=bn−1, annyn=bn

であるので,すべてのiに対して,対角成分aii がnon-zeroと仮定すると（この仮定はU が 正則であることと同値である）

yn = 1 ann

bn, yn−1= 1

an−1n−1

(bn−1−an−1nyn), ...

y1= 1 a11

(b1−a1nyn− · · · −a12y2)

として,yn, yn−1, . . . , y1 の順序で陽的に解くことができる. これを後退代入と呼ぶ.

• したがって,連立一次方程式Ax=bを解くためには, Aに行基本変形を施して U x=b^′ の 形にすればよい. 実際,これは以下のように実行できる. これをGauss の消去法と呼ぶ.

i= 1, . . . , nに対して,以下を繰り返す.

1. aji (j≥i)の中で絶対値最大の要素を探し, それを含む行とi行目を交換し,その後aii

が1となるように i行目を定数倍する.

2. j=i+ 1, . . . , nに対して,非対角成分ajiを0 となるように,j 列目からi列目の定数 倍を引く.

このアルゴリズムは,Aが正則ならば必ず終了する. 実際,これが終了しない可能性があるの は, ステップ１において, {aji}ⁿ_j=i がすべて0となっている場合である. この場合には,i 列 目はi−1列以前の列の一次結合で書けることとなり, Aが正則であるという仮定と矛盾す る. よって, Gaussの消去法はA が正則であれば,必ず終了する.

• Gaussの消去法の計算量はO(n³)であるが,ステップ２のループ回数が掃出し法と比較して

半分となっているので,掃出し法のほぼ半分の計算量となる.

★ LU 分解

• 現実に連立一次方程式を解く際に,同一のA, 異なるb に対して,Ax=b を何度も解く場合 がある. 特に,あるb1に対してAx=b1 を解き,そのxから得られるb2 に対して解を求め る場合もある. このような場合, Gaussの消去法を何度も適用することは, O(n³)の計算を何 度も行うことになり,非常に効率が悪い.

もし,Aに対して,ある下三角行列L, 上三角行列U を用いて,事前にA=LU と分解され ていれば,Ax=b を解く手順は,LU x=b より,

として,ともにO(n²)のアルゴリズムで解を得ることができる. ここで,前進代入とは,後退 代入と同じ手順を下三角行列に適用したものである.

• Aを正則行列としたとき,次の定理が成り立つ.

任意の正則行列Aに対して,Aの行を適当に入れ替えれば,入れ替えたあとの行列を再びA と書いたとき,ある上三角行列U,下三角行列Lが存在して,A=LU と書ける. これをLU 分解と呼ぶ.

この証明はGauss の消去法のアルゴリズムから導かれる. ここでは, Gauss の消去法で枢軸 選択が必要ない（対角成分が 0にはならない）場合のみを考える. この時, ある基本変形行 列の列{Mi}^N_i=1と上三角行列 U が存在して,

MN· · ·M1A=U

と書ける. いま,Mi には行交換の行列が含まれないので,すべてのMi は下三角行列である.

よって,

A= (M₁⁻¹· · ·M_N⁻¹)U

となる. ここで, 下三角行列の逆行列は再び下三角となり,下三角行列の積は下三角なので, M₁⁻¹· · ·M_N⁻¹は下三角行列となる.

• 実際のLU分解は,L= (ℓij),U = (uij)と書いたとき,

aij=

j

X

k=1

ℓikukj, i > j,

aii=

i

X

k=1

ℓikuki=

i

X

k=1

ℓijukj, i=j,

aij=

i

X

k=1

ℓikukj, i < j,

(1)

となるが,n²+n個の未知数{ℓij, uij}に対して, n本の方程式があり,このままではL,U を定めることができない. そこで,

– Lの対角成分を1 にとる(Doolittle Type) – U の対角成分を 1にとる(Crout Type) のいずれか一方の条件をおくことが多い.

• いま, Crout Type (uii= 0)を仮定すると,

aij=

j−1

X

k=1

ℓikukj+ℓij, i≥j,

aij=

i−1

X

k=1

ℓikukj+ℓiiuij, i < j,

となるので,これを書き直すと, ℓij =aij−

j−1

X

k=1

ℓikukj, i≥j,

uij = 1

ℓii aij−

i−1

X

k=1

ℓikukj

!

, i < j,

(2)

となる.

• もっとも標準的なアルゴリズムは,Crout Algorithmと呼ばれる次の方法である. （Crout typeを仮定している.）

i= 1, . . . , nに対して,以下を行う.

1. j= 1, . . . , iに対して,

ℓij=aij−

j−1

X

k=1

ℓikukj

を計算する.

2. j=i+ 1, . . . , nに対して,

uij = 1 ℓii

aij−

i−1

X

k=1

ℓikukj

!

を計算する.

• この方法では（多くのLU分解のアルゴリズムでは）, L,U のデータはAの上に上書きで きる. 逆に,Aの上にL,U を上書きした方が,プログラムが書きやすい.

▼ 反復法

★ Jacobiの反復法

• Aは,n×n正方行列で,その対角成分は0 を含まないと仮定する. この時,連立一次方程式 Ax=b を解くために,以下の反復解法を考えることができる.

x⁽⁰⁾ を適当に与え,

x^(k+1)₁ = 1 a11

b1−a12x^(k)₂ − · · · −a1nx^(k)_n , ...

x^(k+1)_n = 1 ann

bn−an2x^(k)₂ − · · · −ann−1x^(k)_n−1 , これをJacobiの反復法と呼ぶ.

• Aを対角部分D,（対角部分を含まない）上三角部分U,（対角部分を含まない）下三角部 分Lを使ってA=D+L+U と書いたとき, Jacobiの反復法は

x^(k+1)=AJx^(k)+D⁻¹b, AJ =−D⁻¹(L+U) とかきあらわすことができる.

★ Gauss-Seidel法

• Jacobiの反復法で,x1 から順に計算しているとき,すでに計算済みのx^(k+1) の成分を利用

して反復計算を行った方が,収束が速いと予想できる. そのように改良した以下の反復法を Gauss-Seidel 法と呼ぶ.

x^(k+1)₁ = 1 a11

b1−a12x^(k)₂ − · · · −a1nx^(k)_n , ...

x^(k+1)_ℓ = 1 aℓℓ

bℓ−aℓ1x^(k+1)_ℓ − · · · −aℓℓ−1x^(k+1)_ℓ−1 −aℓℓ+1x^(k)_ℓ+1− · · · −ann−1x^(k)_n−1 , ...

x^(k+1)_n = 1 ann

bn−an2x^(k+1)₂ − · · · −ann−1x^(k+1)_n−1 ,

• Gauss-Seidel法は,

(D+L)x^(k+1)=−U x^(k)+b と書けているので,

x^(k+1)=AGSx^(k)+ (D+L)⁻¹b, AGS=−(D+L)⁻¹U

と書きあらわすことができる.

★ 反復法の収束

• 一般に反復計算は,任意の初期値x⁽⁰⁾ に対して収束するわけではないが,B を n×n行列と し,その絶対値最大の固有値ρ(B)がρ(B)<1をみたすならば,反復計算

x^(k+1)=Bx^(k)+c は,任意の初期値に対して収束し,収束値x^∞ は

x^∞=Bx^∞+c

を満たすことが証明できる. そのときの収束の速さはO(ρ(B)^k)となる.

• 具体的なAに関して,AJ,AGS がρ(AJ)<1,ρ(AGS)<1をみたすことを確かめるのは,一 般には容易ではないが,以下の定理が知られている.

– Aが正定値実対称行列ならば,ρ(AGS)<1が成り立つ.

– Aが対角優位,すなわち, すべての行に対して

|aii| ≤X

i6=j

|aij|

が成り立つならばρ(AJ)<1

• １次元の２階微分作用素 d²

dx² のDirichlet境界値問題に対応する行列

A=



















−2 1

1 −2 1

. .. ... ...

1 −2 1

1 −2



















に対して,以下が成り立つ. （行列サイズはN×N とする）

– AJ の固有値は

cos kπ N+ 1

N k=1

となる. よって,ρ(AJ)<1 が成り立つ.

– AGS の固有値は,⌈N/2⌉個が0であり,残りの⌊N/2⌋個は,

cos² kπ N+ 1

⌊N/2⌋

k=1

とな る. よって,ρ(AGS)<1が成り立つ.