数値計算レポート

{G99P057-3 \atop 斎藤卓也}

出題日： 2001 年 6 月 1 日
提出期限： 2001 年 6 月 15 日
提出日： 2001 年 6 月 15 日

1 問題

1. 教科書「数値計算」 P.31 の問題３を解け

2. 性質３と５と６を証明せよ（性質１，２，４は仮定してよい）。

1.1 教科書「数値計算」 P.31 の問題３

1. ガウスの消去法のプログラムを作り、連立一次方程式の近似解を求めてみよ。最初はピボット選択なしで、次に部分ピボット選択をするプログラムを作れ。

2. LU 分解法のプログラムを作り、連立一次方程式の近似解を求めてみよ。

3. LU 分解法を利用して逆行列を求めるプログラムを作れ。

4. 行列

A = [ \begin{matrix}1 & 2 & 3 \\ 2 & 8 & 11 \\ 3 & 22 & 35 \end{matrix} ] の LU 分解を求めよ。

5. ４の結果を用いて

b = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ]

b = [ \begin{matrix}1 \\ 2 \\ 3 \end{matrix} ]

b = [ \begin{matrix}−1 \\ 0 \\ 1 \end{matrix} ]

として Ax = b

を解け。

1.2 性質１

n × n 行列 A が ||A|| < 1 をみたすとする。このとき、 I を単位行列とすると

I + A

は正則行列となる。

1.3 性質２

A を n × n 行列とする。 ||A|| が誘導されたノルムとすると、元のベクトルのノルム || · || に対して

||Ax|| ≤ ||A|| · ||x||

1.4 性質３

A を n × n 行列、 ||A|| を誘導されたノルムとする。 R を n × n 行列、 I を n × n の単位行列とする。

||RA − I|| < 1

ならば、 A は正則であり、

||A^{−1} ≤ \frac{||R||}{1 − ||RA − I||}

が成り立つ。

1.5 性質４

A を n × n 行列、 ||A|| を誘導されたノルムとする。このとき、 λ を A の任意の固有値とすると、

|λ| ≤ ||A||

が成り立つ。

1.6 性質５

（ゲルシュゴリンの定理） A : n × n の（複素）行列とする。このとき、行列 A の固有値は

κ = {λ ∈ C||λ − a_{ii}| ≤ ∑_{k=1(k≠i}^{n} |a_{ik}|}

の中に存在する。

1.7 性質６

u = 2^{−53}

F で normalized 及び denormalized 及び 0 の double の集合を考える。

x, y ∈ F
∘ ∈ +, −, ∗, /

x ⊙ y = (x · y)(1 + δ) |δ| ≤ u

となる。これを利用すると

x, y· ∈ F^{n}

とすると、 s = x · y （実数としての計算）に対して、

near();
t=0;
for (i=1; i\le n; i++)
  t += x_i * y_i;

として計算した t との間に

|s − t| ≤ γ_{n} ∑_{i=1}^{n} |x_{i}y_{i}|, γ_{n} = \frac{nu}{1 − nu}

の関係がある。

2 教科書の問題３の解答

私の場合、柏木先生の課題で「行列ライブラリ」を作成したので、そのライブラリにガウスの法則（ピボット選択あり／なし）をする関数、 LU 分解をする関数、 LU 分解の結果を使い連立一次方程式の近似解を求める関数、 LU 分解法を利用して逆行列を求めるための関数をそのライブラリに追加した。

まず、ガウスの消去法についてのプログラムだが、「行列ライブラリ」にガウスの消去法で連立方程式を解く関数 Gauss() を追加した。

この行列ライブラリには、行列用のメモリ領域を確保する関数や、それを破棄する関数等が全て揃っているので、それらを利用すれば、簡単に Gauss() 関数を使って、連立方程式を解くことが出来る。

また、行列ライブラリファイルは、 lu.h と lu.c の２つに分かれている。 lu.h をインクルードすると、行列演算ライブラリ関数が使えるようになる。

用意されている関数は以下のものである。

• void Gauss(MATRIX* a, double* b, int pivot)

ax = b の形の連立方程式を解き、答えは b に入る。 pivot=1 だと、部分ピボット選択あり、 0 だとピボット選択なしで、ガウスの消去法を使い連立方程式を解く。

• void lu_decompose(MATRIX* u, MATRIX* l)

u で与えられた行列を u と l に LU 分解して返す。 l には u と同じサイズの行列を作成して渡してやる。

• void set_b(double* b, int size, int n)

b に、 n で指定した所だけ１をいれ、後は全部０を入れる関数。 LU 分解を利用して逆行列を求めるときのベクトルを作るときに使う。

• void set_inv(MATRIX* inv, double* b, int n)

n × n 行列 inv の n 列に、 n 次元ベクトル b の内容をコピーする関数。 LU 分解を利用して逆行列を作成するときに使用する。

• void lu_solve(MATRIX* u, MATRIX* l, double* b)

ＬＵ分解の結果を使って方程式を解く関数。 u と l に LU 分解により得られた行列を与え、また、 Ax = b の b に相当するベクトルを b として呼び出すと、 b に方程式を解いた結果を返す。

• int alloc_matrix(MATRIX* m, int x_width, int y_width)

MATRIX* m に、行 (x_width)、列 (y_width) で指定されただけメモリを確保する。

また、 m->x_width, m->y_width に引数で指定された値を設定する。

戻り値：メモリ確保に失敗したら０を、成功したら１を返す。

• void free_matrix(MATRIX* m)

MATRIX* の領域を開放する関数

• int add_matrix(MATRIX* a, MATRIX* b, MATRIX* c)

行列の足し算

c = a + b

戻り値：成功すれば１、失敗すれば０を返す

• int sub_matrix(MATRIX* a, MATRIX* b, MATRIX* c)

行列の引き算

c = a − b

戻り値：成功すれば１、失敗すれば０を返す

• int mul_matrix(MATRIX* a, MATRIX* b, MATRIX* c)

行列の掛け算

c = ab

戻り値：成功すれば１、失敗すれば０を返す

• int t_matrix(MATRIX* a, MATRIX* t)

転置行列を作る関数 t に a の転置行列を入れる戻り値：問題なければ１、問題があり転置行列が作れなければ０

• int naiseki_matrix(MATRIX* a, MATRIX* b, double* c)

行列の内積を求める関数

c =< a|b >

戻り値：問題が無ければ１。引数の *c に内積を入れて戻る問題があり内積が求められないときは０

• void clear_matrix(MATRIX* m)

行列の要素を全部０にする関数

• void print_matrix(MATRIX* m)

行列を表示する

関数の使い方は、以上で述べたとおりである。また、行列本体は、行列の実体が存在する場所を指すポインタ、行数、列数を、各行列で１つにまとめるため、構造体でまとめてある。また、その宣言は、 MATRIX という型が用意されるようにしてあるので

MATRIX* m1;

のようにすればよい。これで、この行列演算ライブラリで使用する行列型 MATRIX のポインタ m1 が宣言される。

MATRIX は、ヘッダファイルで以下のように定義されている。

typedef struct {
    int retu_width;
    int gyou_width;
    double** matrix;
} MATRIX;

2.0.1 問題３−１

プログラムの実行結果は以下のとおりである。これは、問題 3-4 と問題 3-5(1) の行列とベクトルを与えて解かせてみた例である。実際の解答は

[ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ −\frac{37}{24}_{\frac{5}{6}}^{\frac{19}{12}} ]

であるから、きちんと近似値が求められていることが確認できる。

A=
1.000000  2.000000  3.000000
2.000000  8.000000  11.000000
3.000000  22.000000  35.000000
b=
1.000000
0.000000
0.000000

答え

1.583333
-1.541667
0.833333
Press any key to continue

2.0.2 問題３−２

ここでは、以下の方程式を与え、ＬＵ分解によってこの連立一次方程式の近似解を求めてみることにした。

\{ \begin{matrix}2x_{1} + 6x_{2} + 4x_{3} = 8 \\ 3x_{1} + 7x_{2} + 2x_{3} = 6 \\ x_{1} + 2x_{2} + 3x_{3} = 7 \end{matrix}

実行結果は以下の通りである。

A=
2.000000  6.000000  4.000000
3.000000  7.000000  2.000000
1.000000  2.000000  3.000000
b=
8.000000
6.000000
7.000000
L=
1.000000  0.000000  0.000000
1.500000  1.000000  0.000000
0.500000  0.500000  1.000000
U=
2.000000  6.000000  4.000000
0.000000  -2.000000  -4.000000
0.000000  0.000000  3.000000

答え=

3.000000
-1.000000
2.000000
Press any key to continue

この結果より、以下のようにきちんと L、 U 分解が行われており、また、正しく解も得られている。

L = [ \begin{matrix}1 & 0 & 0 \\ 1.5 & 1 & 0 \\ 0.5 & 0.5 & 1 \end{matrix} ]

U = [ \begin{matrix}2 & 6 & 4 \\ 0 & −2 & −4 \\ 0 & 0 & 3 \end{matrix} ]

[ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ \begin{matrix}3 \\ −1 \\ 2 \end{matrix} ]

9 2.0.3 問題３−３

ここでは、以下の行列を与え、ＬＵ分解によってこの行列の逆行列を求めてみることにした。

[ \begin{matrix}1 & 2 & −4 \\ 2 & 1 & −4 \\ 1 & −2 & 1 \end{matrix} ]

これの逆行列は

[ \begin{matrix}−7 & 6 & −4 \\ −6 & 5 & −4 \\ −5 & 4 & −3 \end{matrix} ]

である。実際にプログラムで逆行列を求めてみると、以下のようになり、これは上の逆行列と一致している。

A=
1.000000  2.000000  -4.000000
2.000000  1.000000  -4.000000
1.000000  -2.000000  1.000000
U=
1.000000  2.000000  -4.000000
0.000000  -3.000000  4.000000
0.000000  0.000000  -0.333333
L=
1.000000  0.000000  0.000000
2.000000  1.000000  0.000000
1.000000  1.333333  1.000000

逆行列=

-7.000000  6.000000  -4.000000
-6.000000  5.000000  -4.000000
-5.000000  4.000000  -3.000000
Press any key to continue

2.1 問題３−４

A = [ \begin{matrix}1 & 2 & 3 \\ 2 & 8 & 11 \\ 3 & 22 & 35 \end{matrix} ]

これで、２行＝２行−２×１行、３行＝３行−３×１行として

[ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 16 & 26 \end{matrix} ]

となる。次に、３行＝３行−４×２行より

U = [ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{matrix} ]

また、以上掛けた数より

L = [ \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & 1 \end{matrix} ]

を得る。よって、

A = LU = [ \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & 1 \end{matrix} ] [ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{matrix} ]

また、これを作成した LU 分解をする関数で解かせてみた結果は以下の通りである。当然だが、同じ結果が得られている。

A=
1.000000  2.000000  3.000000
2.000000  8.000000  11.000000
3.000000  22.000000  35.000000
b=
1.000000
0.000000
0.000000
L=
1.000000  0.000000  0.000000
2.000000  1.000000  0.000000
3.000000  4.000000  1.000000

U=
1.000000  2.000000  3.000000
0.000000  4.000000  5.000000
0.000000  0.000000  6.000000
Press any key to continue

2.1.1 問題３−５ (1)

まず、 Ux = y とおいて、 Ly = b を解いて y = L^{−1}b を求める。

Ly = b

より、

[ \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & 1 \end{matrix} ] [ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ] = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ]

これより、まず、

y_{1} = 1

を得る。これを使って、

2y_{2} + y_{2} = 0
y_{2} = −2y_{1} = −2 × 1 = −2

また、

3y_{1} + 4y_{2} + y3 = 0
y_{3} = −3y_{1} − 4y_{2}
= −3 × 1 − 4 × (−2)
= 5

よって

[ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ] = [ \begin{matrix}1 \\ −2 \\ 5 \end{matrix} ]

を得る。次に、 Ux = y から、 x = U^{−1}y を求める。

Ux = y

より、

[ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{matrix} ] [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ \begin{matrix}1 \\ −2 \\ 5 \end{matrix} ]

これも同様にして、

x_{3} = \frac{5}{6}

4x_{2} + 5x_{3} = −2
x_{2} = −\frac{37}{24}
x_{1} + 2x_{2} + 3x_{3} = 1

x_{1} = \frac{19}{12}

よって、

[ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ −\frac{37}{24}_{\frac{19}{12}}^{\frac{5}{6}} ]

また、これを、今回作成した LU 分解で連立一次方程式の近似解を求めるプログラムで解かせてみた結果は以下のとおりである。きちんと近似した解が得られている。

A=
1.000000  2.000000  3.000000
2.000000  8.000000  11.000000
3.000000  22.000000  35.000000
b=
1.000000
0.000000
0.000000
L=
1.000000  0.000000  0.000000
2.000000  1.000000  0.000000
3.000000  4.000000  1.000000
U=
1.000000  2.000000  3.000000
0.000000  4.000000  5.000000
0.000000  0.000000  6.000000

答え=

1.583333
-1.541667
0.833333
Press any key to continue

2.1.2 問題３−５ (2)

問題３−５ (1) と同様にして、 Ly = b より

[ \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & 1 \end{matrix} ] [ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ] = [ \begin{matrix}1 \\ 2 \\ 3 \end{matrix} ]

これより

[ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ] = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ]

次に、 Ux = y より

[ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{matrix} ] [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ]

これより、

[ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ]

を得る。

また、これを、今回作成した LU 分解で連立一次方程式の近似解を求めるプログラムで解かせてみた結果は以下の通りである。同じ解が得られている。

A=
1.000000  2.000000  3.000000
2.000000  8.000000  11.000000
3.000000  22.000000  35.000000
b=
1.000000
2.000000
3.000000
L=
1.000000  0.000000  0.000000
2.000000  1.000000  0.000000
3.000000  4.000000  1.000000
U=
1.000000  2.000000  3.000000
0.000000  4.000000  5.000000

0.000000  0.000000  6.000000

答え=

1.000000
0.000000
0.000000
Press any key to continue

2.1.3 問題３−５ (3)

問題３−５ (1),(2) と同様にして、 Ly = b より

[ \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & 1 \end{matrix} ] [ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ] = [ \begin{matrix}−1 \\ 0 \\ 1 \end{matrix} ]

これより

[ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ] = [ \begin{matrix}−1 \\ 2 \\ −4 \end{matrix} ]

次に、 Ux = y より

[ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{matrix} ] [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ \begin{matrix}−1 \\ 2 \\ −4 \end{matrix} ]

これより、

[ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] = [ −\frac{5}{3} ; \frac{4}{3} ; −\frac{2}{3} ]

を得る。

また、これを、今回作成した LU 分解で連立一次方程式の近似解を求めるプログラムで解かせてみた結果は以下のとおりである。きちんと近似した解が得られている。

A=
1.000000  2.000000  3.000000
2.000000  8.000000  11.000000
3.000000  22.000000  35.000000
b=
-1.000000
0.000000
1.000000
L=
1.000000  0.000000  0.000000
2.000000  1.000000  0.000000
3.000000  4.000000  1.000000
U=
1.000000  2.000000  3.000000
0.000000  4.000000  5.000000

答え=

-1.666667
1.333333
-0.666667
Press any key to continue

3 性質３の証明

3.1 補題

M を n 次の正方行列、 I を単位行列とするとき、 ||M|| < 1 ならば、 (I − M)^{−1} が存在して次のことが成り立つ。

||(I − M)^{−1}|| ≤ \frac{1}{1 − ||M||}

3.2 補題の証明

（証明）

||M|| < 1 より

lim_{n→∞} ||M^{n}|| ≤ lim_{n→∞} ||M||^{n} = 0 これより、 lim_{n→∞} M^{n} = O

また、

(I − M) ∑_{k=0}^{n} M^{k} = (I − M) · \frac{I − M^{n+1}}{I − M} = I − M^{n+1}

以上より

(I − M) ∑_{n=0}^{∞} M^{n} = lim_{n→∞}(I − M) ∑_{k=0}^{n} M^{k} = lim_{n→∞}(I − M^{n+1}) = I

ゆえに、 (I − M)^{−1} の存在が示され、また、 (I − M)^{−1} は

(I − M)^{−1} = ∑_{n=0}^{∞} M^{n}

(I − M)^{−1} = ∑_{n=0}^{∞} M^{n} = ∑_{n=0}^{m} M^{n} + M^{m+1} ∑_{n=0}^{∞} M^{n}

(I − M)^{−1} = ∑_{n=0}^{m} M^{n} + M^{m+1}(I − M)^{−1}
||(I − M)^{−1}|| = || ∑_{n=0}^{m} M^{n} + M^{m+1}(I − M)^{−1}||
≤ || ∑_{n=0}^{m} M^{n}|| + ||M^{m+1}(I − M)^{−1}||
≤ ∑_{n=0}^{m} ||M^{n}|| + ||M^{m+1}||||(I − M)^{−1}||
≤ ∑_{n=0}^{m} ||M||^{n} + ||M||^{m+1}||(I − M)^{−1}||

ゆえに

(1 − ||M||^{m+1})||(I − M)^{−1}|| ≤ ∑_{n=0}^{m} ||M||^{n} = \frac{1 − ||M||^{m+1}}{1 − ||M||}

ここで、 ||M|| < 1 より、 1 − ||M||^{m+1} > 0 だから

||(I − M)^{−1}|| ≤ \frac{1}{1 − ||M||}

（証明終わり）

3.3 問題の証明

（証明） ||RA − I|| = ||I − RA|| < 1 より M = I − RA とすると補題により

(I − M)^{−1} = (RA)^{−1} = A^{−1}R^{−1}

が存在する。

また、 A^{−1}R^{−1} が存在するということは、同時に A^{−1} が存在するということを示していることに他ならないよって、 A は正則である。

次に、

||A^{−1}R^{−1}|| ≤ \frac{1}{1 − ||RA − I||}

両辺に ||R|| を掛けて

||A^{−1}R^{−1}||||R|| ≤ \frac{||R||}{1 − ||RA − I||}

ここで

||A^{−1}|| = ||A^{−1}I|| = ||A^{−1}R^{−1}R|| ≤ ||A^{−1}R^{−1}||||R|| だから

||A^{−1}|| ≤ \frac{||R||}{1 − ||RA − I||}

（証明終わり）

4 性質５の証明

4.1 解答・証明

（証明）

行列 A の固有ベクトルを u とすると Au = λu と表せる。ここで固有ベクトル u の要素の中で絶対値が最大のものを u_{i} とすると、

(λ − a_{ii})u_{i} = ∑_{k=1(k≠i)}^{n} a_{ik}u_{k}

したがって、

|λ − a_{ii}| ≤ ∑_{k=1(k≠i)}^{n} |a_{ik}||\frac{u_{k}}{u_{i}} |

|u_{i}| が最大であるから、 |\frac{u_{k}}{u_{i}} | ≤ 1 である。したがって

|λ − a_{ii}| ≤ ∑_{k=1(k≠i)}^{n} |a_{kj}|

（証明終り）

5 性質６の証明

5.1 問題の証明

まず、問題の for 文を以下のようなものと考える。

for(i = 1; i <= n; i + +)
t_{i} = x_{i} ∗ y_{i};
for(i = 1; i <= n; i + +)
t+ = t_{i};

また、各成分の実数での計算結果を a_{i} = x_{i} ∗ y_{i} とすると、 t_{i} = a_{i}(1 + δ_{i}) である。

したがって、

t_{1} = a_{1}(1 + δ_{1})
t_{2} = (1 + δ_{3}) {a_{1}(1 + δ_{1}) + a_{2}(1 + δ_{2})}
= a_{1}(1 + δ_{1})(1 + δ_{3}) + a_{2}(1 + δ_{2})(1 + δ_{3})
t_{3} = (1 + δ_{5}) {a_{1}(1 + δ_{1})(1 + δ_{3}) + a_{2}(1 + δ_{2})(1 + δ_{3}) + a_{3}(1 + δ_{4})}
= a_{1}(1 + δ_{1})(1 + δ_{3})(1 + δ_{5}) + a_{2}(1 + δ_{2})(1 + δ_{3})(1 + δ_{5}) + a_{3}(1 + δ_{4})(1 + δ_{5})

これを t_{n} まで繰り返すと、

t_{n} = a_{1} \underbrace{(1 + δ_{d}) ⋯ (1 + δ_{j})}_{n 個}+a_{2} \underbrace{(1 + δ_{p}) ⋯ (1 + δ_{j})}_{n 個}+a_{3} \underbrace{(1 + δ_{h}) ⋯ (1 + δ_{j})}_{n−1 個}+ ⋯ + a_{n} \underbrace{(1 + δ_{f})(1 + δ_{j})}_{2 個}

（ただし、 d, p, h, f, j は任意に用いた）

ここで、 1 + δ_{i} の積を上のように書き連ねるのは不便であるから、 1 + δ_{i} の k 個の積を便宜的に次のように書くことにする（厳密に言うと δ の添え字の番号は異なってしまうのだが、内容的には同じことなので、ここでは気にしないことにする）。

∏_{i=1}^{k}(1 + δ_{i}) = 1 + b_{k}

ここで、 b_{k} について以下のことが成り立つ。

|b_{k}| ≤ \frac{ku}{1 − ku}

このことを証明する。（証明） (i)k = 1 のとき

1 + δ_{1} = 1 + b_{1} ⟺ b_{1} = δ_{1}

|b_{1}| = |δ_{1}| ≤ u ≤ \frac{u}{1 − u}

\frac{u}{1 − u} − u = \frac{u^{2}}{1 − u} ≥ 0

⟺ u ≤ \frac{u}{1 − u}

よって、

|b_{1}| ≤ \frac{u}{1 − u}

(ii)k = m のとき

|b_{m}| ≤ \frac{mu}{1 − mu}

が成り立つと仮定する。

1 + b_{m+1} = ∏_{i=1}^{m+1}(1 + δ_{i}) = (1 + δ_{m+1}) ∏_{i=1}^{m}(1 + δ_{i})

⟺ 1 + b_{m+1} = (1 + δ_{m+1})(1 + b_{m})
⟺ b_{m+1} = δ_{m+1} + b_{m}(1 + δ_{m+1})

したがって、

|b_{m+1}| = |δ_{m+1} + b_{m}(1 + δ_{m+1})|
≤ |δ_{m+1}| + |b_{m}||1 + δ_{m+1}|
≤ u + \frac{mu}{1 − mu}(1 + u)
= \frac{(m + 1)u}{1 − mu}

よって、

0 < 1 − (m + 1)u < 1 − mu < 1

だから、

\frac{(m + 1)u}{1 − mu} ≤ \frac{(m + 1)u}{1 − (m + 1)u}

ゆえに

|b_{m+1}| ≤ \frac{(m + 1)u}{1 − (m + 1)u}

これは、 k = m + 1 のときも成り立つことを示している。 (i), (ii) より、任意の k に対して成立する。

（証明終わり）

これを用いると、 t = t_{n} = a_{1}(1 + b_{n}^{(1)}) + a_{2}(1 + b_{n}^{(2)}) + a_{3}(1 + b_{n−1}) + ⋯ + a_{n−1}(1 + b_{3}) + a_{n}(1 + b_{2})

s − t = −(a_{1}b_{n}^{(1)} + a_{2}b_{n}^{(2)} + a_{3}b_{n−1} + ⋯ + a_{n−1}b_{3} + a_{n}b_{2})
|s − t| = |a_{1}b_{n}^{(1)} + a_{2}b_{n}^{(2)} + a_{3}b_{n−1} + ⋯ + a_{n−1}b_{3} + a_{n}b_{2}|
≤ |a_{1}b_{n}^{(1)}| + |a_{2}b_{n}^{(2)}| + |a_{3}b_{n−1}| + ⋯ + |a_{n−1}b_{3}| + |a_{n}b_{2}|
= |a_{1}||b_{n}^{(1)}| + |a_{2}||b_{n}^{(2)}| + |a_{3}||b_{n−1}| + ⋯ + |a_{n−1}||b_{3}| + |a_{n}||b_{2}|

ここで、 |b_{n}| ≤ |b_{n+1}| だから、 b_{n}^{(1)} = b_{n}^{(2)} = b_{n} とすると、

|s − t| ≤ |a_{1}||b_{n}| + |a_{2}||b_{n}| + |a_{3}||b_{n}| + ⋯ + |a_{n−1}||b_{n}| + |a_{n}||b_{n}|
= |b_{n}|(|a_{1}| + |a_{2}| + |a_{3}| + ⋯ + |a_{n−1}| + |a_{n}|)
= |b_{n}| ∑_{i=1}^{n} |a_{i}|
≤ γ_{n} ∑_{i=1}^{n} |a_{i}|

したがって、

|s − t| ≤ γ_{n} ∑_{i=1}^{n} |x_{i}y_{i}| γ_{n} = \frac{nu}{1 − nu} (nu < 1)

である。（証明終わり）

6 作成したプログラム

6.1 lu.h

/*

行列ライブラリ

（ＬＵ分解法、逆行列、ガウスの消去法付き）

  header file
  2001-6-15
  Programmed by Takuya Saito
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct {
    int gyou_width;
    int retu_width;
    double** matrix;
} MATRIX;
/*

関数のプロトタイプ宣言

*/
int alloc_matrix(MATRIX* m, int gyou, int retu);
void free_matrix(MATRIX* m);
void lu_solve(MATRIX* u, MATRIX* l, double* b);
void lu_decompose(MATRIX* u, MATRIX* l);
void disp_matrix(MATRIX* a);
void disp_vector(double* b, int size);
void Gauss(MATRIX* a, double* b, int pivot);
int add_matrix(MATRIX* a, MATRIX* b, MATRIX* c);
int sub_matrix(MATRIX* a, MATRIX* b, MATRIX* c);
int mul_matrix(MATRIX* a, MATRIX* b, MATRIX* c);
void clear_matrix(MATRIX* m);
int t_matrix(MATRIX* a, MATRIX* t);
int naiseki_matrix(MATRIX* a, MATRIX* b, double* c);

void set_b(double* b, int size, int n);
void set_inv(MATRIX* inv, double* b, int n);

6.2 lu.c

/*

行列ライブラリ

（ＬＵ分解法、逆行列、ガウスの消去法付き）

  2001-6-15
  Programmed by Takuya Saito
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "lu.h"

/*

MATRIX* m に、行、列で指定されただけメモリを確保する。

   また、 m->x_width, m->y_width に値を設定する。

（メモリ確保に失敗したら０を、成功したら１を返す）

*/
int alloc_matrix(MATRIX* m, int gyou, int retu)
{
    int i;
    m->gyou_width = gyou;
    m->retu_width = retu;
    if ((m->matrix = (double**) malloc(sizeof(double*) * m->retu_width))
        == NULL)
        return 0;
    if ((m->matrix[0] = (double*) malloc(sizeof(double) * m->retu_width * m->gyou_width))
        == NULL)
        return 0;
    for (i=1; i < m->retu_width; i++) {
        m->matrix[i] = m->matrix[i-1] + m->gyou_width;

}
return 1;
}
/*

MATRIX* の領域を開放する関数

*/
void free_matrix(MATRIX* m)
{
free(m->matrix[0]);
free(m->matrix);
}
/*

ＬＵ分解をする関数

*/
void lu_decompose(MATRIX* u, MATRIX* l)
{
int i, j, k;
double temp;
/*

行列 L を作成

*/
for (i=0; i < l->gyou_width; i++) {
    for (j=0; j < l->retu_width; j++) {
        if (i == j) l->matrix[i][j] = 1;
        else l->matrix[i][j] = 0;
    }
}
/*

前進消去

*/
for (i=0; i < u->retu_width; i++) {
    for (j=i+1; j < u->gyou_width; j++) {
        /* j 行目から、 i 行目の a[j][i] 倍を引く。 */

        temp = u->matrix[j][i] / u->matrix[i][i];
        l->matrix[j][i] = temp;
        for (k=i; k < u->retu_width; k++) {
            u->matrix[j][k] -= temp * u->matrix[i][k];
        }
    }
}
}
/*

ＬＵ分解の結果を使って方程式を解く関数

*/
void lu_solve(MATRIX* u, MATRIX* l, double* b)
{
int i, j;
/*

前進代入

*/
for (i=1; i < l->gyou_width; i++) {
    for (j=0; j < i; j++) {
        b[i] -= l->matrix[i][j]*b[j];
    }
}

/*

後退代入

*/
for (i=u->gyou_width-1; i >=0; i--) {
    for (j=i+1; j < u->retu_width; j++) {
        b[i] -= u->matrix[i][j] * b[j];
    }
    b[i] = b[i] / u->matrix[i][i];
}
}
/*

b に、 n で指定した所だけ１をいれ、後は全部０を入れる関数

*/
void set_b(double* b, int size, int n)
{
    int i;
    for (i=0; i < size; i++) {
        if (i == n)
            b[i] = 1;
        else
            b[i] = 0;
    }
}
/*

ベクトル b の値を inv 行列の n 列目にコピーする関数

*/
void set_inv(MATRIX* inv, double* b, int n)
{
    int i;
    for (i=0; i < inv->retu_width; i++) {
        inv->matrix[i][n] = b[i];
    }
}
/*

ガウスの消去法で連立方程式を解く

*/
void Gauss(MATRIX* a, double* b, int pivot)
{
    int i, j, k, q, q2, w;
    double temp;
    /*

前進消去

    */

for (i=0; i < a->gyou_width; i++) {
/*

a[i][i] 〜 a[n-1][i] で絶対値が最大のものを探し、その行と i 行目とを

入れ替える pivot != 0 ならこれを実行する

*/
if (pivot != 0) {
/*

部分ピボット選択ありのとき

*/
q2 = i;
for (q=i+1; q < a->gyou_width; q++) {
    if (fabs(a->matrix[i][i]) < fabs(a->matrix[q][i]))
        q2 = q;
}
if (a->matrix[q2][i] == 0.0) {
    printf("fault!\n");
    return;
}
if (q2 != i) {
    for (w=i; w < a->retu_width; w++) {
        temp = a->matrix[i][w];
        a->matrix[i][w] = a->matrix[q2][w];
        a->matrix[q2][w] = temp;
    }
    temp = b[i];
    b[i] = b[q2];
    b[q2] = temp;
}
} else {
/*

ピボット選択なしのとき

*/

if (a->matrix[i][i] == 0.0) {
q2 = i;
for (q=i+1; q < a->gyou_width; q++) {
if (a->matrix[q][i] != 0.0) {
    q2 = q;
    break;
}
}
if (q2 == i)
printf("Can't Solve !!!\n");
else {
for (w=i; w < a->retu_width; w++) {
    temp = a->matrix[i][w];
    a->matrix[i][w] = a->matrix[q2][w];
    a->matrix[q2][w] = temp;
}
temp = b[i];
b[i] = b[q2];
b[q2] = temp;
}
}
}
/*
i 行目を a[i][i] で割る
*/
temp = a->matrix[i][i];
for (j=i; j < a->retu_width; j++) {
a->matrix[i][j] /= temp;
}
b[i] /= temp;
for (j=i+1; j < a->gyou_width; j++) {
/*

j 行目から、 i 行目の a[j][i] 倍を引く。

*/
temp = a->matrix[j][i];

        for (k=i; k < a->retu_width; k++) {
            a->matrix[j][k] -= temp * a->matrix[i][k];
        }
        b[j] -=  temp * b[i];
    }
}
/*

後退代入

*/
for (i=a->gyou_width - 2; i >=0; i--) {
    for (j=i+1; j < a->retu_width; j++) {
        b[i] -= a->matrix[i][j] * b[j];
    }
}
}
/*

行列の足し算

c = a + b

戻り値：成功すれば１、失敗すれば０を返す

*/
int add_matrix(MATRIX* a, MATRIX* b, MATRIX* c)
{
int x, y;
if (a->retu_width != b->retu_width || b->retu_width != c->retu_width)
    return 0;
if (a->gyou_width != b->gyou_width || b->gyou_width != c->gyou_width)
    return 0;
for (x = 0; x < a->retu_width; x++) {
    for (y=0; y < a->gyou_width; y++) {
        c->matrix[y][x] = a->matrix[y][x] + b->matrix[y][x];
    }
}

    return 1;
}
/*

行列の引き算

  c = a - b

戻り値：成功すれば１、失敗すれば０を返す

*/
int sub_matrix(MATRIX* a, MATRIX* b, MATRIX* c)
{
    int x, y;
    if (a->retu_width != b->retu_width || b->retu_width != c->retu_width)
        return 0;
    if (a->gyou_width != b->gyou_width || b->gyou_width != c->gyou_width)
        return 0;
    for (x = 0; x < a->retu_width; x++) {
        for (y=0; y < a->gyou_width; y++) {
            c->matrix[y][x] = a->matrix[y][x] - b->matrix[y][x];
        }
    }
    return 1;
}
/*

行列の掛け算

  c = a * b

戻り値：成功すれば１、失敗すれば０を返す

*/
int mul_matrix(MATRIX* a, MATRIX* b, MATRIX* c)
{
    int x, y, z;

/* a の列と b の行が同じか？ */

if (a->retu_width != b->gyou_width)
    return 0;

/* a の行と c の行、 b の列と c の列が同じか？ */

if (a->gyou_width != c->gyou_width || b->retu_width != c->retu_width)
    return 0;
clear_matrix(c);
for (x=0; x < c->retu_width; x++) {
    for (y=0; y < c->gyou_width; y++) {
        for (z=0; z < a->retu_width; z++) {
            c->matrix[y][x]
                += a->matrix[y][z] * b->matrix[z][x];
        }
    }
}
return 1;
}
/*

行列の要素を全部０にする関数

*/
void clear_matrix(MATRIX* m)
{
int x,y;
for (y=0; y < m->gyou_width; y++) {
    for (x=0; x < m->retu_width; x++) {
        m->matrix[y][x] = 0.;
    }
}
}
/*

転置行列を作る関数

t に a の転置行列を入れる

戻り値：問題なければ１、問題があり転置行列が作れなければ０

*/
int t_matrix(MATRIX* a, MATRIX* t)
{
int x, y;
if (a->retu_width != t->gyou_width)
    return 0;
if (a->gyou_width != t->retu_width)
    return 0;
for (x=0; x < t->retu_width; x++) {
    for (y=0; y < t->gyou_width; y++) {
        t->matrix[y][x] = a->matrix[x][y];
    }
}
return 1;
}
/*

行列の内積を求める関数

c = <a|b>

戻り値：問題が無ければ１。引数の *c に内積を入れて戻る問題があり内積が求められないときは０

*/
int naiseki_matrix(MATRIX* a, MATRIX* b, double* c)
{
MATRIX* t;
MATRIX* q;
int i;

/* a と b の行数と列数は同じでなければ内積は求まらない */

if (a->retu_width != b->retu_width || a->gyou_width != b->gyou_width)
    return 0;

/* a と b の列数が同じでなければ内積は求まらない */

if (a->gyou_width != b->gyou_width)
return 0;

/* a の転置行列 t を作る（メモリ領域も確保する) */

if ( (t = (MATRIX*) malloc(sizeof(MATRIX))) == NULL)
return 0;
if (alloc_matrix(t, a->gyou_width, a->retu_width) == 0) {
free(t);
return 0;
}
if (t_matrix(a, t) == 0) {
free_matrix(t);
free(t);
return 0;
}

/* t*b の結果を入れる行列 q を作る（メモリ領域も確保する） */

if ( (q = (MATRIX*) malloc(sizeof(MATRIX))) == NULL) {
free_matrix(t);
free(t);
return 0;
}
if (alloc_matrix(q, t->gyou_width, b->retu_width) == 0) {
free_matrix(t);
free(t);
free(q);
return 0;
}
/* q = t*b */
if (mul_matrix(t, b, q) == 0) {
free_matrix(q);
free_matrix(t);
free(q);
free(t);
return 0;
}

/* c = tr(q) */
*c = 0;
for (i=0; i < q->retu_width; i++) {
    *c += q->matrix[i][i];
}

/* メモリ開放 */

free_matrix(q);
free_matrix(t);
free(q);
free(t);
return 1;
}
/*

行列の成分を表示する関数

*/
void disp_matrix(MATRIX* a)
{
int i,j;
for (i=0; i < a->gyou_width; i++) {
    for (j=0; j < a->retu_width; j++) {
        printf("%f  ", a->matrix[i][j]);
    }
    printf("\n");
}
}
/*

ベクトルの成分を表示する関数

*/
void disp_vector(double* b, int size)
{
int i;
for (i=0; i < size; i++) {

printf("%f\n", b[i]);
}
}

6.3 gauss.c

/*

ガウスの消去法で方程式を解くプログラム

2001-6-15
Programmed by Takuya Saito
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "lu.h"
int main()
{
int size;
MATRIX* a;

/* 方程式の係数 (Ax=b の b) をセットする */

double b[3] = {
    {1.0},
    {0.0},
    {0.0},
};
a = (MATRIX*) malloc( sizeof(MATRIX) );
size = 3;  /* 行列のサイズ指定 */

/* 行列のメモリ領域確保 */

alloc_matrix(a, size, size);

/* 解く方程式の係数（Ax=b の A）を代入する */

a->matrix[0][0] = 1.0;
a->matrix[0][1] = 2.0;
a->matrix[0][2] = 3.0;
a->matrix[1][0] = 2.0;

a->matrix[1][1] = 8.0;
a->matrix[1][2] = 11.0;
a->matrix[2][0] = 3.0;
a->matrix[2][1] = 22.0;
a->matrix[2][2] = 35.0;

/* 解く方程式の係数（行列）を表示 (Ax=b の A) */

printf("A=\n");
disp_matrix(a);
printf("\nb=\n");
disp_vector(b, size);

/* ガウスの消去法で Ax=b を解き、答えを b に代入する */

Gauss(a, b, 1);

/* 解を表示する */

printf("\n 答え\n");
disp_vector(b, size);
free_matrix(a);
free(a);
return 0;
}

6.4 lu.c

/*-------------------------------------------------

ＬＵ分解で連立方程式を解くプログラム

2001-6-15
Programmed by Takuya Saito
--------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "lu.h" /* 行列演算ライブラリ・ヘッダ */

int main()
{
int size;
double* b;
MATRIX* u;
MATRIX* l;
size = 3; /* 行列のサイズ指定 */
u = (MATRIX*) malloc( sizeof(MATRIX) );
l = (MATRIX*) malloc( sizeof(MATRIX) );
b = (double*) malloc( sizeof(double) * size);
/*

行列用メモリ領域確保

*/
alloc_matrix(u, size, size);
alloc_matrix(l, size, size);
/*

解く方程式の係数を配列に代入する

*/
u->matrix[0][0] = 2.0;
u->matrix[0][1] = 6.0;

u->matrix[0][2] = 4.0;
u->matrix[1][0] = 3.0;
u->matrix[1][1] = 7.0;
u->matrix[1][2] = 2.0;
u->matrix[2][0] = 1.0;
u->matrix[2][1] = 2.0;
u->matrix[2][2] = 3.0;
b[0] = 8.0;
b[1] = 6.0;
b[2] = 7.0;
/* 定義された行列表示 (Ax=b の A)*/
printf("A=\n");
disp_matrix(u);

/* 定義された n 次ベクトル表示 (Ax=b の b) */

printf("\nb=\n");
disp_vector(b, size);
/* LU 分解する */
lu_decompose(u, l);
/* L 行列表示 */
printf("\nL=\n");
disp_matrix(l);
/* U 行列表示 */
printf("\nU=\n");
disp_matrix(u);
printf("\n 答え=\n");

/* LU 分解された結果を用い、解を求める */

lu_solve(u, l, b);
/* 解を表示 */
disp_vector(b, size);

/* 行列＆ベクトル用メモリ領域開放 */

free_matrix(u);
free_matrix(l);
free(l);
free(u);
free(b);
return 0;
}

6.5 inverse.c

/*--------------------------------------------------

ＬＵ分解を用いて逆行列を求める

2001-6-15
Programmed by Takuya Saito
--------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include "lu.h"
int main()
{
int i, size;
double* b;
MATRIX* u;
MATRIX* l;
MATRIX* inv;
size = 3; /* 行列のサイズ指定 */
u = (MATRIX*) malloc( sizeof(MATRIX) );
l = (MATRIX*) malloc( sizeof(MATRIX) );
inv = (MATRIX*) malloc( sizeof(MATRIX) );
b = (double*) malloc( sizeof(double) * size);

/* 行列用メモリ領域確保 */

alloc_matrix(u,   size, size);
alloc_matrix(l,   size, size);
alloc_matrix(inv, size, size);

/* 逆行列を求める行列の係数を指定する */

u->matrix[0][0] = 1.0;
u->matrix[0][1] = 2.0;
u->matrix[0][2] = -4.0;
u->matrix[1][0] = 2.0;
u->matrix[1][1] = 1.0;

u->matrix[1][2] = -4.0;
u->matrix[2][0] = 1.0;
u->matrix[2][1] = -2.0;
u->matrix[2][2] = 1.0;

/* 逆行列を求める行列を表示する */

printf("A=\n");
disp_matrix(u);

/* ＬＵ分解する */

lu_decompose(u, l);

/* Ｕ行列表示 */

printf("\nU=\n");
disp_matrix(u);

/* Ｌ行列表示 */

printf("\nL=\n");
disp_matrix(l);
/*

上で作ったＬ、Ｕ行列を使い、逆行列を求める

*/
for (i=0; i < size; i++) {

/* i 次だけ 1 で、後は０のベクトル b を作る */

    set_b(b, size, i);

/* L,U 行列を使い、 LUx=b の解を求め、 b に答えを代入する */

    lu_solve(u, l, b);

/* inv 行列の i 列目にベクトル b の値をコピーする */

    set_inv(inv, b, i);
}
/* 逆行列（inv）表示 */
printf("\n 逆行列=\n");
disp_matrix(inv);
free_matrix(u);

free_matrix(l);
free_matrix(inv);
free(inv);
free(b);
free(l);
free(u);
return 0;
}