2013年度・後期・数理解析・計算機数学３・第１２回 1
● 講義資料
▼ 講義予定
• 連立一次方程式に関するLU分解,および反復法
● 前回の講義のまとめ
★ ２階線形常微分方程式の境界値問題
• u: [0,1]−→Rに対する微分方程式の境界値問題
u′′(x) +cu(x) =f(x), x∈(0,1),
u(0) =u(1) = 0 (1)
を離散変数法を用いて数値的に解くことを考える.
• 定義域の区間[0,1]を N 等分して, xk =kh, (h= 1/N)での厳密解の値u(xk)の近似値を uk と書く. この時,２階微分u′′(x)を
u′′(x)∼ u(x+h)−2u(x) +u(x−h) h2
と差分化する. これは,
u′(x)∼u(x+h)−u(x) h
∼u(x)−u(x−h) h
u′′(x)∼u′(x+h)−u′(x) h
と考えて導出したものである.
• したがって,方程式を差分化すると,
u(xk+1)−2u(xk) +u(xk−1)
h2 +cu(xk) =f(xk) uk+1−2uk+uk−1
h2 +cuk =fk
となる. ここで,この式はk= 1, . . . , N−1 に対して成り立つ式である.
• ここで,境界条件u(0) =u(1) = 1は,
u0=uN = 0 と書き直せるので,解くべき式は
u2−2u1
h2 +cu1=f1, uk+1−2uk+uk−1
h2 +cuk=fk, k= 2, . . . , N−2,
−2uN−1+uN−2
h2 +cuN−2=fN−2
(2)
となり,未知ベクトルU = (uk)に対する連立一次方程式となる.
Jan. 15, 2014, Version: 1.0 [email protected]
2013年度・後期・数理解析・計算機数学３・第１２回 2
• いま,
A=
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
と書くと, (2) は
1
h2A+cE
U =F (3)
と書き直すことができる. したがって,２階線形常微分方程式の境界値問題(1) を離散変数 法で解くことは,連立一次方程式(3)を解くことに帰着できた.
▼ 消去法
★ はき出し法
• 以下,A∈Mn(R),b,x∈Rn として,連立一次方程式 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
Jan. 15, 2014, Version: 1.0 [email protected]
2013年度・後期・数理解析・計算機数学３・第１２回 3
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−ii1)を左からかける」または,「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
Jan. 15, 2014, Version: 1.0 [email protected]
2013年度・後期・数理解析・計算機数学３・第１２回 4
は以下のようにしてO(n2)の計算量で解くことができる. いま,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}nj=i がすべて0となっている場合である. この場合には,i 列 目はi−1列以前の列の一次結合で書けることとなり, Aが正則であるという仮定と矛盾す る. よって, Gaussの消去法はA が正則であれば,必ず終了する.
• Gaussの消去法の計算量はO(n3)であるが,ステップ２のループ回数が掃出し法と比較して
半分となっているので,掃出し法のほぼ半分の計算量となる.
● 実習内容
1. （★★★）連立一次方程式
−2x+ y = 3
x−2y+ z= 4 y−2z= 5
をLU分解を用いて解くプログラムと, Jacobiの反復法, Gauss-Seidelの反復法で解くプロ グラムを書きなさい.
Jan. 15, 2014, Version: 1.0 [email protected]
