● 前回の講義のまとめ
★
LU
分解• LU
分解とは, 正方・正則な行列A
を,上三角行列U ,
下三角行列L
を用いて,A = LU
と分解する方法である.
•
一旦, LU 分解が与えられれば, 連立一次方程式Ax = b
は, 前進代入U x = L
−1b
と, 後退代入x = U
−1y
によって,O(n
2)
回の計算で解くことができる.•
実際のLU
分解の,L
の対角成分を1
に取るアルゴリズムは以下のようなものである.L =
α
11α
21α
22.. . . ..
α
n1· · · · · · α
nn
, U =
β
11β
12· · · β
1nβ
22· · · β
2n. .. .. . β
nn
とおき,
A = LU
を求める.– j = 1, . . . , n
に対して以下を繰り返す.1. i = 1, . . . , j
に対して,β
ij= a
ij− (α
1iβ
1j+ · · · + α
i−1,iβ
ij) 2. i = j + 1, . . . , n
に対してα
ij= 1 β
jj(a
ij− (α
1iβ
1j+ · · · + α
j−1,j−1β
j−1,j)
このアルゴリズムは,「枢軸選択」を行っていない.•
枢軸選択つきのLU
分解のアルゴリズムは以下の通りである.– j = 1, . . . , n
に対して以下を繰り返す.1. i = 1, . . . , j − 1
に対して,β
ij= a
ij− (α
1iβ
1j+ · · · + α
i−1,iβ
ij) 2.
β
jj= a
jj− (α
1iβ
1j+ · · · + α
j−1,iβ
jj)
を計算し,p
j= β
jj とおく.3. i = j + 1, . . . , n
に対してα
ij= (a
ij− (α
1iβ
1j+ · · · + α
j−1,j−1β
j−1,j)
とし,|p
j| < |α
ij|
ならば,i
列目とj
列目を交換し,p
j= α
ij とおく.4. i = j + 1, . . . , n
に対してα
ij をp
j で割る.• LU
分解は,最初に与えられたA
と同一のメモリエリアにLU =
β
11β
12· · · β
1nα
21β
22· · · β
2n.. . . .. ... .. . α
n1· · · · · · β
nn
を上書きすることが可能である.
•
三重対角行列A =
b
1c
10 · · · 0 a
2b
2c
2· · · 0
· · ·
· · · a
n−1b
n−1c
n−1· · · 0 a
nb
n
に対しては, LU分解に関して次の重要な性質がある.
–
三重対角行列のLU
分解,前進・後退代入はO(n)
で実行可能.–
三重対角行列が「対角優位」すなわち|b
i| > |a
i| + |c
i|, i = 1, · · · , n
が成り立つとき, LU分解において,枢軸選択の必要がない.この事実は,次のようにして証明できる. この行列に対する
LU
分解のアルゴリズムはβ
jj= b
j− α
jj−1β
j−1j,
β
j−1j= c
j−1, α
jj−1= a
jβ
j−1j−1(1)
となる. そこで,帰納的に
|β
k−1k−1| > |c
k−1|
を証明すればよい.★ 反復法
•
連立一次方程式Ax = b
を解くもうひとつの方法は,反復計算に基づく方法である.x ∈ R
n に対してx
(k+1)= Φ(x
(k))
により反復計算を行った場合,もし{x
(k)}
が収束するならばΦ
の固定点(x = Φ(x)
を満たす点)に収束する. また,初期条件x
(0) を固定点の十分近くにとった場合,固定点x
におけるΦ
の微分|Φ(x)|
が|Φ(x)| < 1
を満たすならば,その反復計算は収束する. （これは,以前に証明した事実である）
•
反復計算x
(k+1)= Φ(x
(k))
が,ある線形写像A
を用いて,x
(k+1)= Ax
(k)+ b
と書けるならば,|Φ(x)|
は
A
の「作用素ノルム」kAk = max
x∈Rn|Ax|
|x|
と等しい. したがって,A
の絶対値最大の固有値が1
未満であれば,反復計算x
(k+1)= Ax
(k)+ b
は,任意の初期条件に対して収束する. また,そのとき の収束は１次収束となり,|x − x
(k)| ≤ C|λ|
k が成り立つ. ここで,λ
はA
の絶対値最大の固有値である.★
Jacobi
の反復法•
以下では,A
の対角部分からなる行列をD,
上三角部分（対角成分を除く）をU ,
下三角部分（対角 成分を除く）をL
とおく. すなわち,A = D + L + U
である.• Jacobi
の反復法とは,Ax = b
を解くためにx
k= 1
a
kk(b
k− a
k1x
1+ · · · + a
knx
n) k = 1, . . . , n
と考える方法である.•
これは,A
J= D
−1(L + U )
とおき,x
(k+1)= A
Jx
(k)+ D
−1b
によって反復する方法である.★
Gauss-Seidel
の反復法• Jacobi
の反復法とは,Ax = b
を解くためにx
(i+1)k= 1
a
kkb
k− a
k1x
(i+1)1+ · · · + a
k,k−1x
(i+1)k+ a
k,k+1x
(i)k+ · · · + a
knx
(i)nk = 1, . . . , n
と考える方法であり, Jacobiの反復法を,すでに計算が終わった要素を用いて反復するように改良し た方法である.•
x
(k+1)= D
−1(Lx
(k+1)+ U x
(k)) + D
−1b
と書くことができる. すなわち,A
GS= (D
−1− L)
−1U
とおき,x
(k+1)= A
GSx
(k)+ (D − L)
−1b
と反復する方法である.★ 反復法の収束
•
反復法が収束するか否かは,具体的なA
に対してA
J またはA
GSの絶対値最大の固有値を計算する 必要がある.• (N − 1) × (N − 1)
行列A =
−2 1
1 −2 1
.. . . .. ...
1 −2
の場合には
λ(A
J)
k= − cos(kπ/N), λ(A
GS) = cos
2(kπ/N ), or 0
（ただし,
A
GS の0
でない固有値は,すべて相異なる）となる. したがって,この場合には, Jacobiの反復法も
Gauss-Seidel
の反復法もともに収束する. 収束の速度はGauss-Seidel
法の方が幾らか高速である.
•
なお,A
の固有値を見れば明らかなように,A
に対する反復法では, 行列サイズが大きくなるほど収 束の速度は遅くなる.● 実習内容
以下では,
A =
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
とおく.
1. A
の絶対値最大固有値とその固有ベクトルを求めなさい.2. A
の絶対値最小固有値とその固有ベクトルを求めなさい.● レポート問題
★ レポート問題
以下の問題は評価対象のレポート問題です.
【レポート問題
06】
常微分方程式の境界値問題x
′′(t) − x(t) = sin(t), t ∈ [0, 1], x(0) = x(1) = 0
の数値解を
[0, 1]
区間を100
等分して計算し,解のグラフを書きなさい.【レポート問題
07】
A =
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
の絶対値最大の固有値とその固有ベクトルを数値的に求めなさい. 行列のサイズは
20 × 20
とします.これらの問題については,
A: 10
点, C: 0 点,と採点します. （中間的な点をつける可能性あり）★ 注意事項
【締め切り】 締め切りは２００９年８月５日（水）とします.
その他の注意事項については,前回のレポートの注意事項と同じです.
★ レポート問題
(Extra)
以下の問題は評価対象のレポート問題ですが,この問題を「評価の点の満点」にはカウントしません.（つ まり,以下の問題での得点は「ボーナスポイント」となります.）
【レポート問題
E07】
レポート問題07
の行列A
に対して,その絶対値最小の固有値と固有ベクトルを 求めなさい.【レポート問題
E08】
レポート問題07
の行列A
に対して,そのすべての固有値と固有ベクトルを求め なさい.【レポート問題
E09】
A
(2)=
5 −4 1
−4 6 −4 1
1 −4 6 −4 1
. .. ... ... ... ...
1 −4 6 −4 1
1 −4 6 −4
1 −4 5
のすべての固有値と固有ベクトルを求めなさい.（注意：
A
(2)= A
2 であることを使うと,A
の固有 値からA
(2) の固有値を計算することができる. しかし,それを使うのではなく,A
(2) を３重対角行列 に変換して固有値を計算する.）なお,これらの問題についての行列サイズは
20 × 20
程度でかまわない. ただし,「すべての固有値と固有 ベクトル」を求めるときには,その結果が分かりやすくなるように何らかの工夫をすること.問題
