『発展系の数値解析』に加えること

桂田祐史

1997 年, 2009 年 7 月 21 日, 10 月 22 日

(2009 年 9 月 30 日結構大きく書き換える。なかなか収束しないなあ…まあまあ理解している方だと思うが、きちんと書くのは大変だ。誰かさぼらず書いておいてくれたら良いのに、としみじみ思う。言い出しっぺの法則なのか…)

目次

1 前書き 2

2 LU 分解 2 2.1 何のための工夫か . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 LU 分解の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2.1 三角行列の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2.2 三角行列の性質 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.3 LU 分解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.3 LU 分解を用いた連立 1 次方程式の解法 . . . . . . . . . . . . . . . . . . . . . . 4 2.4 LU 分解の計算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.4.1 疎性が保存されるか . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5 LU 分解の実例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.6 連立 1 次方程式以外への応用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.6.1 行列式の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.6.2 実対称行列の正値性・負値性の判定、符号 . . . . . . . . . . . . . . . . 12 2.6.3 2 次形式の平方完成と Gauss の消去法 . . . . . . . . . . . . . . . . . . 12

2.7 LU 分解可能なための条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.8 プログラムの書き方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 行列法による安定性解析 21 3.1 差分解は「等比ベクトル列」 . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Jordan 標準形を用いた「等比マトリックス列」の解析 . . . . . . . . . . . . . 22 3.3 J のスペクトル解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.4 R のスペクトル解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.5 差分解の漸近挙動の解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.6 安定性の条件を満たす場合の差分解の収束証明 . . . . . . . . . . . . . . . . . . 28

A メモ (未整理) 28 A.1 「下三角行列の逆行列は下三角行列」 . . . . . . . . . . . . . . . . . . . . . . 28 A.2 Lagrange の方法による 2 次形式の標準形への変換 . . . . . . . . . . . . . . . . 29 B Changelog 31

1 前書き

『発展系の数値解析』 [2] は、放送大学の講義科目「応用数学」のテキストの内容として用意されたものだが、「応用数学」も終了したことだし、この辺でゼミ用のテキストとしての整備を始めよう、と考えた。

下書きと言うか、とりあえず書いたものは既に存在するものが多い。

1. LU 分解について説明する。

2. 境界条件が Neumann であったり、非同次である場合の手短な導入をする。

  3. プログラム heat1d-i-glsc.c, heat1d-n-glsc.c を読むという項を用意する。

4. 安定性を行列を用いて解析する方法の紹介

5. von Neumann の安定性解析

6. 2 次元の Laplacian の差分近似

(a) 長方形領域における Poisson 方程式

2 LU 分解

『発展系の数値解析』の「Gauss の消去法」の後に読むべきものである。

( 「発展系の数値解析」を読むのに並行して、 LU 分解を勉強するとき迷子にならないように、手短な案内をしよう。 )

2.1 何のための工夫か

熱方程式を陰解法で解くには、連立 1 次方程式を 1 ステップにつき 1 度解かねばならない。トータルでは、膨大な回数解くことになる。この連立 1 次方程式には

1. 係数行列はいつも同じ

2. 係数行列は疎行列 (0 でない成分は対角線の近くにしかない)

という大きな特徴がある。このことを念頭に、なるべく上手に計算する方法を考える (ここで説明する方法を採用せず、素朴なやり方で解こうとすると、解くことがほとんど非現実的なくらい効率が低くなってしまう)。

係数行列がいつも同じということから、逆行列を最初に計算しておくと良いと考えるかもしれないが、実はそれはひどく非効率的になって拙いやり方である。

『発展系の数値解析』で説明した Gauss の消去法を用いれば、係数行列が疎であるという性質を生かした効率的な計算ができるが、毎回行列の変形をしなければならないのが嫌味である。しかし、実際にその行列の変形をしてみれば分かるが、大部分の計算は毎回まったく同じである。うまく工夫すれば計算の大部分が省略できることが想像できるはずである。それを実現する方法が、以下に説明する係数行列の LU 分解である。

忘れないように

[figure]

• 逆行列は使わない (もし使うと、どうしようもないくらい非能率的になる)

• LU 分解は Gauss の消去法を工夫したものとみなせる^{a}

aLU 分解そのものは Gauss の消去法と無関係に定義され、 Gauss の消去法以外の方法で計算することも可能である。 Gauss の消去法は LU 分解を求めるための代表的なアルゴリズムである、というに過ぎない。

[figure]

2.2 LU 分解の定義

LU 分解の定義と基本的な性質のいくつかを (証明を省略して^{1}) 述べる。

2.2.1 三角行列の定義

対角線より上にあるすべての成分が 0 であるような行列を下三角行列と呼ぶ。すなわち L = (ℓ_{ij}) が下三角行列であるとは、

i < j ⟹ ℓ_{ij} = 0

が成り立つことをいう。

特に以下では対角成分がすべて 1 である下三角行列をしばしば用いる。このような行列を単位下三角行列と呼ぶ。

同様に対角線より下にあるすべての成分が 0 である行列を上三角行列とよび、そのうち対角成分がすべて 1 であるようなものを単位上三角行列と呼ぶ。

2.2.2 三角行列の性質

概観を得るために基本的な性質を証明抜きに述べておく (難しくはないから自力で証明にチャレンジしてもよい)。

下三角行列について得られる結果と並行した結果が上三角行列について得られるが、以下では省略する (下三角行列についてのみ説明する)。

n 次下三角行列 L = (ℓ_{ij}) について、

det L = ∏_{i=1}^{n} ℓ_{ii}.

すなわち、下三角行列の行列式は、すべての対角成分の積である。ゆえに、下三角行列が正則であるためには、すべての対角成分が 0 でないことが必要十分である。特に単位下三角行列は正則である。

下三角行列同士の和、差、積、また (存在する場合の) 逆行列は下三角行列である^{2}。正則な下三角行列全体は乗算に関して群をなす。特に単位下三角行列全体はその部分群と

なる。 2.2.3 LU 分解

1大部分は簡単であるから、自分で証明を試みるとよい。

2証明は、例えば桂田 [4] の第 4 章 (2009 年 9 月 24 日現在) に書いてある。一つの証明のアウトラインを書いておく: A = (a_{ij}), B = (b_{ij}) を共に下三角行列で、 AB = I, ∀i a_{ii} ≠ 0 が成り立つとすると、 b_{ij} が簡単に求められる。なお、後の余談 2.1 も見よ。

行列 A に対して、

A = LU, L は下三角行列, U は上三角行列

をみたす L, U が存在するとき、 L と U の組を A の LU 分解とよび、 L と U の組を求めることを A を LU 分解する、という。

以下、後の内容を予告する。

任意に与えられた行列 A に対して、その LU 分解はいつも存在するとは限らない (定理 2.7 で示すように、 A が正則である場合は、 A のすべての主座小行列式^{3}が ≠ 0 であることが LU 分解が存在するための必要十分条件である)。しかし A が正則であれば、適当な置換行列^{4} P を左からかけた P A を LU 分解することができる (P A は A の行ベクトルを入れ替えたもの)。

正則行列 A が LU 分解できるとき、分解の仕方を

(1) L の対角成分はすべて 1 (すなわち L は単位下三角行列)

(2) U の対角成分はすべて 1 (すなわち U は単位上三角行列)

のいずれかに限定すると、分解は一意的に定まる (命題 2.9)。

2.3 LU 分解を用いた連立 1 次方程式の解法

n 次正則行列 A が A = LU と LU 分解されているとき、連立 1 次方程式

Ax = b

は少ない計算量で解くことが出来る。以下、このことを説明する。

Ax = b (⇔ LUx = b)

は

Ly = b, Ux = y

という二つの問題に分解される。

まず Ly = b は

ℓ_{11}y_{1}  = b_{1}
ℓ_{21}y_{1} +ℓ_{22}y_{2}  = b_{2}
ℓ_{31}y_{1} +ℓ_{32}y_{2} +ℓ_{33}y_{3}  = b_{3}
⋮ ... ⋮
ℓ_{n1}y_{1} +ℓ_{n2}y_{2} +ℓ_{n3}y_{3} ⋯ +ℓ_{nn}y_{n} = b_{n}

3A の k 次の主座小行列 (首座小行列とも呼ぶ) とは、 A の最初の k 行 k 列から得られる k 次の正方行列のことをいう。その行列式のことを k 次の主座行列式と呼ぶ。

4知らないという人がいたので、簡単に説明しておく。 i ≠ j のとき、 P_{ij} := ∑_{k≠i,j} E_{kk} + E_{ij} + E_{ji} とおくと、 P_{ij}A は A の第 i 行と第 j 行が入れ替えた行列となる。いくつかの P_{ij} の積として得られる行列 P を置換行列という。

ということであり、これは上から順に

y_{1} = b_{1}/ℓ_{11},
y_{2} = (b_{2} − ℓ_{21}y_{1}) /ℓ_{22},
y_{3} = (b_{3} − ℓ_{31}y_{1} − ℓ_{32}y_{2}) /ℓ_{33},
⋯
y_{i} = (b_{i} − ∑_{j=1}^{i−1} ℓ_{ij}y_{j}) /ℓ_{ii},
⋯
y_{n} = (b_{n} − ∑_{j=1}^{n−1} ℓ_{nj}y_{j}) /ℓ_{nn}

と解くことが出来る (A が正則であると仮定したことから、 L も正則で、すべての i について ℓ_{ii} ≠ 0 が成り立つことに注意)。

これを計算するには、 1 + 2 + ⋯ + n = n(n + 1)/2 回の乗除算で十分である^{5}。

余談 2.1 (実は「下三角行列の逆行列は下三角」の別証になっている) 上の結果から、

y_{1} は b_{1} の線形結合,
y_{2} は b_{1}, b_{2} の線形結合,
⋮
y_{i} は b_{1}, b_{2}, ⋯, b_{i} の線形結合,
⋮
y_{n−1} は b_{1}, b_{2}, ⋯, b_{n−1} の線形結合,
y_{n} は b_{1}, b_{2}, ⋯, b_{n} の線形結合

であることが分かる。ゆえに y = (下三角行列)b. これは L の逆行列が下三角行列であることを示している。

同様に Ux = y は

u_{11}x_{1} ⋯ +u_{1,n−2}x_{n−2} +u_{1,n−1}x_{n−1} +u_{1n}x_{n}  = y_{1},
... ⋮ ⋮
u_{n−2,n−2}x_{n−2} +u_{n−2,n−1}x_{n−1} +u_{n−2,n}x_{n} = y_{n−2},
u_{n−1,n−1}x_{n−1} +u_{n−1,n}x_{n} = y_{n−1},
u_{n,n}x_{n}  = y_{n}

ということである。やはり A が正則という仮定から、 u_{ii} ≠ 0 (i = 1, …, n) が導かれることに注意すると、この連立方程式は、下から順に

5計算にどれくらい手間がかかるかを計るために、基本的な演算の回数を数えるという方法がある。しばしば、四則演算のそれぞれについて数える代りに、乗除算の回数だけを数えて目安にする、という手段が採用される (4 次元ベクトルよりは 1 つの数値が分かりやすい)。もう少し詳しいことが知りたければ、例えば桂田 [4] を見よ。

x_{n} = y_{n}/u_{nn},
x_{n−1} = (y_{n−1} − u_{n−1,n}x_{n}) /u_{n−1,n−1},
x_{n−2} = (y_{n−2} − u_{n−2,n−1}x_{n−1} − u_{n−2,n}x_{n}) /u_{n−2,n−2},
⋮  ⋮
x_{i} = (y_{i} − ∑_{j=i+1}^{n} u_{ij}x_{j}) /u_{ii},
⋮  ⋮
x_{1} = (y_{1} − ∑_{j=2}^{n} u_{1j}x_{j}) /u_{11}

と解くことができる。これも計算するには、 1 + 2 + ⋯ + n = n(n + 1)/2 回の乗除算で十分である。まとめると、 n(n + 1)/2 + n(n + 1)/2 = n(n + 1) 回の乗除算で連立 1 次方程式が解けるこ

とになる^{6}。これは連立 1 次方程式を「普通に」解く場合に、 n^{3} に比例する回数の乗除算が必要なことと比較して、 (n が大きな場合は) かなり少ない回数となる。

普通に解く O(n^{3}) v.s. LU 分解 (が与えられていてそれ) を使った場合 n^{2} 程度

2.4 LU 分解の計算法

実は、 Gauss の消去法は LU 分解を計算していることに相当する。このことを理解すれば、 LU 分解の計算法を知っていることになる。 Gauss の消去法が最後まで実行できるためには、消去の過程で現れる対角成分が 0 にならないことが必要十分であるが、この項では、それが常に満たされるとして議論する (その条件を吟味することは後の 2.7 で行う)。

Gauss の消去法では、「行に関する基本変形」で係数行列を上三角行列に変形したが、「行に関する基本変形」は基本行列を左からかけることに相当する (このことは基本変形を説明してある線形代数の教科書ならば、必ず解説してあるはずであるが、桂田 [4] の第 4 章にも書いておいた)。

n 次正方行列 A = (a_{ij}) の (k, k) 成分 a_{kk} で (i, k) 成分 (i > k) を掃き出すには、第 i 行から、第 k 行の q_{ik} := a_{ik}/a_{kk} 倍を引く、という操作をするわけだが、これは

L_{ik} := I − q_{ik}E_{ik} = ( \begin{matrix}1 & & & & & & & & & & & & 0 \\ & . & . & . \\ & & & & 1 \\ & & & & & . & . & . \\ & & & & −q_{ik} & & & & 1 \\ & & & & & & & & & . & . & . \\ 0 & & & & & & & & & & & & 1 \end{matrix} ) \begin{matrix}⟵ k 行目 \\ ⟵ i 行目 \end{matrix}

6細かい話をすると、 L が単位下三角行列である場合は、割り算の回数が n 回減って、 n^{2} 回の乗除算で解けることになる。逆行列を知っている場合、連立 1 次方程式は n^{2} 回の掛け算で解けるが、それと互角であることが分かる。

を左からかけることで実現される (ここで I は n 次の単位行列, E_{ik} は第 (i, k) 成分のみ 1 で、他の成分は 0 であるような n 次正方行列)。すると、 (k, k) 成分 a_{kk} で、第 k 列の k 行目以降 ((k + 1, k), (k + 2, k), ⋯, (n, k) 成分) を掃き出す操作は、

L_{k} := L_{n,k} ⋯ L_{k+2,k}L_{k+1,k}

を左からかけることで実現されることが分かる。後で逆行列が必要になるので調べておこう。まず

L_{ik}^{−1} = I + q_{ik}E_{ik} = ( \begin{matrix}1 & & & & & & & & & & & & 0 \\ & . & . & . \\ & & & & 1 \\ & & & & & . & . & . \\ & & & & q_{ik} & & & & 1 \\ & & & & & & & & & . & . & . \\ 0 & & & & & & & & & & & & 1 \end{matrix} ) \begin{matrix}⟵ k 行目 \\ ⟵ i 行目 \end{matrix}

である。これは計算で示すことも出来るし^{7}、行列の表す基本変形の意味から考えても明らかである。

次に

(1)  L_{k}^{−1} = L_{k+1,k}^{−1}L_{k+2,k}^{−1} ⋯ L_{n,k}^{−1} = ( \begin{matrix}1 \\ & . & . & . \\ & & & & 1 \\ & & & & q_{k+1,k} & 1 \\ & & & & q_{k+2,k} & & 1 \\ & & & & ⋮ & & & ... \\ & & & & q_{n,k} & & & & 1 \end{matrix} )

である。これも基本変形の意味から明らかである。

上の操作 (L_{k} を左からかける) を k = 1, 2, ⋯, n − 1 の順に行うと、 A の対角線よりも下の部分は掃き出されて、上三角行列になる:

L_{n−1}L_{n−2} ⋯ L_{2}L_{1}A = 上三角行列=: U.

簡単のため

L := L_{n−1}L_{n−2} ⋯ L_{2}L_{1}

とおけば、

(2)  LA = U.

各 L_{k} が単位下三角行列であるから、 L も単位下三角行列である。ゆえに L は正則で、逆行列 L := L^{−1} (= L_{1}^{−1}L_{2}^{−1} ⋯ L_{n−1}^{−1}) も単位下三角行列である。ゆえに (2) に左から L をかけると、

A = LU

となり、 A の LU 分解が得られた。

7E_{ij}E_{kℓ} = δ_{jk}E_{iℓ} という公式が成り立つことに注意すれば、 i ≠ k であれば (I − q_{ik}E_{ik})(I + q_{ik}E_{ik}) = I − q_{ik}^{2}E_{ik}^{2} = I − q_{ik}^{2}δ_{ik}E_{ik} = I − q_{ik}^{2}O = I.

実はより具体的に

(3)  L = ( \begin{matrix}1 & & & & 0 \\ q_{21} & 1 \\ q_{31} & q_{32} & 1 \\ ⋮ & ⋮ & ... & ... \\ q_{n1} & q_{n2} & ⋯ & q_{n,n−1} & 1 \end{matrix} )

であることが分かる (これをもっと前に持って来るのが良いと思うが、はてどうやろうか…)。 (3) の証明 (1) の証明と同様に基本変形の意味を考えても分かるが、ここでは地道な計算で示してみよう。

各 k に対して

L_{1}^{−1} ⋯ L_{k}^{−1} = ( \begin{matrix}1 & & 0 \\ q_{21} & ... & & & 0 \\ ⋮ & ... & 1 \\ ⋮ & & q_{k+1,k} & 1 & & 0 \\ ⋮ & & ⋮ & & ... \\ q_{n1} & ⋯ & q_{n,k} & 0 & & 1 \end{matrix} )

となることを帰納法で示す。

k = 1 のとき成り立つことは明らかである (式 (1) から、 L_{1}^{−1} がこの形をしていることがすぐ分かる)。

次に k のとき成立すると仮定する。

L_{1}^{−1} ⋯ L_{k}^{−1}L_{k+1}^{−1} = ( \begin{matrix}1 & & 0 \\ q_{21} & ... & & & 0 \\ ⋮ & ... & 1 \\ ⋮ & & q_{k+1,k} & 1 & & & 0 \\ ⋮ & & q_{k+2,k} & & 1 \\ ⋮ & & ⋮ & & & ... \\ q_{n1} & ⋯ & q_{n,k} & 0 & & & 1 \end{matrix} ) ( \begin{matrix}1 & & 0 \\ & ... & & & 0 \\ 0 & & 1 \\ & & & 1 & & & 0 \\ & & & q_{k+2,k+1} & 1 \\ & 0 & & ⋮ & & ... \\ & & & q_{n,k+1} & 0 & & 1 \end{matrix} )

= ( \begin{matrix}A & O \\ C & I \end{matrix} ) ( \begin{matrix}I & O \\ O & D \end{matrix} ) = ( \begin{matrix}A & O \\ C & D \end{matrix} )

= ( \begin{matrix}1 & & 0 \\ q_{21} & ... & & & 0 \\ ⋮ & ... & 1 \\ ⋮ & & q_{k+1,k} & 1 & & & 0 \\ ⋮ & & q_{k+2,k} & q_{k+2,k+1} & 1 \\ ⋮ & & ⋮ & ⋮ & & ... \\ q_{n1} & ⋯ & q_{n,k} & q_{n,k+1} & 0 & & 1 \end{matrix} ) .

これは k + 1 のときも成立することを示している。注意 2.1 (実は後で必要のない事実だが良くある誤解を正す)L_{k} は L_{k}^{−1} と似ていて、

L_{k} = ( \begin{matrix}1 \\ & . & . & . \\ & & & & 1 \\ & & & & −q_{k+1,k} & 1 \\ & & & & −q_{k+2,k} & & 1 \\ & & & & ⋮ & & & ... \\ & & & & −q_{n,k} & & & & 1 \end{matrix} )

となる。ところが、

(これは誤り) L = ( \begin{matrix}1 & & & & 0 \\ −q_{21} & 1 \\ −q_{31} & −q_{32} & 1 \\ ⋮ & ⋮ & ... & ... \\ −q_{n1} & −q_{n2} & ⋯ & −q_{n,n−1} & 1 \end{matrix} )

は成り立たない！筆者自身が一時混乱したので (成り立つはずだと勘違いして証明を試みたりした) 参考情報として残しておく。

2.4.1 疎性が保存されるか

A が A = LU と LU 分解可能なとき、 A の疎性は L と U に「遺伝する」、つまり A が疎行列ならば、 L と U も疎行列であることは比較的容易に分かる (より具体的には「バンド幅」 — まだ説明していない概念^{8}だが — は増えない)。

これに対して、 A が疎行列であっても、 A^{−1} が疎行列になるとは限らない (比較的有名であるので例は後回しにする — なお、次の例 2.2 もその例だと言えないこともない)。 L^{−1} も疎行列になるとは限らないことを例で示しておこう。

例 2.2 (L が疎行列であっても L^{−1} はそうとは限らない)下三角行列を係数行列に持つ連立 1 次方程式

( \begin{matrix}1 & & & & 0 \\ a_{2} & 1 \\ & a_{3} & ... \\ & & ... & 1 \\ 0 & & & a_{n} & 1 \end{matrix} ) ( x_{1} ; x_{2} ; x_{3} ; ⋮ ; x_{n} ) = ( b_{1} ; b_{2} ; b_{3} ; ⋮ ; b_{n} )

は、既に説明した手順によって、

x_{1} = b_{1},
x_{2} = b_{2} − a_{2}x_{1},
(4)  x_{3} = b_{3} − a_{3}x_{2},
⋮  ⋮
x_{n} = b_{2} − a_{n}x_{n−1}

8例えば、桂田 [4] などを見よ。

と (簡単かつ少ない手間で) 解ける。 x_{i} を a_{j}, b_{k} のみで表すと

x_{1} = b_{1},
x_{2} = b_{2} − a_{2}x_{1} = b_{2} − a_{2}b_{1},
x_{3} = b_{3} − a_{3}x_{2} = b_{3} − a_{3}b_{2} + a_{3}a_{2}b_{1},
x_{4} = b_{4} − a_{4}x_{3} = b_{4} − a_{4}b_{3} + a_{4}a_{3}b_{2} − a_{4}a_{3}a_{2}b_{1},
⋮  ⋮
x_{n} = b_{2} − a_{n}x_{n−1} = b_{n} − a_{n}b_{n−1} + a_{n}a_{n−1}b_{n−2} − ⋯ + (−1)^{n−1}a_{n}a_{n−1} ⋯ a_{2}b_{1}.

この結果は

( \begin{matrix}1 & & & & 0 \\ a_{2} & 1 \\ & a_{3} & ... \\ & & ... & 1 \\ 0 & & & a_{n} & 1 \end{matrix} )^{−1} = ( \begin{matrix}1 & & & & & & 0 \\ −a_{2} & 1 \\ a_{3}a_{2} & −a_{3} & 1 \\ −a_{4}a_{3}a_{2} & a_{4}a_{3} & −a_{4} & 1 \\ ⋮ & ... & ... & ... & ... \\ ⋮ & & ... & a_{n−1}a_{n−2} & −a_{n−1} & 1 \\ (−1)^{n−1}a_{n} ⋯ a_{2} & ⋯ & & −a_{n}a_{n−1}a_{n−2} & a_{n}a_{n−1} & −a_{n} & 1 \end{matrix} )

を示している。

この例の問題で L を記憶しておいて、 (4) で解を求めるやり方と、 L^{−1} を計算して b = (b_{i}) にかけて解を求めるやり方を比べてみるのは有益である。前者は、 n に比例する計算量ですむ (疎性を有効に活かせた) が、後者は、行列 L^{−1} とベクトル b の掛け算だけで、 n^{2} に比例する計算量となってしまう (疎性を活かせなかった)。

2.5 LU 分解の実例

Gauss の消去法の例としてあげた

( \begin{matrix}2 & 3 & −1 & 5 \\ 4 & 4 & −3 & 3 \\ −2 & 3 & −1 & 1 \end{matrix} ) → ( \begin{matrix}2 & 3 & −1 & 5 \\ 0 & −2 & −1 & −7 \\ 0 & 6 & −2 & 6 \end{matrix} ) → ( \begin{matrix}2 & 3 & −1 & 5 \\ 0 & −2 & −1 & −7 \\ 0 & 0 & −5 & −15 \end{matrix} )

を考えてみよう。もちろん係数行列の変形だけ取り出すと

( \begin{matrix}2 & 3 & −1 \\ 4 & 4 & −3 \\ −2 & 3 & −1 \end{matrix} ) → ( \begin{matrix}2 & 3 & −1 \\ 0 & −2 & −1 \\ 0 & 6 & −2 \end{matrix} ) → ( \begin{matrix}2 & 3 & −1 \\ 0 & −2 & −1 \\ 0 & 0 & −5 \end{matrix} )

となる。

q_{21} = \frac{4}{2} = 2, q_{31} = \frac{−2}{2} = −1,
q_{32} = \frac{6}{−2} = −3

が分かるので、

L = ( \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ −1 & −3 & 1 \end{matrix} )

であり、もちろん (上で見えているように)

U = ( \begin{matrix}2 & 3 & −1 \\ 0 & −2 & −1 \\ 0 & 0 & −5 \end{matrix} )

である。実際

LU = ( \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ −1 & −3 & 1 \end{matrix} ) ( \begin{matrix}2 & 3 & −1 \\ 0 & −2 & −1 \\ 0 & 0 & −5 \end{matrix} ) = ( \begin{matrix}2 & 3 & −1 \\ 4 & 4 & −3 \\ −2 & 3 & −1 \end{matrix} ) = A.

(加筆) §2.4 は読むのが面倒なのか、読めなかった人がいるので、この実例で説明してみる。

L = (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 3 & 1 \end{matrix}) (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{matrix}) (\begin{matrix}1 & 0 & 0 \\ −2 & 1 & 0 \\ 0 & 0 & 1 \end{matrix})

であるから、

L = L^{−1} = (\begin{matrix}1 & 0 & 0 \\ −2 & 1 & 0 \\ 0 & 0 & 1 \end{matrix})^{−1} (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{matrix})^{−1} (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 3 & 1 \end{matrix})^{−1}

= (\begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ −1 & 0 & 1 \end{matrix}) (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & −3 & 1 \end{matrix}) · (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix})

= (\begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ −1 & 0 & 1 \end{matrix}) (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & −3 & 1 \end{matrix})

= (\begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ −1 & −3 & 1 \end{matrix}) = (\begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ −1 & −3 & 1 \end{matrix}) .

2.6 連立 1 次方程式以外への応用

2.6.1 行列式の計算行列式を計算したい場合にも、 LU 分解はしばしば最も効率的な計算法となる。 A が単位下三角行列 L と上三角行列 U の積に分解できたとする:

A = LU, L = ( \begin{matrix}1 & & & & 0 \\ & . & . & . \\ ∗ & & & & 1 \end{matrix} ) , U = ( \begin{matrix}u_{11} & & & & ∗ \\ & . & . & . \\ 0 & & & & u_{nn} \end{matrix} ) .

このとき

det A = det(LU) = det L det U = 1 · ∏_{i=1}^{n} u_{ii} = ∏_{i=1}^{n} u_{ii}.

2.5 の例では、

det A = 2 · (−2) · (−5) = 20.

(補足: 対角線の下の掃き出しをするのに行の交換が必要だった場合。 P A = LU, P は置換行列で det P = (−1)^{ℓ} (ℓ は行交換の回数) ということなので、 det A = (−1)^{ℓ} det U = (−1)^{ℓ} ∏_{i=1}^{n} u_{ii}.)

2.6.2 実対称行列の正値性・負値性の判定、符号

(どうも「符号を授業で習っていない」ことが少なくないようで、困ったものである。正値性・負値性の判定は重要であるから、それだけ切り離して読めるようにすべきかもしれない。 )

実対称行列の符号^{9}(正の固有値の個数と負の固有値の個数の組) を求めるにも、 LU 分解はしばしば有効である。

[figure]

命題 2.3 A が単位下三角行列 L と上三角行列 U の積に分解できたとする:

A = LU, L = ( \begin{matrix}1 & & & & 0 \\ & . & . & . \\ ∗ & & & & 1 \end{matrix} ) , U = ( \begin{matrix}u_{11} & & & & ∗′ \\ & . & . & . \\ 0 & & & & u_{nn} \end{matrix} ) .

このとき

A の符号= “{u_{ii}} のうちの正の個数、負の個数” .
[figure]

証明 (L^{−1})^{T} =: U とおくと、

U^{T} AU = UU

であるから、 A の符号は UU の符号と等しい。ところが、

U = ( \begin{matrix}1 & & & & ∗′′ \\ & . & . & . \\ 0 & & & & 1 \end{matrix} )

であるから、

UU = ( \begin{matrix}u_{11} & & & & ∗′′′ \\ & . & . & . \\ 0 & & & & u_{nn} \end{matrix} )

であり、この行列の符号は {u_{ii}} のうちの正であるものの個数、負であるものの個数の組に他ならない。

次項で見るように、実対称行列が正値であれば、 Gauss の消去法により、必ず LU 分解が可能である。従って、上に書いた変形は必ず可能である。つまり、与えられた実対称行列が正値であるかどうかを判定するには、 (ピボット選択無しの) Gauss の消去法の前進消去法が最後まで実行でき、その結果できる上三角行列の対角成分がすべて正であるかどうかを見れば良い。 (途中で対角成分が 0 になったり負になったりしたら、正値でないことが分かる)。行列が負値であることの判定も同様である。

2.6.3 2 次形式の平方完成と Gauss の消去法

実対称行列の正値性、負値性の判定に、 Gauss の消去法が利用できることに唐突な印象を持ったかも知れないが、実は 2 次式に対しての基本操作である「平方完成」は Gauss の消去法と関係がある。このことを見てみよう。

9符号の定義とか、 Sylvester の慣性律とか、線形代数で習うべきことだと思うが、残念ながら授業では省略されることもあるようである。定番の齋藤 [5], 佐武 [6] などを見よ。

A = (a_{ij}) が実対称行列であるとする。 a_{11} ≠ 0 であれば

a_{11} (x_{1} + ∑_{j=2}^{n} \frac{a_{1j}}{a_{11}} x_{j})^{2} = a_{11}x_{1}^{2} + 2x_{1} ∑_{j=2}^{n} a_{1j}x_{j} + \frac{1}{a_{11}} (∑_{j=2}^{n} a_{1j}x_{j})^{2}

= a_{11}x_{1}^{2} + x_{1} ∑_{j=2}^{n} a_{1j}x_{j} + x_{1} ∑_{i=2}^{n} a_{i1}x_{i} + \frac{1}{a_{11}} (∑_{i=2}^{n} a_{i1}x_{i}) (∑_{j=2}^{n} a_{1j}x_{j})

= a_{11}x_{1}^{2} + ∑_{j=2}^{n} a_{1j}x_{1}x_{j} + ∑_{i=2}^{n} a_{i1}x_{i}x_{1} + ∑_{i,j=2}^{n} \frac{a_{ij}a_{1j}}{a_{11}} x_{i}x_{j}

= a_{11}x_{1}^{2} + ∑_{j=2}^{n} a_{1j}x_{1}x_{j} + ∑_{i=2}^{n} a_{i1}x_{i}x_{1} + ∑_{i,j=2}^{n} a_{ij}x_{i}x_{j} + ∑_{i,j=2}^{n} (\frac{a_{ij}a_{1j}}{a_{11}} − a_{ij}) x_{i}x_{j}

= ∑_{i,j=1}^{n} a_{ij}x_{i}x_{j} + ∑_{i,j=2}^{n} (\frac{a_{ij}a_{1j}}{a_{11}} − a_{ij}) x_{i}x_{j}

であるから、

∑_{i,j=1}^{n} a_{ij}x_{i}x_{j} = a_{11} (x_{1} + ∑_{j=2}^{n} \frac{a_{1j}}{a_{11}} x_{j})^{2} + ∑_{i,j=2}^{n} (a_{ij} − \frac{a_{i1}a_{1j}}{a_{11}} ) x_{i}x_{j}.

これは Gauss の消去法の前進消去過程と同じである。与えられた 2 次形式が、平方完成を繰り返すことで標準形に変換出来る場合には (齋藤 [5], pp. 157–158 に載っている Lagrange の方法の特別な場合)、 Gauss の消去法で計算可能である。実際、実対称行列 A が A = LU と分解できたとき、

d_{i} := u_{ii}, u′_{ij} := \frac{u_{ij}}{d_{i}} ,

D := diag(d_{1}, …, d_{n}), U′ := (u′_{ij})

とおくと、

A = LDU′

となる (LDU 分解)。 x′ = U′x, すなわち

x′_{1} = u_{11}x_{1} + u_{12}x_{2} + ⋯ + u_{1,n−1}x_{n−1} + u_{1n}x_{n},
x′_{2} =  u_{22}x_{2} + ⋯ + u_{2,n−1}x_{n−1} + u_{2n}x_{n},
⋮ ⋮
x′_{n−1} =  u_{n−1,n−1}x_{n−1} + u_{n−1,n}x_{n},
x′_{n} =  u_{n,n}x_{n}

とおくと、

(Ax, x) = d_{1}x′_{1}^{2} + d_{2}x′_{2}^{2} + ⋯ + d_{n}x′_{n}^{2}.

2.7 LU 分解可能なための条件

この項では「正則行列が LU 分解可能であるためには、その行列のすべての主座小行列式が 0 にならないことが必要十分であり、そのとき Gauss の消去法の前進消去過程によって LU 分解が得られる」という定理を証明するのが目標である。記号の約束以下、行列 A に対して、 A^{(j)} で、 j 次の主座小行列を表す:

A^{(j)} = (\begin{matrix}a_{11} & ⋯ & a_{1j} \\ ⋮ & & ⋮ \\ a_{j1} & ⋯ & a_{jj} \end{matrix}) .

蛇足ながら、数値計算言語の MATLAB では、 a という行列に対して、 a(1:j,1:j) という式で j 次の主座小行列を表せる。

[figure]

命題 2.4 (主座小行列式 ≠ 0 の必要性)n 次正方行列 A に対して、 LA = U を満たす下三角行列 L, 上三角行列 U が存在するならば、

∀k ∈ {1, …, n} L^{(k)}A^{(k)} = U^{(k)}, det L^{(k)} det A^{(k)} = det U^{(k)}.

特に正則行列 A が LU 分解を持つならば、

∀k ∈ {1, …, n} det A^{(k)} ≠ 0. [figure]

証明行列を k 行, k 列のところで線を引いてブロックわけする:

A = ( \begin{matrix}A^{(k)} & ∗ \\ ∗ & ∗ \end{matrix} ) , L = ( \begin{matrix}L^{(k)} & O \\ ∗ & ∗ \end{matrix} ) , U = ( \begin{matrix}U^{(k)} & ∗ \\ O & ∗ \end{matrix} ) .

(L の右上のブロック、 U の左下のブロックは零行列ということを主張している。 )

LA = ( \begin{matrix}L^{(k)} & O \\ ∗ & ∗ \end{matrix} ) ( \begin{matrix}A^{(k)} & ∗ \\ ∗ & ∗ \end{matrix} ) = ( \begin{matrix}L^{(k)}A^{(k)} & ∗ \\ ∗ & ∗ \end{matrix} )

であるから、

L^{(k)}A^{(k)} = U^{(k)}, det L^{(k)} det A^{(k)} = det U^{(k)}.

正則 A が LU 分解 A = LU を持つとする。 0 ≠ det A = det L det U より、 L と U も正則で

ある。 det U = ∏_{i=1}^{n} u_{ii} より、 u_{ii} ≠ 0 (i = 1, …, n). L の逆行列を L とすると、 L は下三角で、

LA = U が成り立つから、前半を用いて、

∀k ∈ {1, …, n} det L^{(k)} det A^{(k)} = det U^{(k)} = ∏_{i=1}^{k} u_{ii} ≠ 0.

これから det A^{(k)} ≠ 0.

逆に det A^{(k)} ≠ 0 (k = 1, …, n) ならば A は LU 分解できることを示すには、 Gauss の消去法によって具体的に LU 分解を構成してみせる。

[figure]

補題 2.5 n 次正方行列 A が

∃k ∈ {1, 2, ⋯, n − 1} ∀j ∈ {1, 2, ⋯, k} det A^{(j)} ≠ 0

を満たすならば、 A について Gauss の消去法の前進消去過程は k 段まで実行できる。特に n 次単位下三角行列 L が存在して

LA = ( \begin{matrix}u_{1} & & * \\ & ... & & * \\ 0 & & u_{k} \\ & 0 & & * \end{matrix} )

という形になる。

[figure]

証明 k についての帰納法による。 k = 1 のとき、 0 ≠ det A^{(1)} = a_{11} より、第 1 列の対角線より下を掃き出すことができる。実際

L := ( \begin{matrix}1 & & & 0 \\ −q_{21} & 1 \\ ⋮ & & ... \\ −q_{n1} & 0 & & 1 \end{matrix} ) , q_{j1} := \frac{a_{j1}}{a_{11}} (2 ≤ j ≤ n)

とおくと、

LA = ( \begin{matrix}a_{11} \\ 0 \\ ⋮ & * \\ 0 \end{matrix} ) .

ゆえに k = 1 のとき確かに成り立つ。

k = ℓ のとき成り立つと仮定して、 k = ℓ + 1 とする。 A が det A^{(j)} ≠ 0 (1 ≤ ∀j ≤ k) を満たすと仮定する。特に det A^{(j)} ≠ 0 (1 ≤ ∀j ≤ ℓ) であり、帰納法の仮定より Gauss の消去法の前進消去過程は ℓ 段まで実行できる。すなわち行に関する基本変形を表す下三角の基本行列の積である下三角行列 L が存在して、

LA = ( \begin{matrix}u_{1} & & * \\ & ... & & * \\ 0 & & u_{ℓ} \\ & 0 & & A′ \end{matrix} ) .

これから

L^{(ℓ+1)}A^{(ℓ+1)} = (LA)^{(ℓ+1)} = ( \begin{matrix}u_{1} & & * \\ & ... & & * \\ 0 & & u_{ℓ} \\ 0 & ⋯ & 0 & u_{ℓ+1} \end{matrix} ) , u_{ℓ+1} := A′ の (1, 1) 成分.

このとき

u_{1} ⋯ u_{ℓ}u_{ℓ+1} = det (L^{(ℓ+1)}A^{(ℓ+1)}) = det L^{(ℓ+1)} det A^{(ℓ+1)} = 1 · det A^{(ℓ+1)} = det A^{(ℓ+1)}.

仮定より det A^{(ℓ+1)} ≠ 0 であるから、 u_{ℓ+1} ≠ 0. ゆえに第 (ℓ + 1) 列の対角線より下も掃き出すことができる。すなわち Gauss の消去法の前進消去過程は (ℓ + 1) 段まで実行できる。

[figure]

命題 2.6 (主座小行列式 ≠ 0 の十分性)n 次正方行列 A が

det A^{(k)} ≠ 0 (k = 1, …, n − 1)

を満たすならば、 A について Gauss の消去法の前進消去過程は n − 1 段まで実行できる。すなわち n 次単位下三角行列 L と上三角行列 U が存在して

LA = U.

特に A が正則 (det A^{(n)} ≠ 0) であれば、

A = LU

を満たす単位下三角行列 L と上三角行列 U = (u_{ij}) が存在する。このとき、

det A^{(k)} = det U^{(k)} = ∏_{i=1}^{k} u_{ii} (k = 1, …, n),

u_{11} = det A^{(1)}, u_{kk} = \frac{det A^{(k)}}{det A^{(k−1)}} (k = 2, …, n). [figure]

証明前半は補題 2.5 から明らか。 A が正則であるならば、 L も正則であることに注意して、 L := L^{−1} とおけば A = LU. L は単位下三角であるから det L^{(k)} = 1 であるので、 det A^{(k)} = det U^{(k)} = ∏_{i=1}^{k} u_{ii}.

以上をまとめると、次の定理が得られる。

[figure]

定理 2.7 正則行列 A が LU 分解可能であるためには、 A のすべての主座小行列式が 0 でないことが必要十分である。そのとき Gauss の消去法の前進消去過程を施すことによって、 A の LU 分解が得られる。

[figure]

この定理の、簡単だが応用上重要な系として、次のものがある。

[figure]

系 2.8 A を n 次実対称行列とする。 A が正値であるための必要十分条件は、 A は Gauss の消去法の前進消去過程によって、単位下三角行列 L と上三角行列 U = (u_{ij}) の積として LU 分解でき、 u_{kk} > 0 (k = 1, …, n) が成り立つことである。

[figure]

証明よく知られているように、対称行列 A が正値であるためには

det A^{(k)} > 0 (k = 1, …, n)

が成り立つことが必要十分である。これを思い出せば明らか。

なお、正値対称行列の LU 分解としては、 Cholesky 分解が重要である。これについては、桂田 [1] などを見よ。

LU 分解の一意性についても片付けておこう。

[figure]

命題 2.9 正則行列を単位下三角行列と上三角行列の積に分解する仕方は一意的である。すなわち

A = L_{1}U_{1} = L_{2}U_{2} (L_{1} と L_{2} は単位下三角行列, U_{1} と U_{2} は上三角行列)

ならば、 L_{1} = L_{2} かつ U_{1} = U_{2} が成り立つ。

[figure]

証明 A が正則であることから U_{1} も正則である。したがって

L_{2}^{−1}L_{1} = U_{2}U_{1}^{−1}.

左辺は単位下三角行列の積であるから単位下三角行列であり、右辺は上三角行列の積であるから上三角行列である。ゆえにこの等式の両辺は、単位下三角行列かつ上三角行列であるので、実は単位行列 I に等しい。これから

L_{1} = L_{2}, U_{1} = U_{2}

が得られる。

やり残し「任意の正則行列 A に対して、ある置換行列 P が存在して、 P A は LU 分解できる」を証明したい。

2.8 プログラムの書き方

A = LU (L は下三角行列, U は上三角行列)

と LU 分解できるとき、 L と U を計算機で計算することを考える。

二つの行列 L と U を記憶するには、行列ちょうど 1 個分の n^{2} 個の成分を記憶できれば十分である。つまり

L = (ℓ_{ij}), U = (u_{ij})

とするとき、

i = j ⟹ ℓ_{ij} = 1,  i < j ⟹ ℓ_{ij} = 0,  i > j ⟹ u_{ij} = 0

が分かっているので、それらについては記憶しておく必要がない。上の例で

A = ( \begin{matrix}2 & 3 & −1 \\ 4 & 4 & −3 \\ −2 & 3 & −1 \end{matrix} )

から

L = ( \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ −1 & −3 & 1 \end{matrix} ) , U = ( \begin{matrix}2 & 3 & −1 \\ 0 & −2 & −1 \\ 0 & 0 & −5 \end{matrix} )

を得たいわけだが、 L の対角線より上 (含む対角線) と、 U の対角線から下 (対角線は含まない) は不要なので、結果を

( \begin{matrix}2 & 3 & −1 \\ 2 & −2 & −1 \\ −1 & −3 & −5 \end{matrix} )

という形に圧縮することができる。そこで実際のプログラムでは、次のようなコーディングをする。 C 言語風疑似プログラム

[figure]

    for (k = 1; k < n; k++) {
      /* (k,k) 成分で掃き出す */
      for (i = k + 1; k <= n; k++) {
        /* q_{ik} = a_{ik} / a_{kk} を求める */
        q = a[i][k] / a[k][k];

/* i 行から第 k 行の q 倍を引く (ただし計算は j (≧ k+1) 列のみ実行) */

        for (j = k + 1; j <= n; j++)
          a[i][j] -= q * a[k][j];

/* q_{ik} = L_{ik} を記憶しておく (記憶不要になったメモリーを利用) */

        a[i][k] = q;
      }
    }

[figure]

この事情は、係数行列が疎行列の場合も基本的には変らない。例えば三重対角行列を LU 分解し、それを用いて連立 1 次方程式を解くプログラムを末尾に引用するが、 LU 分解の部分は

三重対角行列の LU 分解

[figure]

  /* 前進消去 (forward elimination) */
  for (k = 1; k < n; k++) {
    al[k + 1] /= ad[k];             /* a[i][k] /= a[k][k] */
    ad[k + 1] -= au[k] * al[k + 1]; /* a[i][j] -= a[k][j] * q  (j は k+1 のみ) */
  }

[figure]

と非常にコンパクトに記述できている (A が三重対角行列の場合は、上三角因子 U の対角線より上の部分はもとの行列 A のそれと同じ (なぜそうなるか考えて確かめよ) という特徴があるため、 j に関するループが一つの代入文に置き換わっている。また下三角因子 L は、対角線と、対角線より 1 つ下のところにだけ非零成分が存在しうる。こういう事情のおかげで、プログラムは極限までコンパクトになっている)。

以下、三重対角行列 A の LU 分解が得られているとして、連立 1 次方程式 Ax = b をどうやって解くか述べる。すでに説明したように、 Ly = b, Ax = y という二つの方程式に分解する。 Ly = b を解く部分は、 ℓ_{ii} = 1 と、三重対角であることに注意すると

y_{1} = b_{1},

そして i = 2, 3, …, n の順に

y_{i} = b_{i} − ℓ_{i,i−1}y_{i−1}

とすればよい。最初に b[] に b が記憶されていれば

[figure]

    for (i = 2; i <= n; i++)
      b[i] -= al[i] * b[i-1];

[figure]

というコードで計算すると b[] に y が計算 (して記憶) される。

次に Ax = y を解く部分は

x_{n} = \frac{y_{n}}{u_{nn}}

としてから、 i = n − 1, n − 2, …, 1 の順に

x_{i} = \frac{(y_{i} − u_{i,i+1}x_{i+1})}{u_{ii}}

とすればよい。 b[] に y が記憶されていれば

[figure]

   b[n] /= ad[n];
   for (i = n-1; i >= 1; i--)
     b[i] = (b[i] - au[i] * b[i+1]) / ad[i];

[figure]

というコードで計算すると b[] に x が記憶 (して記憶) される。

3 重対角行列の LU 分解プログラム

/* trid-lu1.c -- 3 項方程式を Gauss の消去法で解く */
#include "trid-lu1.h"

/* 3 項方程式 (係数行列が三重対角である連立 1 次方程式のこと) Ax=b を解く

*

* 入力

* n: 未知数の個数

*     al,ad,au: 連立 1 次方程式の係数行列
*       (al: 対角線の下側 i.e. 下三角部分 (lower part)
*        ad: 対角線  i.e. 対角部分 (diagonal part)
*        au: 対角線の上側 i.e. 上三角部分 (upper part)

* つまり

*
*         ad[1] au[1]   0   .................. 0
*         al[2] ad[2] au[2]   0  ............. 0
*           0   al[3] ad[3] au[3]  0 ......... 0
*                         ....................
*                         al[n-1] ad[n-1] au[n-1]
*                           0     al[n]   ad[n]
*        al[i] = A_{i,i-1}, ad[i] = A_{i,i}, au[i] = A_{i,i+1},
*        al[1], au[n] は意味がない)
*

* b: 連立 1 次方程式の右辺の既知ベクトル

*        (添字は 1 から。 i.e. b[1],b[2],...,b[n] にデータが入っている。 )

* 出力

* al,ad,au: 入力した係数行列を LU 分解したもの

* b: 連立 1 次方程式の解

* 能書き

* 一度 call すると係数行列を LU 分解したものが返されるので、

* 以後は同じ係数行列に関する連立 1 次方程式を解くために、

*  関数 trisol1() が使える。

* 注意

* Gauss の消去法を用いているが、ピボットの選択等はしていな

* いので、ピボットの選択をしていないので、係数行列が正定値である

* などの適切な条件がない場合は結果が保証できない。

*/
void trid1(int n, double *al, double *ad, double *au, double *b)

{
trilu1(n,al,ad,au);
trisol1(n,al,ad,au,b);
}
/* 三重対角行列の LU 分解 (pivoting なし) */
void trilu1(int n, double *al, double *ad, double *au)
{
int i;
/* 前進消去 (forward elimination) */
for (i = 1; i < n; i++) {
  al[i + 1] /= ad[i];
  ad[i + 1] -= au[i] * al[i + 1];
}
}

/* LU 分解済みの三重対角行列を係数に持つ 3 項方程式を解く */

void trisol1(int n, double *al, double *ad, double *au, double *b)
{
int i;
/* 前進消去 (forward elimination) */
for (i = 1; i < n; i++) b[i + 1] -= b[i] * al[i + 1];
/* 後退代入 (backward substitution) */
b[n] /= ad[n];
for (i = n - 1; i >= 1; i--) b[i] = (b[i] - au[i] * b[i + 1]) / ad[i];
}

これを使って、

A = ( \begin{matrix}2 & 3 & 0 \\ 4 & 4 & −3 \\ 0 & 3 & −1 \end{matrix} )

が

A = LU, L = ( \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 0 & −\frac{3}{2} & 1 \end{matrix} ) , U = ( \begin{matrix}2 & 3 & 0 \\ 0 & −2 & −3 \\ 0 & 0 & −\frac{11}{2} \end{matrix} ) ,

と LU 分解されることを確かめてみよう。

[figure]

/*
* test-trid-lu1.c
*   gcc -o test-trid-lu1 test-trid-lu1.c trid-lu1.c
*/
#include <stdio.h>
#include "trid-lu1.h"
#define N 3
int main(void)
{
int n;
double al[N+1], ad[N+1], au[N+1], b[N+1];
n = 3;
/* 行列 A の設定 */
ad[1] = 2; au[1] = 3;
al[2] = 4; ad[2] = 4; au[2] = -3;
           al[3] = 3; ad[3] = -1;
/* A の LU 分解 */
trilu1(n, al, ad, au);

/* LU 分解の結果の表示 */

printf("%f %f %f\n", ad[1], au[1], 0.0);
printf("%f %f %f\n", al[2], ad[2], au[2]);
printf("%f %f %f\n", 0.0,   al[3], ad[3]);
/* x=(1,2,3) に対応する b=A x */
b[1] = 8; b[2] = 3; b[3] = 3;
/* A x = b を解く */
trisol1(n, al, ad, au, b);

/* 解を表示する */

printf("%f %f %f\n", b[1], b[2], b[3]);
return 0;
}

[figure]
[figure]

oyabun% gcc -o test-trid-lu1 test-trid-lu1.c trid-lu1.c
oyabun% ./test-trid-lu1
2.000000 3.000000 0.000000
2.000000 -2.000000 -3.000000
0.000000 -1.500000 -5.500000
1.000000 2.000000 3.000000
oyabun%

[figure]

3 行列法による安定性解析

行列のスペクトルを調べて安定性を解析する方法について、既に桂田 [3] の中でまとめてあるが、書き直す必要を感じていて、その下書きとして以下の説明を書く。 [3] と適当にマージして完成させる予定。

以下この節では、任意の n ∈ N に対して、 n 次の単位行列を I_{n} で、 (n ≥ 2 のとき) 対角線の両隣の成分のみが 1 で他の成分が 0 である n 次正方行列を J_{n} で表す:

I_{n} := ( \begin{matrix}1 & 0 & & 0 \\ 0 & 1 & ... \\ & ... & ... & 0 \\ 0 & & 0 & 1 \end{matrix} ) , J_{n} := ( \begin{matrix}0 & 1 \\ 1 & 0 & 1 \\ & ... & ... & ... \\ & & 1 & 0 & 1 \\ & & & 1 & 0 \end{matrix} ) .

3.1 差分解は「等比ベクトル列」

同次 Dirichlet 境界条件下の熱方程式の初期値境界値問題の θ 法による差分方程式は、

U⃗^{n} := (U_{1}^{n}, ⋯, U_{N−1}^{n})^{T} ,
A = A_{N−1} := (1 + 2θλ)I_{N−1} − θλJ_{N−1},
B = B_{N−1} := [1 − 2(1 − θ)λ] I_{N−1} + (1 − θ)λJ_{N−1}

とおくと、

(5)  AU⃗^{n+1} = BU⃗^{n}

と書ける。そこで

R = R_{N−1} := (A_{N−1})^{−1} B_{N−1}

とおくと、

(6)  U⃗^{n+1} = RU⃗^{n}.

ゆえに

(7)  U⃗^{n} = R^{n}U⃗^{0}.

(6) や (7) は、実際の計算には役に立たない (A, B は疎行列であるので、連立 1 次方程式 (5) の解として U⃗^{n+1} を計算するのは非常に効率的にできるのに対して、 R を計算して (これは密行列) それを用いて U⃗^{n+1} を計算するのは非常に非効率的になってしまう)。しかし、安定性の解析のような理論的な考察には大いに役立つ。

3.2 Jordan 標準形を用いた「等比マトリックス列」の解析

(行列の列を「行列列」というのは何か気が引けるので、マトリックス列ということにする。 )

[figure]

命題 3.1 (等比数列行列版の解析)A ∈ M(n; C) として、 A のスペクトル半径を r(A) と書くとき、以下の (1)–(4) が成り立つ。

(1) lim_{m→∞} A^{m} = O であるための必要十分条件は r(A) < 1.

(2) lim_{m→∞} A^{m} が存在するための必要十分条件は、次の (a), (b) が成り立つこと。

(a) r(A) ≤ 1.

(b) |λ| = 1 となる固有値 λ があるとき、実は λ = 1 で、それは A の最小多項式の単根である。

(3) {A^{m}}_{m∈N} が有界であるための必要十分条件は、次の (a), (b) が成り立つこと。

(a) r(A) ≤ 1.

(b) |λ| = 1 となる固有値 λ は A の最小多項式の単根である。

[figure]

証明与えられた A に対して、適当な P ∈ GL(n; C) が存在して、

P^{−1}AP = J

となる。ただし J は Jordan 行列である:

J = ( \begin{aligned}λ_{1} & 1 \\ & λ_{1} & 1 \\ & & ... & ... \\ & & & λ_{1} & 1 \\ & & & & λ_{1} \end{aligned} ; \begin{aligned}λ_{2} & 1 \\ & λ_{2} & 1 \\ & & ... & ... \\ & & & λ_{2} & 1 \\ & & & & λ_{2} \end{aligned} ; ... ; \begin{aligned}λ_{r} & 1 \\ & λ_{r} & 1 \\ & & ... & ... \\ & & & λ_{r} & 1 \\ & & & & λ_{r} \end{aligned} )

= J(λ_{1}, n_{1}) ⊕ J(λ_{2}, n_{2}) ⊕ ⋯ ⊕ J(λ_{r}, n_{r}).

ただし、 λ ∈ C, n ∈ N に対して、 J(λ, n) を次式で定義する。

J(λ, n) := ( \begin{matrix}λ & 1 & & 0 \\ & ... & ... \\ & & λ & 1 \\ 0 & & & λ \end{matrix} ) ∈ M(n; C).

任意の自然数 m に対して A^{m} = P J^{m}P^{−1} であるから、

lim_{m→∞} A^{m} = O ⟺ lim_{m→∞} J^{m} = O,

lim_{m→∞} A^{m} が存在 ⟺ lim_{m→∞} J^{m} が存在,

{A^{m}}_{m} が有界 ⟺ {J^{m}}_{m} が有界

が分かる。さらに

J^{m} = J(λ_{1}, n_{1})^{m} ⊕ J(λ_{2}, n_{2})^{m} ⊕ ⋯ ⊕ J(λ_{r}, n_{r})^{m}

であるから、 J^{m} を調べるには、 J(λ, n)^{m} を調べればよい。 m に関する数学的帰納法により、容易に

J(λ, n)^{m} = ( \begin{matrix}λ^{m} & \binom{m}{1}λ^{m−1} & \binom{m}{2}λ^{m−2} & ⋯ & \binom{m}{n−1}λ^{m−n+1} \\ & λ^{m} & \binom{m}{1}λ^{m−1} & ... & ⋮ \\ & & λ^{m} & ... & \binom{m}{2}λ^{m−2} \\ & & & ... & \binom{m}{1}λ^{m−1} \\ 0 & & & & λ^{m} \end{matrix} )

であることが分かる。これから容易に^{10}

|λ| < 1 ⟺ lim_{m→∞} J(λ, n)^{m} = O.

ゆえに

r(A) < 1 ⟺ lim_{m→∞} A^{m} = O.

一方、明らかに

|λ| > 1 ⟹ lim_{m→∞} ∥J(λ, n)^{m}∥ = ∞ (特に J(λ, n)^{m} は収束しない)

であるから、

r(A) > 1 ⟹ lim_{m→∞} ∥A^{m}∥ = ∞ (特に A^{m} は収束しない).

|λ| = 1 の場合を考えよう。 J(λ, n)^{m} の対角成分 λ^{m} は有界である。 J(λ, n)^{m} そのものが有界であるためには、 n = 1 が必要十分である。

また m → ∞ のとき J(λ, n)^{m} が収束するためには、まず対角成分 λ^{m} が収束する必要がある。そのための条件は λ = 1. そして J(λ, n)^{m} は有界でなければならないから、上に述べたように n = 1 が必要である。逆に λ = 1 かつ n = 1 であるとき、

J(λ, n) = J(1, 1) = 1, J(λ, n)^{m} = 1^{m} = 1.

これは定数なので、もちろん m → ∞ のとき収束する。

最小多項式の単根であるとは、対応する Jordan 細胞のサイズが 1 (つまりただの複素数) であることに他ならない。

3.3 J のスペクトル解析

行列 J = J_{N−1} の固有値、固有ベクトルは次のように完全に求まる。

[figure]

命題 3.2 任意の N ≥ 3 に対して、

2 cos nπh, v⃗_{n} = (sin nπx_{1}, ⋯, sin nπx_{N−1})^{T} (n = 1, 2, ⋯, N − 1)

は J_{N−1} ∈ M_{N−1}(R) の固有値、固有ベクトルである。ただし

h = \frac{1}{N} , x_{j} = jh (j = 0, 1, ⋯, N). [figure]

この命題の証明を直接行うことも出来るが^{11}、それ自身明快な意味を持つ次の命題を足掛かりに証明することにする。

10念のため書いておくと、 λ^{m} が 0 に収束するために |λ| < 1 が必要で、逆にこのとき任意の k に対して \binom{m}{k}λ^{m−k} → 0 (m → ∞) が分かるので、 J^{m} → O となるための十分条件である。

11 さわりの部分は sin \frac{(j − 1)nπ}{N} + sin \frac{(j + 1)nπ}{N} = 2 cos \frac{nπ}{N} sin \frac{jnπ}{N}.

[figure]

補題 3.3 N ≥ 3 なる任意の N ∈ N に対して、

(2 sin \frac{nπh}{2})^{2} , v⃗_{n} = (sin nπx_{1}, ⋯, sin nπx_{N−1})^{T} (n = 1, 2, ⋯, N − 1)

は

L_{N−1} := 2I_{N−1} − J_{N−1} = ( \begin{matrix}2 & −1 & & & 0 \\ −1 & 2 & −1 \\ & ... & ... & ... \\ & & −1 & 2 & −1 \\ 0 & & & −1 & 2 \end{matrix} )

の固有値、固有ベクトルである。ただし

h = \frac{1}{N} , x_{j} = jh (j = 0, 1, ⋯, N). [figure]

注意 3.4 1 次元 Laplacian の固有値問題

(EP)  \{ \begin{matrix}−u′′ = λu in (0, 1), \\ u(0) = u(1) = 0 \end{matrix}

を差分法で「すなおに」離散化すると

\frac{1}{h^{2}} L_{N−1}U⃗ = λU⃗

という行列の固有値問題が現れる。 (EP) の固有値・固有関数は

n^{2}π^{2}, sin nπx (n = 1, 2, ⋯)

で与えられることを思い出すと、上の補題は

(EP) を差分法ですなおに離散化して得られた固有値問題の固有関数は、もとと同じく sin nπx である

ということを主張していることになる。このことを理解すると、上の補題の主張の後半 (v⃗_{n} が固有ベクトルになること) は覚えやすいであろう。なお、番号 n が小さいときに、対応する固有値も良い近似になっている。このこともチェックしておこう。

証明 n ∈ N に対して、

v_{n}(x) := sin nπx,
v_{nj} := v_{n}(x_{j}) = sin nπx_{j} = sin nπjh = Im (e^{inπjh}) (j = 0, 1, ⋯, N),
v⃗_{n} := (v_{n1}, v_{n2}, ⋯, v_{n,N−1})^{T}

とおく。 v_{n0} = v_{n,N} = 0 に注意しておく。任意の j ∈ {1, 2, …, N − 1} について

L_{N−1}v⃗_{n} の第 j 成分 = 2v_{nj} − v_{n,j−1} − v_{n,j+1}
= Im (2e^{inπjh} − e^{inπ(j−1)h} − e^{inπ(j+1)h})
= Im [e^{inπjh} (2 − e^{−inπh} − e^{inπh})]
= Im [e^{inπjh} · 2 (1 − cos nπh)] .

ここで

2(1 − cos nπh) = 2 · 2 sin^{2} \frac{nπh}{2} = 4 sin^{2} \frac{nπh}{2} =: λ_{n}

は実数であることに注意すると

L_{N−1}v⃗_{n} の第 j 成分= λ_{n} Im e^{inπjh} = λ_{n}v_{nj}.

すなわち

L_{N−1}v⃗_{n} = λ_{n}v⃗_{n}.

ところで

0 < \frac{1 · πh}{2} < \frac{2πh}{2} < ⋯ < \frac{(N − 1)πh}{2} < \frac{π}{2}

であるから

0 < λ_{1} < λ_{2} < ⋯ < λ_{N−1} < 4.

特に {λ_{n}}_{1≤n≤N−1} は互いに相異なる。個数を考えると、これが L_{N−1} の固有値全体となる。

命題 3.2 の証明 J_{N−1} = 2I_{N−1} − L_{N−1} であるから、

μ_{n} := 2 − λ_{n} = 2 − 4 sin^{2} \frac{nπh}{2} = 2 (1 − 2 sin^{2} \frac{nπh}{2}) = 2 cos nπh

とおくと、

J_{N−1}v⃗_{n} = (2I_{N−1} − L_{N−1})v⃗_{n} = 2I_{N−1}v⃗_{n} − L_{N−1}v⃗_{n} = 2v⃗_{n} − λ_{n}v⃗_{n} = (2 − λ_{n})v⃗_{n}
= μ_{n}v⃗_{n} (n = 1, 2, ⋯, N − 1).

3.4 R のスペクトル解析

f を

(8)  f(μ) := \frac{[1 − 2(1 − θ)λ] + (1 − θ)λμ}{(1 + 2θλ) − θλμ} = 1 + \frac{λ(μ − 2)}{1 + θλ(2 − μ)}

で定めると

R_{N−1} = f(J_{N−1})

となる。 J_{N−1} の固有値、固有ベクトルは前項で調べたようにμ_{n} = 2 cos nπh, v⃗_{n} (n = 1, 2, ⋯, N− 1) で与えられるから、 R_{N−1} の固有値、固有ベクトルは f(μ_{n}), v⃗_{n} (n = 1, 2, ⋯, N − 1) で与えられることが分かる (Frobenius の定理の一般化、桂田 [3] の第 4 章にある「有理式版フロベニウスの定理」 )。

[figure]

補題 3.5 0 ≤ θ ≤ 1, λ > 0 とする。 f を (8) で定義するとき、次の (1)–(3) が成り立つ。

(1) f は [−2, 2] で狭義単調増加である。

(2) f(2) = 1.

(3) f(−2) ≥ −1 ⇔ θ ≥ \frac{1}{2} or (θ < \frac{1}{2} and λ ≤ \frac{1}{2(1 − 2θ)}).
[figure]

証明 (1) まず μ ≤ 2 とするとき、 1 + θλ(2 − μ) ≥ 1 であることに注意する。

f′(μ) = \frac{[1 + θλ(2 − μ)] · λ − (−θλ) · λ(μ − 2)}{[1 + θλ(2 − μ)]^{2}} = \frac{λ}{[1 + θλ(2 − μ)]^{2}} > 0

であるから、 f は [−2, 2] で狭義単調増加である。 (2) は単純な計算で分かる。 (3) については、

f(−2) = \frac{1 + 4θλ − 4λ}{1 + 4θλ}

であるから、

f(−2) − (−1) = f(−2) + 1 = \frac{1 + 4θλ − 4λ + 1 + 4θλ}{1 + 4θλ} = \frac{2 [1 + 2λ(2θ − 1)]}{1 + 4θλ}.

この分母は正であるから、

f(−2) ≥ −1 ⇔ 1+2λ(2θ −1) ≥ 0 ⇔ θ ≥ \frac{1}{2} or (θ < \frac{1}{2} and λ ≤ \frac{1}{2(1 − 2θ)}) .

さて、 μ_{n} = 2 cos nπh = 2 cos \frac{nπ}{N} より

2 > μ_{1} > μ_{2} > ⋯ > μ_{N−1} > −2

であるから、上の補題の (1), (2) を用いると、

1 = f(2) = f(μ_{1}) > f(μ_{2}) > ⋯ > f(μ_{N−1}) > f(−2).

補題の (3) より、

θ ≥ \frac{1}{2} or (θ < \frac{1}{2} and λ ≤ \frac{1}{2(1 − 2θ)}) ⟹ |f(μ_{n})| < 1 (n = 1, 2, ⋯, N − 1).

さらに、もしこの仮定が成り立たない、すなわち^{12}

θ < \frac{1}{2} and λ > \frac{1}{2(1 − 2θ)}

ならば、 f(−2) < −1 となるので、

μ_{N−1} = −μ_{1} = −2 cos \frac{π}{N} ↓ −2 (N → ∞)

に注意すると、十分大きな N に対して

f(μ_{N−1}) < −1 ゆえに |f(μ_{N−1})| > 1

となることが分かる。まとめると次の定理を得る。

[figure]

定理 3.6 0 ≤ θ ≤ 1, λ > 0 とするとき、次の (1), (2), (3) が成り立つ。

(1) θ, λ が条件

(9)  θ ≥ \frac{1}{2} or (θ < \frac{1}{2} and λ ≤ \frac{1}{2(1 − 2θ)})

を満たすならば、 N ≥ 3 なる任意の N ∈ N に対して、 R_{N−1} のスペクトル半径は 1 より小さい。

(2) θ, λ が条件 (9) を満足しないとき、十分大きい N を取ると、 R_{N−1} のスペクトル半径は 1 より大きい。

(3) 差分スキームが安定であるためには、 (9) が必要十分である。

[figure]

12ド・モルガン則など、論理の法則を使って計算すると、機械的に計算できる。

3.5 差分解の漸近挙動の解析

前項の議論から、 R_{N−1} のスペクトル半径が 1 より大きいとき、絶対値最大の固有値は f(μ_{N−1}) であることがわかる。このとき、大抵の初期値 U⃗^{0} に対して、 n が大きいとき、 U⃗^{n} は v⃗_{N−1} 成分がメジャーになる (行列の固有値・固有ベクトルを求めるアルゴリズムである冪乗法の原理と同じ)。 v_{N−1}(x) = sin(N − 1)πx は [0, 1] に山と谷が合わせて N − 1 個ある関数であるから、 v⃗_{N−1} = (v_{N−1}(x_{1}), ⋯, v_{N−1}(x_{N−1}))^{T} はジグザグしているグラフを持つことが想像できるであろう (きちんとしたステートメントにして証明できるか?)。それゆえ、不安定なスキームによる差分解のグラフがジグザグしていることは納得できるであろう。

ちゃんと命題の形で書くこと。ジグザグしていることも、良く出て来ることなので、もっときちんと書いておきたい。

問題安定の場合、 f(μ_{1}) と f(μ_{N−1}) のいずれの絶対値が大きいか？つねに |f(μ_{1})| > |f(μ_{N−1})| か、それとも |f(μ_{1})| < |f(μ_{N−1})| となることがあるか？ — 状況が良い場合は |f(μ_{1})| > |f(μ_{N−1})| となるが、 |f(μ_{1})| < |f(μ_{N−1})| となることもありそうである。差分解のグラフは n が大きいところで、前者の場合は sin カーブの半周期分だが、後者の場合はジグザグしていることになる。

3.6 安定性の条件を満たす場合の差分解の収束証明

多分、最大値原理を満たす場合の差分解の収束証明の真似をして、差分解の厳密解への収束定理が得られると思われる。卒業研究レポートの課題に適当？

この文書とはかなり書き方が違うが、田端 [7] に載っている。

A メモ (未整理)

A.1 「下三角行列の逆行列は下三角行列」

A = (a_{ij}) を下三角行列とする。下三角行列 B = (b_{ij}) で AB = I を満たすものが存在することを示す (そうすれば一般に得られている逆行列の一意性から OK)。

∑_{k=1}^{n} a_{ik}b_{kj} = δ_{ij} (1 ≤ i ≤ n, 1 ≤ j ≤ n).

これは i < j ならば両辺 0 でつねに成立している。 i ≤ j とするとき、

∑_{k=j}^{i} a_{ik}b_{kj} = δ_{ij}.

これは

∑_{k=j}^{i−1} a_{ik}b_{kj} + a_{ii}b_{ij} = δ_{ij}.

と書き直せる。 ∑ の部分は既知、 a_{ii} ≠ 0 であることに注意すると、 b_{ij} が求まることが分かる。

A.2 Lagrange の方法による 2 次形式の標準形への変換

次の 2 つは齋藤 [5] で Lagrange の方法として紹介されている例。

例 A.1 A = (\begin{matrix}1 & 0 & 2 \\ 0 & 1 & 2 \\ 2 & 2 & −1 \end{matrix}) を係数行列とする 2 次形式 (Ax, x) = x^{2} + y^{2} − z^{2} + 4xz + 4yz は、

いわゆる平方完成によって、

(Ax, x) = (x + 2z)^{2} + y^{2} − 5z^{2} + 4yz = (x + 2z)^{2} + (y + 2z)^{2} − (3z)^{2}.

ゆえに

x′ := x + 2z, y′ := y + 2z, z′ := 3z

とおくと、

(Ax, x) = x^{′2} + y^{′2} − z^{′2}

と標準形に変換される。変換行列 P は、

x = P x′, x := (\begin{matrix}x \\ y \\ z \end{matrix}) , x′ := (\begin{matrix}x′ \\ y′ \\ z′ \end{matrix})

を満たす行列であるから

P^{−1} = (\begin{matrix}1 & 0 & 2 \\ 0 & 1 & 2 \\ 0 & 0 & 3 \end{matrix}) .

逆変換は

z = \frac{z′}{3} , y = y′ − 2z = y′ − \frac{2}{3}z′, x = x′ − 2z = x′ − \frac{2}{3}z′

であるから、

P = (\begin{matrix}1 & 0 & −\frac{2}{3} \\ 0 & 1 & −\frac{2}{3} \\ 0 & 0 & \frac{1}{3} \end{matrix} ) .

実際

P^{T} AP = (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & −1 \end{matrix}) .

次に LU 分解を使って見直す。容易に

A = LU, L = (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 2 & 2 & 1 \end{matrix}) , U = (\begin{matrix}1 & 0 & 2 \\ 0 & 1 & 2 \\ 0 & 0 & −9 \end{matrix})

と分かる。 U を見るだけで

(Ax, x) = (x + z)^{2} + (y + 2z)^{2} − 9z^{2}

と分かる (L は不要らしい)。例 A.2 A = (\begin{matrix}0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{matrix}) を係数行列とする 2 次形式

(Ax, x) = 2xy + 2yz.
X = x + y, Y = x − y

として、

(Ax, x) = \frac{1}{2}X^{2} − \frac{1}{2}Y^{2} + (X − Y)z

= \frac{1}{2}(X + z)^{2} − \frac{1}{2}Y^{2} − Y z − \frac{1}{2}z^{2}

= \frac{1}{2}(X + z)^{2} − \frac{1}{2}(Y + z)^{2}

= \frac{1}{2}(x + y + z)^{2} − \frac{1}{2}(x − y + z)^{2}.

符号は (1, 1).

(Ax, x) = (\frac{x + y + z}{\sqrt{2}} )^{2} − (\frac{x − y − z}{\sqrt{2}} )^{2}

であるから、

P^{−1} = \frac{1}{\sqrt{2}} (\begin{matrix}1 & 1 & 1 \\ 1 & −1 & 1 \\ 0 & 0 & 1 \end{matrix}) .

例 A.3 A = (\begin{matrix}3 & 0 & 0 \\ 0 & −1 & 3 \\ 0 & 3 & −1 \end{matrix}) とすると、

U = (\begin{matrix}3 & 0 & 0 \\ 0 & −1 & 3 \\ 0 & 0 & 8 \end{matrix}) = (\begin{matrix}3 & 0 & 0 \\ 0 & −1 & 0 \\ 0 & 0 & 8 \end{matrix}) (\begin{matrix}1 & 0 & 0 \\ 0 & 1 & −3 \\ 0 & 0 & 8 \end{matrix})

と分かるので、 (Ax, x) = 3x^{2} − y^{2} − z^{2} + 6yz = 3x^{2} − (y − 3z)^{2} + 8z^{2} = (\sqrt{3}x)^{2} − (y − 3z)^{2} + (\sqrt{2}z)^{2} .

例 A.4 A = (\begin{matrix}1 & \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & 1 & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{2} & 1 \end{matrix}) とするとき、

U = (\begin{matrix}1 & \frac{1}{2} & \frac{1}{2} \\ 0 & \frac{3}{4} & \frac{1}{4} \\ 0 & 0 & \frac{2}{3} \end{matrix} ) = (\begin{matrix}1 & 0 & 0 \\ 0 & \frac{3}{4} & 0 \\ 0 & 0 & \frac{2}{3} \end{matrix}) (\begin{matrix}1 & \frac{1}{2} & \frac{1}{2} \\ 0 & 1 & \frac{1}{3} \\ 0 & 0 & 1 \end{matrix}) .

これから

(Ax, x) = x^{2} + y^{2} + z^{2} + xy + yz + zx

= (x + \frac{1}{2}y + \frac{1}{2}z)^{2} + \frac{3}{4} (y + \frac{1}{3}z)^{2} + \frac{2}{3}z^{2}.

例 A.5 A = (\begin{matrix}1 & 1 & −1 & 1 \\ 1 & 4 & 2 & 1 \\ −1 & 2 & 4 & 0 \\ 1 & 1 & 0 & −1 \end{matrix}) は

A → A′ = (\begin{matrix}1 & 1 & −1 & 1 \\ 0 & 3 & 3 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & −2 \end{matrix})

まで掃き出しが出来る。これは

(Ax, x) = (x + y − z + u)^{2} + 3(y + z)^{2} + 2zu − 2u^{2}

ということを表している。 z と u を取り替えて (3 列と 4 列を入れ替え、 3 行と 4 行を入れ替える)、掃き出しを続行すると

A′′ = (\begin{matrix}1 & 1 & 1 & −1 \\ 0 & 3 & 0 & 3 \\ 0 & 0 & −2 & 1 \\ 0 & 0 & 1 & 0 \end{matrix}) → A′′′ = (\begin{matrix}1 & 1 & 1 & −1 \\ 0 & 3 & 0 & 3 \\ 0 & 0 & −2 & 1 \\ 0 & 0 & 0 & \frac{1}{2} \end{matrix}) = (\begin{matrix}1 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & −2 & 0 \\ 0 & 0 & 0 & \frac{1}{2} \end{matrix}) (\begin{matrix}1 & 1 & 1 & −1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & −\frac{1}{2} \\ 0 & 0 & 0 & 1 \end{matrix})

ゆえに

(Ax, x) = (x + y + u − z)^{2} + 3(y + z)^{2} − 2 (u − \frac{1}{2}z)^{2} + \frac{1}{2}z^{2}.

例 A.6 A = (\begin{matrix}0 & 1 & 1 & 1 \\ 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \end{matrix})

B Changelog

• 2012/7/17 下三角行列の定義を間違えていた！まだこんなのが残っていたとは。

• 2012/7/22 2 次形式の平方完成の話を書き加える。まだ書きかけなので本文には入れない。

参考文献

[1] 桂田祐史. Cholesky 分解ノート. http://nalab.mind.meiji.ac.jp/~mk/labo/text/
   cholesky.pdf.
[2] 桂田祐史. 発展系の数値解析. http://nalab.mind.meiji.ac.jp/~mk/labo/text/
   heat-fdm-0.pdf, 1997 年〜.

[3] 桂田祐史. 熱方程式に対する差分法 I — 区間における熱方程式 —. http://nalab.mind.

   meiji.ac.jp/~mk/labo/text/heat-fdm-1.pdf, 1998 年〜.
[4] 桂田祐史. 連立 1 次方程式 i. http://nalab.mind.meiji.ac.jp/~mk/labo/text/
   linear-eq-1.pdf, 2002 年〜.

[5] 齋藤^{さいとう}正彦^{まさひこ}. 線型代数入門. 東京大学出版会, 1966.

[6] 佐武一郎. 線型代数学. 裳華房, 1958, 1974. もともと『行列と行列式』という書名であったのを、テンソル代数の章を加筆した機会に改題した。

[7] 田端正久. 偏微分方程式の数値解法. 岩波書店, 2010. もともとは岩波講座応用数学の「微分方程式の数値解法 II」 (1994) であった。