目次

1 問題 3

1.1 問題 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.1 問題 3-1 . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.2 問題 3-2 . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.3 問題 3-3 . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.4 問題 3-4 . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.5 問題 3-5 . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2 問題 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 問題 1 の問題 3-1 の解答および考察 6 2.1 Gauss の消去法の概要 . . . . . . . . . . . . . . . . . . . . . . 6 2.2 プログラムへの実現の概要 . . . . . . . . . . . . . . . . . . . . 8

2.2.1 前進消去法段階について . . . . . . . . . . . . . . . . . 9 2.2.2 ピボット選択について . . . . . . . . . . . . . . . . . . 9 2.2.3 後置代入段階について . . . . . . . . . . . . . . . . . . 10

2.3 プログラム 1 （ピボット選択なし） . . . . . . . . . . . . . . . 11 2.4 プログラム 2 （部分ピボット選択あり） . . . . . . . . . . . . . 13 2.5 プログラムの解説 . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.5.1 プログラムの最初の部分について . . . . . . . . . . . . 16 2.5.2 main 関数について . . . . . . . . . . . . . . . . . . . . 16 2.5.3 連立方程式の入力と表示について . . . . . . . . . . . . 17 2.5.4 前進消去段階について . . . . . . . . . . . . . . . . . . 18 2.5.5 絶対値の計算について（プログラム 2 のみ） . . . . . . 19 2.5.6 後退代入段階について . . . . . . . . . . . . . . . . . . 19

2.6 プログラム 1 の実行結果（ピボット選択なし） . . . . . . . . . 20 2.7 考察 (ピボット選択なし) . . . . . . . . . . . . . . . . . . . . . 21 2.8 プログラム 2 の実行結果（部分ピボット選択あり） . . . . . . 21 2.9 考察 (部分ピボット選択あり) . . . . . . . . . . . . . . . . . . 22 2.10 補足 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.10.1 前進消去の工夫 . . . . . . . . . . . . . . . . . . . . . . 22 2.10.2 絶対値の計算（プログラム 2 のみ） . . . . . . . . . . . 23

3 問題 1 の問題 3-2 の解答および考察 23

3.1 LU 分解の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1.1 三角行列について . . . . . . . . . . . . . . . . . . . . . 23 3.1.2 LU 分解について . . . . . . . . . . . . . . . . . . . . . 24 3.1.3 クラウト法について . . . . . . . . . . . . . . . . . . . 24 3.1.4 上記のアルゴリズムの解説 . . . . . . . . . . . . . . . . 25 3.2 LU 分解の連立一次方程式への応用 . . . . . . . . . . . . . . . 26 3.3 プログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4 プログラムの解説 . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.4.1 プログラムの最初の部分について . . . . . . . . . . . . 30 3.4.2 LU 分解を行う関数 decomp について . . . . . . . . . . 31 3.4.3 連立一次方程式の解の算出を行う関数 solve について . 32 3.4.4 main 関数について . . . . . . . . . . . . . . . . . . . . 33

3.5 実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.5.1 問題 3-5-1 に対応する実行結果 . . . . . . . . . . . . . 34 3.5.2 問題 3-5-2 に対応する実行結果 . . . . . . . . . . . . . 36 3.5.3 問題 3-5-3 に対応する実行結果 . . . . . . . . . . . . . 37

4 問題 1 の問題 3-3 の解答および考察 38 4.1 LU 分解の逆行列への応用 . . . . . . . . . . . . . . . . . . . . 38 4.2 プログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.3 プログラムの解説 . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.4 実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5 問題 1 の問題 3-4 の解答および考察 44

6 問題 1 の問題 3-5 の解答および考察 45 6.1 問題 1 の問題 3-5-1 の解答および考察 . . . . . . . . . . . . . . 45 6.2 問題 1 の問題 3-5-2 の解答および考察 . . . . . . . . . . . . . . 47 6.3 問題 1 の問題 3-5-3 の解答および考察 . . . . . . . . . . . . . . 49

7 問題 2 の解答および考察 51

7.1 性質 3 の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 7.1.1 性質 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 7.1.2 証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

7.2 性質 5 の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 7.2.1 性質 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 7.2.2 証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

7.3 性質 6 の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 7.3.1 性質 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 7.3.2 証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

8 感想 55

数値計算:第 3 回レポート

{G99P143-0 \atop 森野岳宏}

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

1 問題

1.1 問題 1

教科書「数値計算」 (p.31) の問題 3(1〜5) を解け。問題 3 の内容は以下の通り。

1.1.1 問題 3-1

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

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

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

1.1.4 問題 3-4

行列

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

の LU 分解を求めよ。

1.1.5 問題 3-5

問題 3-4 の結果を用いて、

1.

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

2.

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

3.

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

として、

Ax = b

を解け。

1.2 問題 2

性質 3,5,6 を証明せよ。（ただし、性質 1,2,4 は仮定してよい）性質 1,2,3,4,5,6 は以下の通り。

• 性質 1 n × n 行列 A が ||A|| < 1 を満たすとする。このとき、 I を単位行列とすると、 I + A は正則行列となる。

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

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

が成り立つ。

(||A|| = sup_{||x||=1} ||Ax||)

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

||RA − 1|| < 1

ならば A は正則であり、

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

が成り立つ。

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

|λ| ≤ ||A||

が成り立つ。

• 性質 5 A:n × n 複素行列とする。このとき、行列 A の個有値は、

K = {λ ∈ C||λ − a_{ij}| ≤ ∑_{k≠i^{k=1}}^{n} |a_{ik}|}

の中に存在する。

• 性質 6 IEEE754 の下で浮動小数点演算の下で議論する。最近点での丸め (nearest) のモードで考える。

u = 2^{−54}

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

x, y ∈ F, 0 ∈ +, −, ∗, /
x ∘ y = (x ∘ y)(1 + δ)
|δ| ≤ u

となる。これを利用すると、 x, y ∈ F^{n} とすると、 s = xy(実数としての計算) に対して、

near();
t=0;
for(i=1;i<=n;i++)
    t+=xi*yi;

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

|x − t| ≤ δ_{n} ∑_{i=1}^{n} |x_{i}y_{i}|

δ_{n} = \frac{nu}{1 − nu}

2 問題 1 の問題 3-1 の解答および考察

2.1 Gauss の消去法の概要

Gauss の消去法の流れを簡単な例をもとに解説する。まず、以下のような簡単な式（3 つの変数 x,y,z と 3 つの方程式を考える。 )

x + 3 − 3y − 4z = 8
x + y − 2z = 2
−x − 2y + 5z = −1

我々の目標は、これらの方程式を同時に満たす変数の値を求めることである。場合によっては、解がないことや（例えば、 x + y = 1,x + y = 2 のように 2 つの方程式が矛盾する時）、解がいくつもあることもある（例えば 2 つの方程式が同じとか、方程式の数より変数の個数が多い時）。方程式の個数と変数の個数は等しいものとし、会があるときは一意な解を求める。 3 つ以上の変数の場合も扱いやすくするために添え字を使って変数名をつける。

x_{1} + 3x_{2} − 4x_{3} = 8
x_{1} + x_{2} − 2x_{3} = 2
−x_{1} − 2x_{2} + 5x_{3} = −1

変数を何回も繰り返して書くのを避けるために、連立方程式を行列で表すと都合がよい。上の方程式は次の行列の方程式と同じである。

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

これらの方程式に対して解を変えないで実行できる操作がいくつかある。 • 方程式の入れ替え：明らかに、方程式を書き下す順番は解に影響しな

い。行列の表現では、この操作は行列（と右辺ベクトル）の行の入れ替えに対応する。

• 変数名の付け替え：この操作は、行列の上では列の入れ替えに対応する。（列 i と j を入れ替えた時は変数 x_{i} と x_{j} も入れ替えなけらばならない。）

• 式の定数倍：行列（と右辺と右辺ベクトル）上のある行を定数倍する。

• 2 つの式の和： 2 つの式の和をとり、一方の式をその和と入れ替える。

これらの操作、とくに最後の操作が解に影響を与えないということを納得するには少し考える必要がある。例えば、上の方程式系で最初の 2 つの方程式の差で 2 番目の式を置き換えると同値な方程式系

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

が得られる。これによって 2 番目の式から x_{1} が消去されることに注意されたい。同様に、第 1 式を第 3 式に加えることによって第 3 式から x_{1} を消去できる。

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

ここで、変数 x_{1} は第１式以外のすべての式から消去されている。このような系統的のやり方を進めると、もとの方程式系と同じ解をもち、もっと解きやすいものに変換できる。この例では、あと 1 ステップ必要である。それは上記の操作の 2 つを組み合わせ、第２式から第３式の 2 倍を引いたものを第３式とすることである。これによって対角要素より下の要素は全てゼロとなる。この方程式系はごく簡単にとける。我々の例は

x_{1} + 3x_{2} − 4x_{3} = 8
2x_{2} − 2x_{3} = 6
−4x_{3} = −8

となる。第 3 式は直ぐに解ける： x_{3} = 2 となる。この値を第 2 式に代入すると x_{2} が得られる。

2x_{2} − 4 = 6
x_{2} = 5

同様にこれら 2 つの値を第 1 式に代入すると、

x_{1} − 15 − 8 = 8
x_{1} = 1

から x_{1} が得られ方程式の解が得られる。この例は Gauss の消去法の 2 つの基本的段階を示している。第 1 段階は前進消去法段階で、系統的に変数を消去していき、対角要素より下の要素が全てゼロになるように変換する。この処理は三角化と呼ばれることもある。第 2 段階は後退代入段階で、第 1 段階で得られた三角行列を使って変数の値を計算する。

2.2 プログラムへの実現の概要

一般に、 N 個の未知数をもつ N 個の方程式からなる系を解く。

a_{11}x_{1} + a_{12}x_{2} + ⋯ + a_{1N}x_{N} = b_{1}
a_{21}x_{1} + a_{22}x_{2} + ⋯ + a_{2N}x_{N} = b_{2}
⋮
a_{N1}x_{1} + a_{N2}x_{2} + ⋯ + a_{NN}x_{N} = b_{N}

これらの方程式は行列を使うと 1 つの式

[ \begin{matrix}a_{11} & a_{12} & ⋯ & a_{1N} \\ a_{21} & a_{22} & ⋯ & a_{2N} \\ ⋮ & & ... \\ a_{N1} & a_{N2} & ⋯ & a_{NN} \end{matrix} ] [ x_{1} ; x_{2} ; ⋮ ; x_{N} ] = [ b_{1} ; b_{2} ; ⋮ ; b_{N} ]

すなわち、 Ax = b のように書ける。ただし、 A は係数行列、 x は変数、 b は方程式の右辺を表す。 b の要素は A の行と一緒に扱われるので、 b を A の第 (N + 1) 列とみなし N × (N + 1) の行列を用いて両者を保持すると便利である。 2.2.1 前進消去法段階について

前進消去法段階は次のようにまとめられる。はじめに、第 1 式に適当な数を掛けたものを第 1 式以外の各式に対して加えて、第 1 式以外から第 1 変数を除去する。次に、第 2 式に適当な数をかけたものを第 3 式から N 式のそれぞれに対して加えて、第 1,2 式以外の各式から第 2 変数を消去する。次にはじめの 3 式以外のすべての式から第 3 変数を消去するというようにする。第 i 変数を第 j 式 (i + 1 ≤ j ≤ N) から消去する時は、第 i 式を a_{ji}/a_{ii} 倍したものを大 j 式から引く。この処理は次のコードのようにもっと簡潔に書ける。

for(j=i+1;j<=N;j++)
    t=A[j][i]/A[i][i];
for(k=N+1;k>=i;k--)
A[j][k]-=A[i][k]*A[j][k]/A[i][i];

このコードは 3 重のループからなっており、全体の計算時間は N^{3} に比例する。 3 番目のループは、同じ行の他の要素の計算をする前に A[j][i] を壊してしまわないようにするために、逆方向に進んでいる。

2.2.2 ピボット選択について

上記のコードがこれで完全であると考えるのは単純すぎる。 A[i][i] はゼロかもしれないし、ゼロによる割り算が発生するかもしれない。しかし、これはその外側のループで i 行を適当な（i + 1 と N の間の）行と入れ替えて、 A[i][i] がゼロでないようにすれば解決する。そのような行が見つからない時には、その行列は特異 (singular) である。すなわち、一意な解が存在しない。上のコードに第 i 列中でゼロでない要素を持つ行で下の方の行を見つけてその行と i 行を入れ替えるコードを追加しなければならない。第 i 列の対角要素より下の要素を消去するのに実際に使われる要素 A[i][i] は桁軸 (pivot) ・ピボットと呼ばれる。

実際には、第 i 列で非ゼロ要素をもつ行を探すだけでなくもう少し手間をかける方がよい。第 i 列で絶対値が最大の要素をもつ行を使うのが最もよい。その理由は行のスケーリング（定数倍）に使われるピボットが非常に小さいと、第 i 変数を第 j 式 (i + 1 ≤ j ≤ N) から消去する時に使われるスケーリング因子 A[i][i] が非常にに大きくなる。実際、それは係数 A[j][i] が小さく見えるほどにまで大きくなり、 A[j][i] に値が”丸めの誤差”のため歪んでしまうこともある。

簡単に述べると、実際の表現として普通にしようされている浮動少数点表示では絶対値が非常に違う加算や引算は正確に実行されない。小さなピボットを使うとそのようなことが発生しやすくなる。第 i 列の i + 1 行から N 行の間の最大の要素を使うとスケーリング因子は 1 以下となり、この種のエラーを防止できる。大きな要素を見つけるために第 i 列より先も探してみてはと思うかもしれないが、そのような余分なことをしなくても正確な答えが得られることが知られている。

Gauss の消去法の前進段階に対する次のコードはこの処理をそのままプログラムしたものである。 1 から N までの各 i に対して、 i 列を下に向かって調べ（i 行より下の中で）最大の要素を見つける。その要素を含む行を i 行と入れ替えた後に、前とまったく同様に i + 1 から N の式から第 i 変数を消去する。

void eliminate(){
int i,j,k,max;
double t;
for(i=1; i<=N; i++){
max=i;
for(j=i+1;j<=N;j++)
    if(abs(A[j][i]) > abs(A[max][i])) max=j;
for(k=i;k<=N+1;k++){
t=A[i][k];
  A[i][k]=A[max][k];
          A[max][k]=t;
}
for(j=i+1;j<=N;j++)
    t=A[j][i]/A[i][i];
for(k=N+1;k>=i;k--)
A[j][k]-=A[i][k]*t;
}

ピボット A[i][i] を使って第 i 式以外の全ての式から（(i + 1) 式から N 式までだけでなく）第 i 変数を消去するアルゴリズムもある。この処理は完全ピボット選択と呼ばれている。前進消去ではこの一部分えお実行するだけのなので部分ピボット選択と呼ばれる。

2.2.3 後置代入段階について

前置消去が終了すると、配列 A の対角要素より下の要素は全てゼロとなっており、後置代入段階に進む。そのコードはもっと簡単に以下のようになる。

void substitute(){
    int i,j,k;
    double t;
for(j=N;j>=1;j--){
        t=0.0;
for(k=j+1;k<=N;k++) t+=A[j][k]*x[k];
x[j]=(A[j][N+1]-t)/A[j][j];
}

手続き eliminste に続いて手続き substutute を呼ぶとサイズ N の配列 x に解が得られる。

2.3 プログラム 1 （ピボット選択なし）

以下に作成したプログラムを載せる。

#include <stdio.h>
#define MAX  30
double A[MAX][MAX+1];
double x[MAX];
int N;
void equation();
void eliminate();
void substitute();

/* main 関数 */
int main() {
equation();
    eliminate();
substitute();
return 0;
}

/* 方程式の入力と表示 */

void equation(){
    int i,j;

/* 変数の数の入力 */

    while(N<MAX){
printf("How many variables ?\n");
printf("Don't over 30 variables !\n");
    printf("Please input:");scanf("%d",&N);
    printf("\n");
    if(N>MAX){
        printf("Over 30 variables !!\n");
    printf("Please input again !\n\n");
        N = 0;
    }
else{
        break;
    }
 }

/* 配列の初期化 */

    for(i=1;i<=N;i++){
for(j=1;j<=N+1;j++){

A[i][j]=0.0;
}
}
    for(i=1;i<=N+1;i++){
x[i]=0.0;
}

/* 方程式の入力 */

printf("Make equations as many as number of variables !\n");
for(i=1; i<=N; i++){
for(j=1;j<=N;j++){
            printf("A[%d%d] = ", i, j );
            scanf("%lf",&A[i][j]);
}
        printf("B[%d] = ", i);
        scanf("%lf", &A[i][N+1]);
}
    printf("\n");

/* 方程式の出力 */

    printf("****** These are %d equations ! ******\n",N);
for(i=1; i<=N; i++){
 printf("%d equation :",i);
     for(j=1;j<=N+1;j++){
 printf("%10f",A[i][j]);
     if(j==N) printf(" x%d = ",j);
 else if(j==N+1)  printf("\n");
         else  printf(" x%d + ",j);
 }
}
printf("\n");
}

/* 前進消去 */
void eliminate(){
int i,j,k;

/* 前進消去の計算 */

for(i=1; i<=N; i++){
for(j=i+1;j<=N;j++)
for(k=N+1;k>=i;k--)
A[j][k]-=A[i][k]*A[j][i]/A[i][i];
}

/* 前進消去の結果出力 */

    printf("****** Gaussian elimination ******\n");
for(i=1;i<=N;i++){
    printf("| ");
for(j=1;j<=N;j++){
printf("%10f ",A[i][j]);
}

printf(" | | x%d |",i);
        printf(" | %10f |",A[i][N+1]);
printf("\n");
}
printf("\n");
}

/* 後退代入 */
void substitute(){
    int i,j,k;
    double t;

/* 後退代入の計算 */

for(j=N;j>=1;j--){
        t=0.0;
for(k=j+1;k<=N;k++) t+=A[j][k]*x[k];
x[j]=(A[j][N+1]-t)/A[j][j];
}
    /* 解の出力 */
    printf("****** The Solution *******\n");
for(i=1;i<=N;i++){
        printf("x%d = %f \n",i,x[i]);
}
}

2.4 プログラム 2 （部分ピボット選択あり）

以下に作成したプログラムを載せる。

#include <stdio.h>
#define MAX  30
double A[MAX][MAX+1];
double x[MAX];
int N;
void equation();
void eliminate();
double abs(double x);
void substitute();

/* main 関数 */
int main() {
equation();
    eliminate();
substitute();
return 0;

}

/* 方程式の入力と表示 */

void equation(){
    int i,j;

/* 変数の数の入力 */

    while(N<MAX){
printf("How many variables ?\n");
printf("Don't over 30 variables !\n");
    printf("Please input:");scanf("%d",&N);
    printf("\n");
    if(N>MAX){
        printf("Over 30 variables !!\n");
    printf("Please input again !\n\n");
        N = 0;
    }
else{
        break;
    }
 }

/* 配列の初期化 */

    for(i=1;i<=N;i++){
for(j=1;j<=N+1;j++){
A[i][j]=0.0;
}
}
    for(i=1;i<=N+1;i++){
x[i]=0.0;
}

/* 方程式の入力 */

printf("Make equations as many as number of variables !\n");
for(i=1; i<=N; i++){
for(j=1;j<=N;j++){
            printf("A[%d%d] = ", i, j );
            scanf("%lf",&A[i][j]);
}
        printf("B[%d] = ", i);
        scanf("%lf", &A[i][N+1]);
}
    printf("\n");

/* 方程式の出力 */

    printf("****** These are %d equations ! ******\n",N);
for(i=1; i<=N; i++){
 printf("%d equation :",i);
     for(j=1;j<=N+1;j++){
 printf("%10f",A[i][j]);
     if(j==N) printf(" x%d = ",j);
 else if(j==N+1)  printf("\n");

         else  printf(" x%d + ",j);
 }
}
printf("\n");
}

/* 前進消去 */
void eliminate(){
int i,j,k,max;
double t;

/* 前進消去の計算 */

for(i=1; i<=N; i++){
max=i;
for(j=i+1;j<=N;j++)
    if(abs(A[j][i]) > abs(A[max][i])) max=j;
for(k=i;k<=N+1;k++){
t=A[i][k];
  A[i][k]=A[max][k];
          A[max][k]=t;
}
for(j=i+1;j<=N;j++)
for(k=N+1;k>=i;k--)
A[j][k]-=A[i][k]*A[j][i]/A[i][i];
}

/* 前進消去の結果出力 */

    printf("****** Gaussian elimination ******\n");
for(i=1;i<=N;i++){
    printf("| ");
for(j=1;j<=N;j++){
printf("%10f ",A[i][j]);
}
printf(" | | x%d |",i);
        printf(" | %10f |",A[i][N+1]);
printf("\n");
}
printf("\n");
}

/* 絶対値の計算*/

double abs(double x){
double y;
if(x >= 0.0) y=x;
else y=-x;
return y;
}

/* 後退代入 */
void substitute(){
    int i,j,k;
    double t;

/* 後退代入の計算 */

for(j=N;j>=1;j--){
        t=0.0;
for(k=j+1;k<=N;k++) t+=A[j][k]*x[k];
x[j]=(A[j][N+1]-t)/A[j][j];
}
    /* 解の出力 */
    printf("****** The Solution *******\n");
for(i=1;i<=N;i++){
        printf("x%d = %f \n",i,x[i]);
}
}

2.5 プログラムの解説

ここではプログラムの個々の作業について解説していく。

2.5.1 プログラムの最初の部分について

まず、配列の上限を決めておく。ここでは、「＃ define MAX 30」として、全ての関数の中で定数として MAX が使えるようにした。また、連立方程式の配列も全ての関数で使えるように、 2 次元配列 A[MAX][MAX+1] を定義し、また、解の配列 x[MAX] を両方とも double 型で宣言した。そして、定義した 4 つの関数について、それぞれ呼び出される前に宣言をしておかないとならないので、最初に宣言している。また、連立方程式の変数の数を記憶しておくための変数として、 N を定義しておく。

2.5.2 main 関数について

まず、関数 equation を呼び出して、連立方程式の変数の数を入力し、 2 次元配列 A[MAX][MAX+1]、配列 x[MAX] を初期化し、連立方程式の入出力を行なう。次に、関数 eliminate を呼び出して、ガウスの消去法の前進消去段階を行い、その結果を出力する。最後に、関数 substitute を呼び出して、後退代入段階の計算を行い、解を出力する。 2.5.3 連立方程式の入力と表示について

まず、最初に連立方程式の変数の数の入力を行なう。ここで、入力された変数の数が最初に設定した 30(MAX) を超えていたら、エラーメッセージを出力して、もう 1 度、入力させるようにしている。そのために、 while(N<MAX) として、その中で if(N>MAX) とすることで、変数の数が 30 を超えていたら N=0 にして、 while 文でまわるようにしてある。そして、もし値が 30 以下なら、 break でその while 文を抜けるようにしてある。

次に、配列の初期化であるが、 N が入力されたので、 A[1][1]〜A[N][N+1] までの配列と x[1]〜x[N] を初期化してそれぞれの値を 0 にする。ここで、配列 A[1][1]〜A[1][N] が、第 1 方程式の左辺のそれぞれの係数値であり、 A[1][N+1] が第 1 方程式の右辺となっている。また、 x[1]〜x[N] は、求める解である。

そして、連立方程式の入力を行なう。ここで、連立方程式の変数の数 N が入力されたので、「Make equations as many as number of variables !」と表示させて、変数の数と同じだけ方程式を作るように示してある。まず、一番外側の for 文で A[i][j] の i を増やしている。そして、その内側の for 文で A[i][j] の j を増やしている。つまり、第 1 方程式の係数の入力を行なってから、第 2 方程式の係数の入力を行なうようにしている。入力は「 scanf(”％ lf”, ＆ A[i][j]);」で、入力された値を配列にいれている。それぞれ、入力された値を配列にしまい、画面にはどの配列にしまわれたのかがわかるように、その配列を表示させている。そして、 j が N までいったら（つまり方程式の左辺の係数の入力が終わったら）、内側の for 文を抜けて、方程式の右辺の入力を行なう。ここで、「printf(”B[％ d] = ”, i);」として、表示としては、第 1 方程式の右辺だったら、「B[1]= 」というようにしている。これは、実際に入力された値がしまわれるのは、 A[i][N+1] であるが、表示としては、 B[i] とするほうが、「右辺」ということがわかりやすくなるので、そのようにしてある。

そして、入力された連立方程式を表示する。入力されたときは、配列として表示されていたので、さらに見やすくするために、きちんと連立方程式の形にして出力させることにした。まず、 A[i][j] を表示し、 j¡N なら「x(j)+」を表示する。これを内側の for 文で j を増やしながら、繰り返す。そして、 j==N となったら、「x(j)=」とし、その後に A[i][N+1] を表示させて右辺とする。そして、今度は j=N+1 となったら、改行を行なう。これを外側の for 文で i を増やしながら、繰り返し、 N 個の方程式を表示させている。 2.5.4 前進消去段階について

一番外側の for 文では、 i を増やしている。つまり、定数倍されて引く式を for 文でまわして決めている。また、部分ピボット選択を行なっているプログラム 2 においては、 i 列を下に向かって調べ、最大の要素を見つけたらそれを含む方程式を、その時点での第 i 列の方程式とする。そして、その入れ替えその次の for 文で行なう。以上の操作を次のようにして行なっている。

max=i;
for(j=i+1;j<=N;j++)
  if(abs(A[j][i]) > abs(A[max][i])) max=j;
for(k=i;k<=N+1;k++){
  t=A[i][k];
 A[i][k]=A[max][k];
  A[max][k]=t;

前進消去の基本的な形は以下の通りになる。

for(i=1; i<=N; i++)
  for(j=i+1;j<=N;j++)
            t = A[i][j] / A[j][j];
for(k=N+1;k>=i;k--)
A[j][k]-=A[i][k]*t;

はじめに、第 1 方程式をそれ以外のそれぞれの方程式に対応して、定数倍したものを、第 1 式以外の各式にから引いて、第 1 式以外から第 1 変数を消去する。次に、第 2 方程式を同じように定数倍したものを第 3 方程式から第 N 方程式のそれぞれから引いて、第 1、 2 以外の各方程式から第 2 変数を消去する。次に、はじめの 3 式以外のすべての式から第 3 変数を消去するというようにする。つまり、このプログラムからいうと、第 i 変数を第 j 式から消去するときは、第 i 式を A[j][i]/A[i][i] 倍したものを第 j 式から引くのである。これは、例えば、 A[2][1] ÷ A[1][1] の商を t とする。この t は前述のスケーリング因子である。それぞれの第 1 方程式の要素を定数倍し、それぞれに対応する第 2 方程式から引くのである。 A[2][2]-A[1][2]*t、 A[2][3]-A[1][3]*t、というようにする。そしてその結果の値を A[2][2]、 A[2][3] にそれぞれ入れる。同様に、 A[3][1] ÷ A[1][1] の商を t とし、 A[3][2]=A[3][2]-A[1][2]*t、 A[3][3]=A[3][3]- A[1][3]*t、というようにする。つまり、一番外側の for 文で i を増やすことで、定数倍して引く方程式が第 i 方程式だということをを決める。そして、「 t = A[i][j] / A[j][j];」でスケーリング因子を定める。次にその内側の for 文では、 j を増やすことで、引かれる方程式が第 j 方程式であることを決めている。つまり、第 1 引数を消す際に、第 1 方程式を定数倍したものを、第 2 方程式から、第 3 方程式から、・・・というように、引くが、その引かれる方程式を順に決めていくのである。そして、最後にもっとも内側の for 文では、それぞれの方程式の中で、左から順に、行列でいったら、列を増やして、引き算を行なうようにしているのである。

前進消去が終わったら、その結果を出力する。前進消去の結果、上三角行列ができる。 A[i][j] として、外側で i を増やす for 文をまわして行を増やす。そして、内側で j を増やす for 文で列を増やす。まず、 j が N まで出力して左辺を表示し、その後に解となる x[i] を出力する。そして最後に A[i][N+1] で、右辺を表示する。これをそれぞれの方程式 (行) について、繰り返して、表示させる。

2.5.5 絶対値の計算について（プログラム 2 のみ）

ここでは、渡された引数が正だったらそのまま値として返し、負だったらそれに-をつけて正にしてから値として返す。そのようにして、絶対値の計算を実現している。

2.5.6 後退代入段階について

後退代入法の計算は次の通りである。

for(j=N;j>=1;j--){
    t=0.0;
for(k=j+1;k<=N;k++) t+=A[j][k]*x[k];
x[j]=(A[j][N+1]-t)/A[j][j];
}

前進消去を行なった結果、上三角行列ができている。従って、第 N 方程式では、左辺の変数は 1 つなので、右辺と割れば、変数 x[N] がでてくる。すると、 x[N] をその上の第 N-1 方程式に代入すれば、第 N-1 方程式も変数は 1 つとなるのである。このように、代入を繰り返して全ての変数を求める。従って、外側の for 文では j を N から j–としている。つまり、第 N 方程式から順に下に向かうようにしている。そして、内側の方程式では第 j 方程式において、右辺から引くものを計算し、それを t としている。例えば、第 3 方程式に変数 x[3] の求め方は、以下の通りである。

x[3]=(A[3][N+1]-(A[3][N]*x[N]+A[3][N-1]*x[N-1]+ … + A[3][4]*x[3]))/A[3][3]

従って、内側の for 文では、この右辺から引く部分を計算して、それを t としている。最後に右辺から引いたものを求める変数の係数で割ってあげればよい。

解の出力は for 文でまわして、 x[1]〜x[N] までを順に出力させていけばよい。

2.6 プログラム 1 の実行結果（ピボット選択なし）

1. 0 除算の起こらないプログラムの実行

How many variables ?
Don't over 30 variables !
Please input:3
Make equations as many as number of variables !
A[11] = 5.2
A[12] = 1.2
A[13] = 0.3
B[1] = 4.0
A[21] = 2.2
A[22] = -1.9
A[23] = -0.01
B[2] = 2.0
A[31] = -3.0
A[32] = -1.1
A[33] = 3.5
B[3] = 5.0
****** These are 3 equations ! ******
1 equation :  5.200000 x1 +   1.200000 x2 +   0.300000 x3 =   4.000000
2 equation :  2.200000 x1 +  -1.900000 x2 +  -0.010000 x3 =   2.000000
3 equation : -3.000000 x1 +  -1.100000 x2 +   3.500000 x3 =   5.000000
****** Gaussian elimination ******
|   5.200000   1.200000   0.300000  | | x1 | |   4.000000 |
|   0.000000  -2.407692  -0.136923  | | x2 | |   0.307692 |
|   0.000000   0.000000   3.696262  | | x3 | |   7.255591 |
****** The Solution *******
x1 = 0.711236
x2 = -0.239427
x3 = 1.962954
Press any key to continue

2. 0 除算の起こるプログラムの実行

How many variables ?
Don't over 30 variables !
Please input:3
Make equations as many as number of variables !
A[11] = 0
A[12] = 1

A[13] = 3
B[1] = 5
A[21] = 2
A[22] = 4
A[23] = 6
B[2] = 8
A[31] = 3
A[32] = 5
A[33] = 7
B[3] = 9
****** These are 3 equations ! ******
1 equation :  0.000000 x1 +   1.000000 x2 +   3.000000 x3 =   5.000000
2 equation :  2.000000 x1 +   4.000000 x2 +   6.000000 x3 =   8.000000
3 equation :  3.000000 x1 +   5.000000 x2 +   7.000000 x3 =   9.000000
****** Gaussian elimination ******
|   0.000000   1.000000   3.000000  | | x1 | |   5.000000 |
|  -1.#IND00  -1.#INF00  -1.#INF00  | | x2 | |  -1.#INF00 |
|  -1.#IND00  -1.#IND00  -1.#IND00  | | x3 | |  -1.#IND00 |
****** The Solution *******
x1 = -1.#IND00
x2 = -1.#IND00
x3 = -1.#IND00
Press any key to continue

2.7 考察 (ピボット選択なし)

ピボット選択をしない場合、第 1 方程式の第 1 変数の係数、つまり、 A[1][1] = 0 だった場合、 A[1][1]/A[2][1] という段階で、 0 除算が発生してしまっている。従って、第 i 列の対角要素より下の要素を消去するのに実際に使われる要素 A[i][i]、枢軸要素は 0 になってはならない。最初の実行結果では、その枢軸要素 (ピボット) は全て、 A[i][i]! = 0 だったので、前進消去はきちんと動き、近似解を求めることができた。しかし、 2 番目の実行結果では、 0 除算が起こってしまい、前進消去の結果の段階で値にエラーが生じ、解を求める結果としてもエラーとなってしまっている。これを解消するためには、方程式の順序を適当に入れ換える、変数の順序を適当に入れ換える、ことにより、 0 でない係数がピボットの位置に来るようにする。この際、丸め誤差の影響を小さくするには、なるべく絶対値の大きな係数をピボットに使う方がよい。前述したように、このような操作をピボット選択という。

2.8 プログラム 2 の実行結果（部分ピボット選択あり）

ピボット選択がない場合に 0 除算が起こってしまった例で実行した。

How many variables ?

Don't over 30 variables !
Please input:3
Make equations as many as number of variables !
A[11] = 0
A[12] = 1
A[13] = 3
B[1] = 5
A[21] = 2
A[22] = 4
A[23] = 6
B[2] = 8
A[31] = 3
A[32] = 5
A[33] = 7
B[3] = 9
****** These are 3 equations ! ******
1 equation :  0.000000 x1 +   1.000000 x2 +   3.000000 x3 =   5.000000
2 equation :  2.000000 x1 +   4.000000 x2 +   6.000000 x3 =   8.000000
3 equation :  3.000000 x1 +   5.000000 x2 +   7.000000 x3 =   9.000000
****** Gaussian elimination ******
|   3.000000   5.000000   7.000000  | | x1 | |   9.000000 |
|   0.000000   1.000000   3.000000  | | x2 | |   5.000000 |
|   0.000000   0.000000  -0.666667  | | x3 | |  -1.333333 |
****** The Solution *******
x1 = -0.000000
x2 = -1.000000
x3 = 2.000000
Press any key to continue

2.9 考察 (部分ピボット選択あり)

実行結果からもわかるように、ピボット選択がないプログラムでは 0 除算が起こってしまい、結果がエラーとなってしまった連立方程式でも、ピボット選択を使ったプログラムでは、きちんと近似解を求めることができた。

2.10 補足

2.10.1 前進消去の工夫

前進消去の計算プログラムを最初は以下のようにしていた。

   for(j=i+1;j<=N;j++)
      for(k=N+1;k>=i;k--)
         A[j][k]-=A[i][k]*A[j][i]/A[i][i]

ここで、「A[j][k]-=A[i][k]*A[j][i]/A[i][i]」の部分では、「A[j][i]/A[i][i]」が内側の for 文の中で値が変わらないのにも関わらず、同じ計算を for 文がまわるたびに行なっていた。従って、「A[j][i]/A[i][i]」を内側の for 文の外側で計算し、 t として、それを使うことで計算量を減らしている。その結果プログラムは次のようになった。

  for(j=i+1;j<=N;j++)
            t = A[i][j] / A[j][j];
for(k=N+1;k>=i;k--)
A[j][k]-=A[i][k]*t;

2.10.2 絶対値の計算（プログラム 2 のみ）

絶対値の計算を最初は標準ライブラリ関数を使って、行わせるために、「♯ include ⟨stdio.h⟩」としようとした。しかし、標準ライブラリ関数として使える、関数 abs は int 型データの絶対値を返すものであった。従って、ここで、使っている値は全て double 型であるので、標準ライブラリ関数は使わず、自分で double 型の絶対値を計算する関数 abs を作った。

3 問題 1 の問題 3-2 の解答および考察

3.1 LU 分解の概要

3.1.1 三角行列について

正方行列の一種で、主対角線の片側が全部ゼロであるような行列のことを三角行列といい、

• 非ゼロ要素が対角の左下半分にだけあるものを下三角行列

• 非ゼロ要素が対角の右上半分にだけあるものを上三角行列という。

（下三角行列の例）

[ \begin{matrix}3 & 0 & 0 \\ 7 & 2 & 0 \\ 1 & 8 & 6 \end{matrix} ]

（上三角行列の例）

[ \begin{matrix}2 & 9 & 3 \\ 0 & 5 & 1 \\ 0 & 0 & 6 \end{matrix} ]

3.1.2 LU 分解について

任意の正則行列 A は、下三角行列と上三角行列の積の形に分解できる。すなわち、

• 下三角行列： L

• 上三角行列： U

が存在して、

A = LU

と書ける。これを A の LU 分解という。厳密にいうと「あらかじめ適当な行交換をしておけば」という条件のもとでの話であるが、とにかく三角行列の積に分解できるのである。

3.1.3 クラウト法について

LU 分解の解法は、大きく分けて二通りある。一つはガウスの消去法とほぼ同じ要領で計算する方法、もう一つはクラウト (Crout) 法と呼ばれる、積和の計算を中心とするアルゴリズムである。クラウト法は以前から手回し計算機や電卓に適した解法として知られていた方法で、コンピュータでも、レジスタをうまく使えば非常に効率よく計算できることがある。

基本的な計算手順は次のとおり。

l_{i0} = a_{i0} (i = 0, 1, ⋯, n − 1)
l_{0j} = 0 (j = 1, 2, ⋯, n − 1)
u_{00} = 1
u_{0j} = a_{0j}/l_{00} (j = 1, 2, ⋯, n − 1)
u_{i0} = 0 (i = 1, 2, ⋯, n − 1)

i = 1, 2, ⋯, n − 1 の順に

l_{ii} = a_{ii} − ∑_{k=0}^{i−1} l_{ik}u_{ki}
u_{ii} = 1

j = i + 1, ⋯, n − 1 について

l_{ji} = a_{ji} − ∑_{k=0}^{i−1} l_{jk}u_{ki}
u_{ji} = 0

l_{ij} = (a_{ij} − ∑_{k=0}^{i−1} l_{jk}u_{kj}) /l_{ii}
u_{ij} = 0

ここで計算される l,u がそれぞれ L,U の要素になる。すなわち、

L = [ \begin{matrix}l_{00} & 0 & ⋯ & 0 \\ l_{10} & l_{11} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ l_{m0} & l_{m1} & ⋯ & l_{mm} \end{matrix} ]

U = [ \begin{matrix}u_{00} & u_{01} & ⋯ & u_{0m} \\ 0 & u_{11} & ⋯ & u_{1m} \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & l_{mm} \end{matrix} ]

となっている（ただし m = n − 1）。 3.1.4 上記のアルゴリズムの解説

行列 A が A = LU と分解されるとすれば、

[ \begin{matrix}a_{00} & a_{01} & ⋯ & a_{0m} \\ a_{10} & a_{11} & ⋯ & a_{1m} \\ ⋮ & ⋮ & ... & ⋮ \\ a_{m0} & a_{m1} & ⋯ & a_{mm} \end{matrix} ] = [ \begin{matrix}l_{00} & 0 & ⋯ & 0 \\ l_{10} & l_{11} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ l_{m0} & l_{m1} & ⋯ & l_{mm} \end{matrix} ]·[ \begin{matrix}u_{00} & u_{01} & ⋯ & u_{0m} \\ 0 & u_{11} & ⋯ & u_{1m} \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & l_{mm} \end{matrix} ]

（ただし m = n − 1）であるから行列の定義から、要素別に、

l_{00}u_{00} = a_{00}
l_{i0}u_{00} = a_{i0}(i = 0, 1, ⋯, m)
l_{00}u_{0j} = a_{0j}(j = 0, 1, ⋯, m)
l_{10}u_{01} + l_{11}u_{11} = a_{11}
l_{i0}u_{01} + l_{i1}u_{11} = a_{i1}(i = 0, 1, ⋯, m)
       ⋮

などの関係が成立しなければならばならない。

このような関係式が、すべての a_{ij} について一つずつできるから、全部で n^{2} 個の「方程式」がある。それに対し「未知数」は

• 0 でない l_{ij} が\frac{n(n+1)}{2} 個

• 0 でない u_{ij} が\frac{n(n+1)}{2} 個

で、合計 n(n + 1) 個あり、したがって未知数の方が n 個多い。これは「未知数の内の n 個を自由に決めてよい（n 個の自由度がある）」ということを意味する。そこで、

u_{00} = u_{11} = ⋯ = u_{mm} = 1

にとることにして、簡単な式から順に解いていくとクラウト法の計算式が得られる。

3.2 LU 分解の連立一次方程式への応用

連立一次方程式

Ax = b  (1)

を解く場合、もし係数行列 A の LU 分解

A = L · U  (2)

ができていれば、

L · Ux = b  (3)

を解けばよいわけであるから、

y = Ux  (4)

と置いて、まず方程式

Ly = b  (5)

を解き、その解を右辺にして、方程式

Ux = y  (6)

を解くことにより、もとの方程式 (1) の解 x を得ることができる。ここで、式 (5) は

[ \begin{matrix}l_{00} & 0 & ⋯ & 0 \\ l_{10} & l_{11} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ l_{m0} & l_{m1} & ⋯ & l_{mm} \end{matrix} ] · [ y_{0} ; y_{1} ; ⋮ ; y_{m} ] = [ b_{0} ; b_{1} ; ⋮ ; b_{m} ]

また式 (6) は

[ \begin{matrix}u_{00} & u_{01} & ⋯ & u_{0m} \\ 0 & u_{11} & ⋯ & u_{1m} \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & l_{mm} \end{matrix} ] · [ y_{0} ; y_{1} ; ⋮ ; y_{m} ] = [ b_{0} ; b_{1} ; ⋮ ; b_{m} ]

の形をしているので、容易に解ける（ただし、 m = n − 1）。

すなわち、これの右辺の処理のアルゴリズムは以下のようになる。

y_{0} = b_{0}/l_{00}

i = 1, 2, ⋯, m の順に、

y_{i} = (b_{i} − ∑_{j=0}^{i−1} l_{ij}y_{j}) /l_{ii}

x_{m} = y_{m}/u_{mm} i = m − 1, m − 2, ⋯, 0 の順に

x_{i} = (y_{i} − ∑_{j=i+1}^{m} u_{ij}x_{j}) /u_{ii}

また、最後の式は、一応 u_{ii} で割っているが、前記のように

u_{ii} = 1

としてあれば、「/u_{ii}」を省略することができる。

3.3 プログラム

以下に作成した作成したプログラムを載せる。

#include <stdio.h>
#define MAX  30
/*A=LU の分解を行う関数 */
decomp(a,l,u,n)
double   a[MAX+1][MAX+1];
double   l[MAX+1][MAX+1];
double   u[MAX+1][MAX+1];

int n; /*変数の個数 */

{
    int      i,j,k;
    double  p,q,s;
    for (j=1 ; j<=n ; ++j){

l[1][j]=0.0; /*L と U の１行目の成分決

定 */
       u[1][j]=a[1][j];

}

    p=a[1][1];
    l[1][1]=1.0;

    for(i=2 ; i<=n ; ++i){

l[i][1]=a[i][1]/p; /*L と U の１列目の成分決

定 */
       u[i][1]=0.0;

}

    for(i=2 ; i<=n ; ++i){         /*L と U の対角成分の決定 */
       s=a[i][i];
       for(k=1 ; k<j ; ++k){
          s-=l[i][k]*u[k][i];

}

       l[i][i]=1.0;
       u[i][i]=q=s;
       for(j=i+1 ; j<=n ; ++j){
          s=a[i][j];
          for(k=1 ; k<j ; ++k){    /*L と U の 2,3, ・・行目の

成分決定*/

             s-=l[i][k]*u[k][j];

}

          l[i][j]=0.0;
          u[i][j]=s;
          s=a[j][i];
          for(k=1 ; k<j ; ++k){    /*L と U の 2,3, ・・列目の

成分決定 */

             s-=l[j][k]*u[k][i];

}

          l[j][i]=s/q;
          u[j][i]=0.0;

} }

/* LU 分解の出力 */

     printf("**** 下三角行列‘Ｌ’ ****\n");
     for(i=1 ; i<=n ; i++){
        for(j=1 ; j<=n ; j++){
           printf("%10f ",l[i][j]);

}

        printf("\n");

}

     printf("\n");
     printf("**** 上三角行列‘Ｕ’ ****\n");
     for(i=1 ; i<=n ; i++){
        for(j=1 ; j<=n ; j++){
           printf("%10f ",u[i][j]);

}

        printf("\n");

}

     printf("\n");
}

/*前進代入、後退代入による解の算出 */

solve(l,u,x,b,n)
double  l[MAX+1][MAX+1];

double  u[MAX+1][MAX+1];
double  x[MAX+1],b[MAX+1];
int    n;
{
     int   i,j,m;
     double  s,y[MAX+1];
     y[1]=b[1]/l[1][1];
     for(i=2 ; i<=n ; ++i){      /* 前進代入 */

s=b[i];

        for(j=1 ; j<i ; ++j){
           s-=l[i][j]*y[j];

}

        y[i]=s/l[i][i];

}

m=n;

     x[m]=y[m]/u[m][m];
     for(i=m-1 ; i>=1 ; --i){     /* 後退代入 */

s=y[i];

        for(j=i+1 ; j<=n ; ++j){
           s-=u[i][j]*x[j];

}

        x[i]=s/u[i][i];

}

/* 解の出力 */

     printf("****** The Solution *******\n");
          for(i=1;i<=n;i++){
         printf("x%d = %f \n",i,x[i]);

}

}
/* main 関数 */
int main(){
double   a[MAX+1][MAX+1],x[MAX+1],b[MAX+1];
double   l[MAX+1][MAX+1],u[MAX+1][MAX+1];
int      n,i,j;

/* 変数の数の入力 */

     while(n<MAX){

printf("変数の個数を入力してください。 \n");

       printf("30 個が限度です。 \n");
       printf("どうぞ:");scanf("%d",&n);
       printf("\n");
         if(n>MAX){

printf("変数が多すぎます。 \n");

         printf("もう 1 度どうぞ。 \n\n");

n = 0; }else{ break; } } /* 配列の初期化 */

     for(i=1 ; i<=n ; i++){
        for(j=1 ; j<=n ; j++){
           a[i][j]=0.0;
           l[i][j]=0.0;
           u[i][j]=0.0;

} }

     for(i=1 ; i<=n ; i++){
        x[i]=0.0;
        b[i]=0.0;

}

/* 左辺の入力 */ printf("左辺の係数を入力してください\n");

     for(i=1 ; i<=n ; i++){
        for(j=1 ; j<=n ; j++){
           printf("a[%d%d] = ", i, j );
           scanf("%lf",&a[i][j]);

} }

     printf("\n");

/* 右辺の入力 */ printf("右辺の値を入力してください\n");

     for(i=1 ; i<=n ; i++){
              printf("b[%d] = ", i );
              scanf("%lf",&b[i]);

}

     printf("\n");

decomp(a,l,u,n); /*左辺処理の呼び出し */

solve(l,u,x,b,n); /*前進、後退代入の呼び

出し*/
     return 0;
}

3.4 プログラムの解説

ここでは、プログラムの個々の作業について解説していく。

3.4.1 プログラムの最初の部分について

まず、配列の上限を決めておく。ここでは、「＃ define MAX 30」として、全ての関数の中で定数として MAX が使えるようにした。 3.4.2 LU 分解を行う関数 decomp について

まず、 LU 分解における 2 次元配列 a[MAX+1][MAX+1],l[MAX+1][MAX+1],u[MAX+1][MAX+1] をそれぞれ double 型で定義した。ここで、配列の要素の個数が MAX + 1 となっているのは、行列の行数・引数をそれぞれ 1, 2, ⋯, MAX としたいために、配列の要素の引数もそれに対応して 1, 2, ⋯, MAX とさせるためである。まず、

for (j=1 ; j<=n ; ++j){
       l[1][j]=0.0;
       u[1][j]=a[1][j];

では、 L の 1 行目をすべて 0 に設定している。これにより行列の乗算を考ると、 U の 1 行目は A の 1 行目と等しくなる。次に、

p=a[1][1];
    l[1][1]=1.0;
    for(i=2 ; i<=n ; ++i){
       l[i][1]=a[i][1]/p;
       u[i][1]=0.0;

}

では、 U の 1 列目をすべて 0 に設定している。また、 l_{11} は 1 としている。ここで、 a_{21} = l_{21}u_{21} より、 l_{21} = \frac{a_{21}}{u_{11}} となる。しかし、 U の 1 行目は A の 1

行目と等しくなっているので、 l_{21} = \frac{a_{21}}{p} （a_{11} = p より）となる。これを for

文で繰り返し行い、 L の 1 列目が求まる。次に

for(i=2 ; i<=n ; ++i){         /*L と U の対角成分の決定 */
       s=a[i][i];
       for(k=1 ; k<j ; ++k){
          s-=l[i][k]*u[k][i];

}

       l[i][i]=1.0;
       u[i][i]=q=s;
       for(j=i+1 ; j<=n ; ++j){
          s=a[i][j];
          for(k=1 ; k<j ; ++k){    /*L と U の 2,3, ・・行目の

成分決定*/

             s-=l[i][k]*u[k][j];

}

          l[i][j]=0.0;
          u[i][j]=s;
          s=a[j][i];
          for(k=1 ; k<j ; ++k){    /*L と U の 2,3, ・・列目の

成分決定 */

             s-=l[j][k]*u[k][i];

}

          l[j][i]=s/q;
          u[j][i]=0.0;

} }

では、まず U の対角成分を求める。ここで、 a_{22} = l_{21}u_{12} + l_{22}u_{22} より、 l_{22}u_{22} = a_{22} − l_{21}u_{12} となる。プログラムでは、まず a_{22} の値を変数 s に入れておいて、上式のような減算をした結果を再び s としている。ここで、 l_{22} を 1 とすることで u_{22} の値を q(q=s) としている。つぎに、 u_{23}, l_{32} を求める。 a_{23} = l_{21}u_{13} + l_{22}u_{23} より、 u_{23} = a_{23} − l_{21}u_{13}(l_{22} = 1 より) となる。プログラムでは、まず、 a_{23} の値を s に入れておいて、上式のような減算をした結果を再び s とし、その値を u_{23} の値としている。また、 l_{23} の値を 0 としている。

次に、 a_{32} = l_{31}u_{12} + l_{32}u_{22} より、 l_{32}u_{22} = a_{32} − l_{31}u_{12} となる。プログラムでは、まず、 a_{23} の値を s に入れておいて、上式のような減算をした結果を再び s としている。この s をさきほど求めた u_{22} の値である q で割ることで l_{32} の値を求めている。また、 u_{32} の値を 0 としている。そして、この作業を繰り返し U の 2 行目、 L の 2 列目を求める。そして次に U の対角成文である u_{33} と L の対角成文である l_{33} を求め、同じような作業を繰り返し、 U の 3 行目、 L の 3 列目を求める。

つまり、 u_{ii} と l_{ii}(i = 2, 3, ⋯, n) の値を求め、 U の i 行目・ i 列目、 L の i 行目・ i 列目という風に求めていっているのである。こうして、 LU 分解に対応する行列 L と U が、 2 次元配列 l と u に収められるのである。

そして最後に、その 2 次元配列 l と u をそれぞれ「下三角行列’L’」「上三角行列’U’」という文字とともに出力させている。

3.4.3 連立一次方程式の解の算出を行う関数 solve について

ここでは、前述の「LU 分解の連立一次方程式への応用」で述べた作業を実現している。まず、プログラムにおいて、

y[1]=b[1]/l[1][1];
     for(i=2 ; i<=n ; ++i){      /* 前進代入 */

s=b[i];

        for(j=1 ; j<i ; ++j){
           s-=l[i][j]*y[j];

}

        y[i]=s/l[i][i];

}

が LU 分解の連立一次方程式への応用で示した式 (5) の、

Ly = b

に対応しており、ここで y の値を求めている。配列 l は LU 分解による下三角行列に対応しているので、 y_{1} の値は、 b_{1}/l_{11} として簡単に求められる。次に y_{2} の値というのは、 l_{21}y_{1} + l_{22}y_{2} = b_{2} より、 y_{2} = (b_{2} − l_{21}y_{1})/l_{22} として求められるので、この値を s として、 y_{2} に代入している。次に y_{3} の値というのは、 l_{31}y_{1} + l_{32}y_{2} + l_{33}y_{3} = b_{3} より、 y3 = (b_{3} − l_{31}y_{1} − l_{32}y_{2})/l_{33} として求められるので、この値を s として、 y_{3} に代入している。これを繰り返し、 y_{i}(i = 1, 2, ⋯, n) の値を求めているのである。

次に、プログラムにおいて

m=n;
     x[m]=y[m]/u[m][m];
     for(i=m-1 ; i>=1 ; --i){     /* 後退代入 */

s=y[i];

        for(j=i+1 ; j<=n ; ++j){
           s-=u[i][j]*x[j];

}

        x[i]=s/u[i][i];

}

は LU 分解の連立一次方程式への応用で示した式 (6) の、

Ux = y

に対応しており、ここで求めたい解である x の値を求めている。まず、前記のプログラムとの区別を明確にするために、 n の値を m に代入して、これを使っている。次に、配列 u は LU 分解による上三角行列に対応しているので、 x_{m} の値は、 y_{m}/u_{mm} として簡単に求められる。次に x_{m−1} の値というのは、 u_{m−1m−1}x_{m−1} + u_{m−1m}x_{m} = y_{m−1} より、 x_{m−1} = (y_{m−1} − u_{m−1m}x_{m})/u_{m−1m−1} として求められるので、この値を s として、 x_{m−1} に代入している。次に y_{m−2} の値というのは、 u_{m−2m−2}x_{m−2}+u_{m−2m−1}y_{m−1}+ u_{m−2m}x_{m} = y_{m−2} より、 x_{m−2} = (y_{m−2}−u_{m−2m−1}y_{m−1}+u_{m−2m}x_{m})/u_{m−2m−2} として求められるので、この値を s として、 x_{m−2} に代入している。これを繰り返し、 x_{i}(i = m, m − 1, ⋯, 1) の値を求めているのである。

そして最後に解を収めた配列 x を for 文で回して出力させている。

3.4.4 main 関数について

まず、最初に連立方程式の変数の数の入力を行なう。ここで、入力された変数の数が最初に設定した 30(MAX) を超えていたら、エラーメッセージを出力して、もう 1 度、入力させるようにしている。そのために、 while(n<MAX) として、その中で if(n>MAX) とすることで、変数の数が 30 を超えていたら N=0 にして、 while 文でまわるようにしてある。そして、もし値が 30 以下なら、 break でその while 文を抜けるようにしてある。

次に、配列の初期化であるが、 n が入力されたので、配列 a[1][1]〜a[n][n],l[1][1] 〜l[n][n],u[1][1]〜u[n][n] と配列 x[1]〜x[n],b[1]〜b[n] を初期化してそれぞれの値を 0 にする。ここで、配列 a[1][1]〜a[1][n] が、第 1 方程式の左辺のそれぞれの係数値であり、配列 l[1][1]〜l[n][n],u[1][1]〜u[n][n] が配列 a[1][1]〜a[1][n] に対応する行列 A を LU 分解したときの行列 L,U に対応する配列である。そして配列 b[1]〜[n] が第 1 方程式の右辺となっている。また、配列 x[1]〜x[n] は、求める解である。

そして、連立方程式の入力を行なう。ここで、連立方程式の変数の数 n が入力されたので、「Make equations as many as number of variables !」と表示させて、変数の数と同じだけ方程式を作るように示してある。まず、一番外側の for 文で a[i][j] の i を増やしている。そして、その内側の for 文で a[i][j] の j を増やしている。つまり、第 1 方程式の係数の入力を行なってから、第 2 方程式の係数の入力を行なうようにしている。入力は「 scanf(”％ lf”, ＆ a[i][j]);」で、入力された値を配列にいれている。次に、この作業同様に、連立方程式の右辺も配列 b[i] にいれている。

そして、 LU 分解を行う関数 decomp と、連立一次方程式の解の算出を行う関数 solve を今定義してきた配列を引数として呼び出してそれぞれの関数に作業を委ね、解を出力させている。

3.5 実行結果

実行結果は、問題 1 の問題 3-4,3-5 に対応させるため、その問題における配列を用いる。

ここで、問題 3-4,3-5 における行列 A とは以下のものである。

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

3.5.1 問題 3-5-1 に対応する実行結果

ここで、問題 3-5-1 は

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

として、 Ax = b を解けというもので、その解は、

x = [ −\frac{37}{24}_{\frac{5}{6}}^{\frac{19}{12}} ]

である。そして、プログラムにおける実行結果は以下のようになった。

How many variables ?
Don't over 30 variables !
Please input:3
Make equations as many as number of variables !
A[11] = 1
A[12] = 2
A[13] = 3
A[21] = 2
A[22] = 8
A[23] = 11
A[31] = 3
A[32] = 22
A[33] = 35
B[1] = 1
B[2] = 0
B[3] = 0
****** These are 3 equations ! ******
1 equation :  1.000000 x1 +   2.000000 x2 +   3.000000 x3 =   1.000000
2 equation :  2.000000 x1 +   8.000000 x2 +  11.000000 x3 =   0.000000
3 equation :  3.000000 x1 +  22.000000 x2 +  35.000000 x3 =   0.000000
**** under-triangle matrix 'L' ****
  1.000000   0.000000   0.000000
  2.000000   1.000000   0.000000
  3.000000   4.000000   1.000000
**** upper-triangle matrix 'U' ****
  1.000000   2.000000   3.000000
  0.000000   4.000000   5.000000
  0.000000   0.000000   6.000000
****** The Solution *******
x1 = 1.583333
x2 = -1.541667
x3 = 0.833333
Press any key to continue

これは、前述の解における近似解と一致している。ゆえに、正しい実行結果が得られたことが判定できた。 3.5.2 問題 3-5-2 に対応する実行結果

ここで、問題 3-5-2 は

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

として、

Ax = b

を解けというもので、その解は、

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

である。そして、プログラムにおける実行結果は以下のようになった。

How many variables ?
Don't over 30 variables !
Please input:3
Make equations as many as number of variables !
A[11] = 1
A[12] = 2
A[13] = 3
A[21] = 2
A[22] = 8
A[23] = 11
A[31] = 3
A[32] = 22
A[33] = 35
B[1] = 1
B[2] = 2
B[3] = 3
****** These are 3 equations ! ******
1 equation :  1.000000 x1 +   2.000000 x2 +   3.000000 x3 =   1.000000
2 equation :  2.000000 x1 +   8.000000 x2 +  11.000000 x3 =   2.000000
3 equation :  3.000000 x1 +  22.000000 x2 +  35.000000 x3 =   3.000000
**** under-triangle matrix 'L' ****
  1.000000   0.000000   0.000000
  2.000000   1.000000   0.000000
  3.000000   4.000000   1.000000
**** upper-triangle matrix 'U' ****
  1.000000   2.000000   3.000000
  0.000000   4.000000   5.000000
  0.000000   0.000000   6.000000
****** The Solution *******
x1 = 1.000000

x2 = 0.000000
x3 = 0.000000
Press any key to continue

これは、前述の解における近似解と一致している。ゆえに、正しい実行結果が得られたことが判定できた。

3.5.3 問題 3-5-3 に対応する実行結果

ここで、問題 3-5-3 は

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

として、

Ax = b

を解けというもので、その解は、

x = [ −\frac{19}{12} ; \frac{37}{24} ; −\frac{5}{6} ]

である。そして、プログラムにおける実行結果は以下のようになった。

How many variables ?
Don't over 30 variables !
Please input:3
Make equations as many as number of variables !
A[11] = 1
A[12] = 2
A[13] = 3
A[21] = 2
A[22] = 8
A[23] = 11
A[31] = 3
A[32] = 22
A[33] = 35
B[1] = -1
B[2] = 0
B[3] = 1
****** These are 3 equations ! ******
1 equation :  1.000000 x1 +   2.000000 x2 +   3.000000 x3 =  -1.000000
2 equation :  2.000000 x1 +   8.000000 x2 +  11.000000 x3 =   0.000000
3 equation :  3.000000 x1 +  22.000000 x2 +  35.000000 x3 =   1.000000
**** under-triangle matrix 'L' ****
  1.000000   0.000000   0.000000
  2.000000   1.000000   0.000000

  3.000000   4.000000   1.000000
**** upper-triangle matrix 'U' ****
  1.000000   2.000000   3.000000
  0.000000   4.000000   5.000000
  0.000000   0.000000   6.000000
****** The Solution *******
x1 = -1.666667
x2 = 1.333333
x3 = -0.666667
Press any key to continue

これは、前述の解における近似解と一致している。ゆえに、正しい実行結果が得られたことが判定できた。

4 問題 1 の問題 3-3 の解答および考察

4.1 LU 分解の逆行列への応用

連立一次方程式に近似解を求めるときには、逆行列を求めるのは得策ではない。しかし、逆行列を求めることが必要な場合も生じる。ここでは、行列 A の逆行列を計算することを考える。

行列 A の逆行列を X = A^{−1} とする。ここで、

X = [ x_{1}, x_{2}, ⋯, x_{n} ]

とする。ただし、 x_{i} は逆行列 X の第 i 列を表す n 次元ベクトルとする。 e_{i} を第 i 成文が 1 でその他の成文が 0 の n 次元ベクトルとするとき、

Ax_{i} = e_{i}

であるから、 A の LU 分解を求め、上式を i = 1, 2, ⋯, n と n 回解けば逆行列が求められる事になる。

4.2 プログラム

以下に、作成したプログラムを載せる。

#include <stdio.h>
#define MAX  30

decomp(a,l,u,n)
double   a[MAX+1][MAX+1];

double   l[MAX+1][MAX+1];
double   u[MAX+1][MAX+1];
int     n;
{
int i,j,k;
double p,q,s;
for (j=1 ; j<=n ; ++j){
l[1][j]=0.0;
u[1][j]=a[1][j];
}
p=a[1][1];
l[1][1]=1.0;
for(i=2 ; i<=n ; ++i){
l[i][1]=a[i][1]/p;
u[i][1]=0.0;
}
for(i=2 ; i<=n ; ++i){
s=a[i][i];
for(k=1 ; k<j ; ++k){
s-=l[i][k]*u[k][i];
}
l[i][i]=1.0;
u[i][i]=q=s;
for(j=i+1 ; j<=n ; ++j){
s=a[i][j];
for(k=1 ; k<j ; ++k){
s-=l[i][k]*u[k][j];
}
l[i][j]=0.0;
u[i][j]=s;
s=a[j][i];
for(k=1 ; k<j ; ++k){
s-=l[j][k]*u[k][i];
}
l[j][i]=s/q;
u[j][i]=0.0;
}
}
 /* LU 分解の出力 */
printf("**** 下三角行列‘Ｌ’ ****\n");
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
printf("%10f ",l[i][j]);
}
printf("\n");
}
printf("\n");
printf("**** 上三角行列‘Ｕ’ ****\n");
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){

printf("%10f ",u[i][j]);
}
printf("\n");
}
printf("\n");
}
solve(l,u,x,b,n)
double  l[MAX+1][MAX+1];
double  u[MAX+1][MAX+1];
double  x[MAX+1],b[MAX+1];
int n;
{
int i,j,m;
double s,y[MAX];
y[1]=b[1]/l[1][1];
for(i=2 ; i<=n ; ++i){
s=b[i];
for(j=1 ; j<i ; ++j){
s-=l[i][j]*y[j];
}
y[i]=s/l[i][i];
}
m=n;
x[m]=y[m]/u[m][m];
for(i=m-1 ; i>=1 ; --i){
s=y[i];
for(j=i+1 ; j<=n ; ++j){
s-=u[i][j]*x[j];
}
x[i]=s/u[i][i];
}
}

int main(){
double   a[MAX+1][MAX+1],x[MAX+1],b[MAX+1],A[MAX+1][MAX+1];
double   l[MAX+1][MAX+1],u[MAX+1][MAX+1];
int      n,i,j,k;
while(n<MAX){

printf("何行何列（正則）の逆行列を求めたいのか入力してください。 \n");

printf("30 が限度です。 \n");
    printf("どうぞ:");scanf("%d",&n);
    printf("\n");
    if(n>MAX){
      printf("数が多すぎます。 \n");
  printf("もう 1 度どうぞ。 \n\n");
      n = 0;
    }else{
      break;
    }

  }

    /* 初期化 */
    for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
a[i][j]=0.0;
l[i][j]=0.0;
u[i][j]=0.0;
A[i][j]=0.0;
}
}
for(i=1 ; i<=n ; i++){
x[i]=0.0;
b[i]=0.0;
}

/* 行列の入力 */ printf("逆行列を求めたい行列を入力してください\n");

for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
            printf("a[%d%d] = ", i, j );
            scanf("%lf",&a[i][j]);
}
}
    printf("\n");

    decomp(a,l,u,n);
for(k=1 ; k<=n ; k++){
for(i=1 ; i<=n ; i++){
if(i==k) b[i]=1;
else     b[i]=0;
}
solve(l,u,x,b,n);
j=k;
for(i=1 ; i<=n ; i++){
A[i][j]=x[i];
}
  }

/* 逆行列の出力*/

printf("**** 逆行列‘Ａ’ ****\n");
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
printf("%10f ",A[i][j]);
}
printf("\n");

}
return 0;
 }

4.3 プログラムの解説

LU 分解を用いての逆行列を求めるこのプログラムは、前述の LU 分解を用いて連立一次方程式を解くプログラムにおいて、 main 関数を改良するだけでよい。逆行列を求めるアルゴリズム自体はすでに述べてあるので、ここでは、 main 関数を改良した部分についてだけ述べることにする。

さて、 main 関数に新たに改良を加えた部分というのは以下の部分である。

for(k=1 ; k<=n ; k++){
for(i=1 ; i<=n ; i++){
if(i==k) b[i]=1;
else     b[i]=0;
}
solve(l,u,x,b,n);
j=k;
for(i=1 ; i<=n ; i++){
A[i][j]=x[i];
}
  }

ここでは、まず最初の 2 つの for 文で、単位行列における 1 列目 2 列目 ⋯n 列目を、連立一次方程式における右辺に対応している配列 b に入れている。そして、連立一次方程式を解く関数 solve を呼び出し、その解を得る。ここで、その解というのは、右辺の配列 b に入っている行列が単位行列の i 列目 (i = 1, 2, ⋯ n) である場合、逆行列における i 列目 (i = 1, 2, ⋯ n) に対応している。故に、その得られた解を、逆行列を収める配列 A に入れていけばよい。これを次の for 文で実現している。

4.4 実行結果

ここでも、問題 3-4,3-5 の行列 A を用いて、それに対する逆行列を求めることにする。ここで、その行列 A とは、

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

であり、その逆行列 A_{−1} は計算で求めると、

A_{−1} = [ \begin{matrix}\frac{19}{12} & −\frac{1}{6} & −\frac{1}{12} \\ −\frac{37}{24} & \frac{13}{12} & −\frac{5}{24} \\ \frac{5}{6} & −\frac{2}{3} & \frac{1}{6} \end{matrix} ]

となる。そして、プログラムにおける実行結果は以下のようになった。

How many dementions ?
Don't over 30 !
Please input:3
  Input Matrix  !
A[11] = 1
A[12] = 2
A[13] = 3
A[21] = 2
A[22] = 8
A[23] = 11
A[31] = 3
A[32] = 22
A[33] = 35
****** you inputed this matrix ******
  1.000000  2.000000  3.000000  0.000000
  2.000000  8.000000 11.000000  0.000000
  3.000000 22.000000 35.000000  0.000000
***** under-triangle matrix 'L' *****
  1.000000   0.000000   0.000000
  2.000000   1.000000   0.000000
  3.000000   4.000000   1.000000
***** upper-triangle matrix 'U' *****
  1.000000   2.000000   3.000000
  0.000000   4.000000   5.000000
  0.000000   0.000000   6.000000
****** inverse matrix ******
  1.583333  -0.166667  -0.083333
 -1.541667   1.083333  -0.208333
  0.833333  -0.666667   0.166667
Press any key to continue

これは、前述の解における近似解と一致している。ゆえに、正しい実行結果が得られたことが判定できた。

5 問題 1 の問題 3-4 の解答および考察

クラウト法を用いて実際に LU 分解を試みる。結果は以下のとおり。

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

ここで、

L = [ \begin{matrix}l_{11} & l_{12} & l_{13} \\ l_{21} & l_{22} & l_{23} \\ l_{31} & l_{32} & l_{33} \end{matrix} ]

U = [ \begin{matrix}u_{11} & u_{12} & u_{13} \\ u_{21} & u_{22} & u_{23} \\ u_{31} & u_{32} & u_{33} \end{matrix} ]

とすると、

A = [ \begin{matrix}1 & 2 & 3 \\ 2 & 8 & 11 \\ 3 & 22 & 35 \end{matrix} ]
= [ \begin{matrix}l_{11} & l_{12} & l_{13} \\ l_{21} & l_{22} & l_{23} \\ l_{31} & l_{32} & l_{33} \end{matrix} ] · [ \begin{matrix}u_{11} & u_{12} & u_{13} \\ u_{21} & u_{22} & u_{23} \\ u_{31} & u_{32} & u_{33} \end{matrix} ]
= [ \begin{matrix}0 & 0 & 0 \\ l_{21} & l_{22} & l_{23} \\ l_{31} & l_{32} & l_{33} \end{matrix} ] · [ \begin{matrix}1 & 2 & 3 \\ u_{21} & u_{22} & u_{23} \\ u_{31} & u_{32} & u_{33} \end{matrix} ]
= [ \begin{matrix}1 & 0 & 0 \\ l_{21} & l_{22} & l_{23} \\ l_{31} & l_{32} & l_{33} \end{matrix} ] · [ \begin{matrix}1 & 2 & 3 \\ u_{21} & u_{22} & u_{23} \\ u_{31} & u_{32} & u_{33} \end{matrix} ]
= [ \begin{matrix}1 & 0 & 0 \\ 2 & l_{22} & l_{23} \\ 3 & l_{32} & l_{33} \end{matrix} ] · [ \begin{matrix}1 & 2 & 3 \\ 0 & u_{22} & u_{23} \\ 0 & u_{32} & u_{33} \end{matrix} ]
= [ \begin{matrix}1 & 0 & 0 \\ 2 & 1 & l_{23} \\ 3 & l_{32} & l_{33} \end{matrix} ] · [ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & u_{23} \\ 0 & u_{32} & u_{33} \end{matrix} ]
= [ \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & l_{32} & l_{33} \end{matrix} ] · [ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & u_{32} & u_{33} \end{matrix} ]

= [ \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & l_{33} \end{matrix} ] · [ \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & u_{33} \end{matrix} ]
= [ \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} ]

故に、

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

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

となり、解を得る。

6 問題 1 の問題 3-5 の解答および考察

6.1 問題 1 の問題 3-5-1 の解答および考察

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

として、

Ax = b  (7)

を解く。

ここで問題 3-4 より、行列 A はすでに LU 分解されている。故に、

A = LU  (8)

から、式 (8) を式 (7) に代入して

Ax = L(UX) = b  (9)

を得る。故に、式 (9) から、

UX = y  (10)

とおいて

Ly = b  (11)

を解き、

y = L^{−1}b  (12)

を求め、次に

Ux = y  (13)

から

x = U^{−1}y  (14)

を求めればよい。

ここで、

y = [ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ]  (15)

とすると、式 (11) より、

[ \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} ]  (16)

となり、

\begin{matrix}y_{1} = 1 \\ 2y_{1} + y_{2} = 0 \\ 3y_{1} + 4y_{2} + y_{3} = 0 \end{matrix}  (17)

が得られ、これを解いて、

\begin{matrix}y_{1} = 1 \\ y_{2} = −2 \\ y_{3} = 5 \end{matrix}  (18)

となり、

y = [ \begin{matrix}1 \\ −2 \\ 5 \end{matrix} ]  (19)

が得られる。

ここで、

x = [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ]  (20)

とすると、式 (10) と式 (19) より、

[ \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} ]  (21)

を得る。故に、

\begin{matrix}x_{1} + 2x_{2} + 3x_{3} = 1 \\ 2x_{1} + 5x_{2} = −2 \\ 6x_{3} = 5 \end{matrix}  (22)

が得られ、これを解いて、

\begin{matrix}x_{1} = \frac{19}{12} \\ x_{2} = −\frac{37}{24} \\ x_{3} = \frac{5}{6} \end{matrix}  (23)

となり、

x = [ −\frac{37}{24}_{\frac{5}{6}}^{\frac{19}{12}} ]  (24)

解を得る。

6.2 問題 1 の問題 3-5-2 の解答および考察

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

として、

Ax = b  (25)

を解く。

ここで問題 3-4 より、行列 A はすでに LU 分解されている。故に、

A = LU  (26)

から、式 (25) を式 (26) に代入して

Ax = L(UX) = b  (27)

を得る。故に、式 (27) から、

UX = y  (28)

とおいて

Ly = b  (29)

を解き、

y = L^{−1}b  (30)

を求め、次に

Ux = y  (31)

から

x = U^{−1}y  (32)

を求めればよい。

ここで、

y = [ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ]  (33)

とすると、式 (29) より、

[ \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} ]  (34)

となり、

\begin{matrix}y_{1} = 1 \\ 2y_{1} + y_{2} = 2 \\ 3y_{1} + 4y_{2} + y_{3} = 3 \end{matrix}  (35)

が得られ、これを解いて、

\begin{matrix}y_{1} = 1 \\ y_{2} = 0 \\ y_{3} = 0 \end{matrix}  (36)

となり、

y = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ]  (37)

が得られる。

ここで、

x = [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ]  (38)

とすると、式 (28) と式 (37) より、

[ \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} ]  (39)

を得る。故に、

\begin{matrix}x_{1} + 2x_{2} + 3x_{3} = 1 \\ 2x_{1} + 5x_{2} = 0 \\ 6x_{3} = 0 \end{matrix}  (40)

が得られ、これを解いて、

\begin{matrix}x_{1} = 1 \\ x_{2} = 0 \\ x_{3} = 0 \end{matrix}  (41)

となり、

x = [ \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ]  (42)

解を得る。

6.3 問題 1 の問題 3-5-3 の解答および考察

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

として、

Ax = b  (43)

を解く。

ここで問題 3-4 より、行列 A はすでに LU 分解されている。故に、

A = LU  (44)

から、式 (44) を式 (43) に代入して

Ax = L(UX) = b  (45)

を得る。故に、式 (45) から、

UX = y  (46)

とおいて

Ly = b  (47)

を解き、

y = L^{−1}b  (48)

を求め、次に

Ux = y  (49)

から

x = U^{−1}y  (50)

を求めればよい。

ここで、

y = [ \begin{matrix}y_{1} \\ y_{2} \\ y_{3} \end{matrix} ]  (51)

とすると、式 (47) より、

[ \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} ]  (52)

となり、

\begin{matrix}y_{1} = −1 \\ 2y_{1} + y_{2} = 0 \\ 3y_{1} + 4y_{2} + y_{3} = 1 \end{matrix}  (53)

が得られ、これを解いて、

\begin{matrix}y_{1} = −1 \\ y_{2} = 2 \\ y_{3} = −5 \end{matrix}  (54)

となり、

y = [ \begin{matrix}−1 \\ 2 \\ −5 \end{matrix} ]  (55)

が得られる。

ここで、

x = [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ]  (56)

とすると、式 (46) と式 (55) より、

[ \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} ]  (57)

を得る。故に、

\begin{matrix}x_{1} + 2x_{2} + 3x_{3} = −1 \\ 2x_{1} + 5x_{2} = 2 \\ 6x_{3} = −5 \end{matrix}  (58)

が得られ、これを解いて、

\begin{matrix}x_{1} = −\frac{19}{12} \\ x_{2} = \frac{37}{24} \\ x_{3} = −\frac{5}{6} \end{matrix}  (59)

となり、

x = [ −\frac{19}{12} ; \frac{37}{24} ; −\frac{5}{6} ]  (60)

解を得る。

7 問題 2 の解答および考察

7.1 性質 3 の証明

7.1.1 性質 3

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

∥RA − I∥ < 1

ならば A は正則であり、このとき

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

が成り立つ。

7.1.2 証明補題として「バナッハの摂動定理」を証明し、これを利用する。

バナッハの摂動定理 X,Y をバナッハ空間、 A ： X → Y を有界線形作用素で有界な逆作用素 A^{−1} ： Y → X をもつとする。このとき、有界線形作用素 B ： X → Y が条件

∥A^{−1}B∥_{L(X,X)} ≤ 1 を満たすとき、有界線形作用素 A + B ： X → Y は有界な逆作用素 (A + B)^{−1} ： Y → X をもち

∥(A + B)^{−1}∥_{L(X,X)} ≤ \frac{∥A^{−1}∥_{L(X,X)}}{1 − ∥A^{−1}B∥_{L(X,X)}}

となる。

補題の証明 y を Y から任意に選び固定しておく。 T ： X → X を x ∈ X に対し

T x = A^{−1}(y − Bx)

と定める。 x_{1} と x_{2} を X の任意の要素として、

∥T x_{1} − T x_{2}∥_{X} = ∥A^{−1}B(x_{1} − x_{2})∥_{X} ≤ ∥A^{−1}B∥_{L(X,X)}∥x_{1} − x_{2}∥_{X}

より、 T ： X → X が X 上の縮小写像となることが分かる。したがって、縮小写像原理より

T x^{∗} = x^{∗}

をみたす X の要素 x^{∗} がただ 1 つ定まる。この x^{∗} は

A(T x^{∗} − x^{∗}) = y − Bx^{∗} − Ax^{∗} = 0 (∵ T x = A^{−1}(y − Bx))

より (A + B)x^{∗} = y を満たすことがわかる。以上をまとめると、 y ∈ Y を選ぶと (A + B)x^{∗} = y を満たす x^{∗} ∈ X がただ 1 つ選べるから、 (A + B)^{−1} ： Y → X の存在がわかり、

x^{∗} = (A + B)^{−1}y

となる。

∥x^{∗}∥_{X} = ∥T x^{∗}∥_{X} = ∥A^{−1}y − A^{−1}Bx^{∗}∥_{X}
≤ ∥A^{−1}∥_{L(X,X)}∥y∥_{Y} + ∥A^{−1}B∥_{L(X,X)}∥x^{∗}∥_{X}

より

∥(A + B)^{−1}y∥_{X} ≤ \frac{∥A^{−1}∥_{L(X,X)}}{1 − ∥A^{−1}B∥_{L(X,X)}} ∥y∥_{Y}

を得る。これから、求める評価式を得る。（証明終）

では、これを使って本題を証明しよう。上のバナッハの摂動定理において、 A^{−1} = R, B = R^{−1} − A とすると、

∥\{A + (R^{−1} − A)\}^{−1}∥ ≤ \frac{∥R∥}{1 − ∥R(R^{−1} − A)∥}

∥(R^{−1})^{−1}∥ ≤ \frac{∥R∥}{1 − ∥I − RA∥}
∥R∥ ≤ \frac{∥R∥}{1 − ∥RA − I∥}
∥A^{−1}∥ ≤ \frac{∥R∥}{1 − ∥RA − I∥}

が得られる。

7.2 性質 5 の証明

7.2.1 性質 5

A ： n × n （複素）行列とする。このとき行列 A の固有値は、

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

の中に存在する。（ヂルシュゴリの定理）

7.2.2 証明

Ax=λx、 x ≠ 0 とする。 x の成分中、絶対値が最大のものの 1 つを x_{i} とする。 x_{i} ≠ 0 である。 Ax=λx の第 i 成分を考えると、

∑_{k=1}^{n} a_{ik}x_{k} = λx_{i}

すなわち

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

である。この絶対値をとると

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

である。 |x_{i}| > 0 だからこれで割れば、 A の固有値 λ は

|λ − a_{ii}| ≤ ∑_{k=1、 k≠i}^{n}|a_{ik}|

を満たす。

7.3 性質 6 の証明

7.3.1 性質 6

IEEE754 の下で浮動小数点演算の下で議論する。最近点での丸め (nearest) のモードで考える。

u = 2^{−54}

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

x, y ∈ F, 0 ∈ +, −, ∗, /
x ∘ y = (x ∘ y)(1 + δ)
|δ| ≤ u

となる。これを利用すると、 x, y ∈ F^{n} とすると、 s = xy(実数としての計算) に対して、

near();
t=0;
for(i=1;i<=n;i++)
    t+=xi*yi;

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

|x − t| ≤ δ_{n} ∑_{i=1}^{n} |x_{i}y_{i}|

δ_{n} = \frac{nu}{1 − nu}

7.3.2 証明

|s − t| = |(x_{1}y_{1} − x_{m1}y_{m1}) + … + (x_{n}y_{n} − x_{mn}ymn)|
≤ |x_{1}y_{1} − x_{m1}y_{m1}| + … + |x_{n}y_{n} − x_{mn}ymn|
= |x_{1}y_{1} − (1 + δ)x_{1}y_{1}| + … + |x_{n}y_{n} − (1 + δ)x_{n}y_{n}|
= |δ||x_{1}y_{1}| + … + |δ||x_{n}y_{n}|

= |δ| ∑_{i=1}^{n} |x_{i}y_{i}|

≤ u ∑_{i=1}^{n} |x_{i}y_{i}|

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

\frac{nu}{1 − nu} − u = \frac{nu − u + nu^{2}}{1 − nu}

= \frac{(n − 1)u + nu^{2}}{1 − nu}

n ≤ 1 より、

\frac{nu}{1 − nu} > u

8 感想

今回は課題の量がとてつもなく膨大なもので、すべてやり遂げるのに、徹夜の日々が続いた。内容的には、ガウスの消去法と LU 分解は比較的簡単に出来たが、性質の証明に関しては、一体何をやっているのかさえ分からず、悪戦苦闘の連続であった。また、今回は図書館で様々な文献を調べて、色々なプログラムを参考にした。しかしどれも基本となるアルゴリズムは一緒で、みんなと同じではつまらないと思い、自分で考えようとしたがなかなか思いつかなく、断念した。やはり、長年使われてきたアルゴリズムには安易に改良できる「抜け目」などないのであろう。全体を通して、教科書の問題の方は完璧に近い理解が得られたと思う。しかし証明の方は正直、不足している部分があり、完璧な理解は得られなかった。これから勉強して理解を深めたいと思う。

とにかく、今回の課題は質、量ともにとてつもないもので、他の課題も合わせ、死にそうな日々が続いた。何とか提出にこぎつけたが、しっかり復習をして、試験のために知識を抜け目のないものにしていきたい。