● 講義資料
▼ 講義予定
• 連立一次方程式の数値解法
– 掃出し法
– Gaussの消去法
– LU分解
● 前回の講義のまとめ
★ 2階線形常微分方程式の境界値問題
• u: [0,1]−→Rに対する微分方程式の境界値問題
u′′(x) +p(x)u(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 +p(xk)u(xk) =f(xk) uk+1−2uk+uk−1
h2 +pkuk =fk
となる. ここで,この式はk= 1, . . . , N−1 に対して成り立つ式である.
• ここで,境界条件u(0) =u(1) = 1は,
u0=uN = 0
と書き直せるので,解くべき式は u2−2u1
h2 +p1u1=f1, uk+1−2uk+uk−1
h2 +pkuk=fk, k= 2, . . . , N−2,
−2uN−1+uN−2
h2 +pN−2uN−2=fN−2
(2)
となり,未知ベクトルU = (uk)に対する連立一次方程式となる.
• いま,
A=
−2 1
1 −2 1
. .. ... ...
1 −2 1
1 −2
, P =
p1 0
0 p2 0
. .. ... . ..
0 pN−2 0 0 pN−1
,
と書くと, (2) は
1 h2A+P
U =F (3)
と書き直すことができる. したがって,2階線形常微分方程式の境界値問題(1) を離散変数 法で解くことは,連立一次方程式(3)を解くことに帰着できた.
• 以下では,特にp(x) =cの場合を考える. すると,連立一次方程式は 1
h2A+cI
U =F (4)
と書くことができ, 連立一次方程式(4)の可解性は, 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 に対して,連 立一次方程式(4)の解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に対して(4)は可解となる. 特 に,c <0ならば(4)は可解となる.
★ 連立一次方程式の数値解法
• 以下,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. (★★★)連立一次方程式
2x+ 2y+ z= 3 x+ y+ 2z= 4 y+ z= 5
をGauss-Jordanの消去法(掃出し法)で解くプログラムを書きなさい.
2. (★★)行列
A=
3 12 9
2 5 4
−1 4 2
の逆行列を求めなさい.
3. (★★★)連立一次方程式
3x+ 12y+ 9z= 3 2x+ 5 y+ 4z= 4
−x+ 3 y+ 2z= 5
をGaussの消去法で解くプログラムを書きなさい.
4. (★★★)連立一次方程式
2x+ 2y+ z= 3 x+ y+ 2z= 4 y+ z= 5
をGaussの消去法で解くプログラムを書きなさい.
5. (★★★)行列
A=
3 12 9
2 5 4
−1 4 2
のLU分解を(種々の方法で)求め,それを用いて, b=t(3,4,5),b=t(3,4,−5)に対して, 連立一次方程式Ax=bを解くプログラムを書きなさい.