● 講義資料
▼ 講義予定
• 連立一次方程式の数値解法
– Gaussの消去法
– LU分解
● 前回の講義のまとめ
★ 2階線形常微分方程式の境界値問題
• 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 と書く. この時,2階微分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)に対する連立一次方程式となる.
• いま,
A=
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
と書くと, (2) は
1
h2A+cE
U =F (3)
と書き直すことができる. したがって,2階線形常微分方程式の境界値問題(1) を離散変数 法で解くことは,連立一次方程式(3)を解くことに帰着できた.
• 連立一次方程式(3)の可解性は,h−2A+cI のすべての固有値が0であるか否かで判定可能 である. いま, Aの固有値(と固有ベクトル)は,
λk=−4 sin2(kπ
2N), Vk=
sin(πkN) sin(2πkN )
... sin((N−N2)πk) sin((N−N1)πk)
=
sin(jkπN )
である. したがって,h−2A+cI の固有値は,N h= 1を使えば, Λk,h=c−4 sin2(kπ
2N) =c−4 sin2(kπh 2 )
と書ける. よって,与えられたh >0に対して, Λk,h6= 0となるとき,任意のF に対して,連 立一次方程式(3)の解U を求めることができる.
• いま,g(h) =h−1sin(h)は,単調減少(hを小さくするほど大きくなる)かつ, limh→0g(h) = 1 を使うと,h∈(0, π)において
4 sin2(kπh/2)
h2 >(kπ)2 が成り立つ. したがって,h∈(0, π)ならば,
Λk,h=c−4 sin2(kπh
2 )< c−(kπ)2< c−π2 が成り立つ.
• 離散変数法で数値解を構成するときには,h >0 は十分小さくとるので, c < π2
が成り立てば, Λk,h<0となり,任意の(十分小さな)h >0に対して(3)は可解となる. 特 に,c <0ならば(3)は可解となる.
★ 連立一次方程式の数値解法
• 以下,A∈Mn(R),b,x∈Rn として,連立一次方程式 Ax=b
を数値的に解くことを考える. (A が正則でない場合には, 「正則でない」ことが判定でき ればよいことにする.)
• 連立一次方程式を数値的に解く方法としては,大きく分けて,消去法と反復法がある. はじめ に,消去法の中でももっとも基本的な「掃出し法」(「Gauss-Jordanの消去法」)を考える.
• 掃出し法とは, 「行基本変形」を繰り返して,行列 Aを単位行列に変形する方法である. こ こで,「行基本変形」は,以下の3種類の変形である.
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−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 であるときには,プログラムを停止する.
● 実習内容
1. (★★★)連立一次方程式
3x+ 12y+ 9z= 3 2x+ 5 y+ 4z= 4
−x+ 3 y+ 2z= 5
をGaussの消去法で解くプログラムを書きなさい.
2. (★★★)連立一次方程式
2x+ 2y+ z= 3 x+ y+ 2z= 4 y+ z= 5
をGaussの消去法で解くプログラムを書きなさい.
3. (★★★)行列
A=
3 12 9
2 5 4
−1 4 2
の LU分解を求め, それを用いて, b =t(3,4,5), b =t(3,4,−5) に対して, 連立一次方程式 Ax=b を解くプログラムを書きなさい.
4. (★★)行列
A=
2 2 1
1 1 2
0 1 1
の LU分解を求め, それを用いて, b =t(3,4,5), b =t(3,4,−5) に対して, 連立一次方程式 Ax=b を解くプログラムを書きなさい.
5. (★★★)定数係数2階線形常微分方程式の境界値問題
u′′(x) = sin(πx), u(0) =u(1) = 0, u: [0,1]−→R
を,適切に刻み幅をきめて, Gaussの消去法(LU分解), Jacobiの反復法, Gauss-Seidel 法, SOR法を用いて解きなさい.