目次

1 教科書「数値計算」 p31 の問題 3(1〜5) を解け。 1

1.1 問題 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2 ガウスの消去法のプログラム (ピボット選択なし) . . . . . . . . . . . . . . . . . . . 1 1.1.3 プログラムの解説 (ピボット選択なし) . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.4 実行結果 (ピボット選択なし) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.5 考察 (ピボット選択なし) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.6 ガウスの消去法のプログラム (ピボット選択あり) . . . . . . . . . . . . . . . . . . . 7 1.1.7 プログラムの解説 (ピボット選択あり) . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.1.8 実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.1.9 考察 (ピボット選択あり) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.2 問題 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.1 問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.2 作成したプログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.3 プログラムの解説 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.2.4 実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.2.5 考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

1.3 問題 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.3.1 問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.3.2 作成したプログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.3.3 プログラムの解説 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.3.4 実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.3.5 考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

1.4 問題 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.4.1 問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.4.2 机上による計算結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.4.3 プログラムによる実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 1.4.4 考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

1.5 問題 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.5.1 問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.5.2 机上による計算結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.5.3 プログラムによる実行結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 1.5.4 考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

2 性質 3 と 5 と 6 を証明せよ。 (性質 1,2,4 は仮定してよい) 41

2.1 性質 3 の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.1.1 性質 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.1.2 証明 (性質 3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2.2 性質 5 の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.2.1 性質 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.2.2 証明 (性質 5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

2.3 性質 6 の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.3.1 性質 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.3.2 証明 (性質 6) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3 感想 45

数値計算（大石教授）第 3 回レポート

{G99P119-8 \atop 福田浩章}

出題日:2000-06-01 提出期限:2000-06-15 提出日:2000-06-15

1 教科書「数値計算」 p31 の問題 3(1〜5) を解け。

1.1 問題 1

1.1.1 問題

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

1.1.2 ガウスの消去法のプログラム (ピボット選択なし)

#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]);
}
}

1.1.3 プログラムの解説 (ピボット選択なし)

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

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

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 個の方程式を表示させるのである。

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

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][j];;

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

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

1.1.4 実行結果 (ピボット選択なし) 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

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

ピボット選択をしない場合、第 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 でない係数がピボットの位置に来るようにする。この際、丸め誤差の影響を小さくするには、なるべく絶対値の大きな係数をピボットに使う方がよい。このような操作をピボット選択という。以下のプログラムではピボット選択の機能が追加されている。

1.1.6 ガウスの消去法のプログラム (ピボット選択あり)

#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]);
}
}

1.1.7 プログラムの解説 (ピボット選択あり)

ピボット選択以外の方程式の入力と表示、前進消去の基本形、後退代入、については最初のプログラムと同じなので、ここでは解説を省く。ここでは、主にピボット選択の機能が追加されたので、その部分のプログラムについて解説する。まず、前進消去のプログラムを以下のように変更した。

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];
}

前進消去の基本的な形は変わっていない。次の部分で絶対値が最大のものを求めている。

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

ここで、 i 行目の次の行から j が始まり、 for 文で回し、 A[j][i] と A[max][i] の絶対値をとって、比較し、その値の大きい方の行を max としている。ここで、絶対値の計算として関数 abs を呼び出している。これは、自分で作った double 型の絶対値を求める関数であり、以下の通りである。

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

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

そして、次に、第 i 列において絶対値が最大なものが見つかったら、今度は第 i 行と、その行、つまり第 max 行のそれぞれの要素を入れ替えている。それが以下の部分である。

for(k=i;k<=N+1;k++){
t=A[i][k];
  A[i][k]=A[max][k];
          A[max][k]=t;
}

ここで、値の入れ換えをするので、 double 型の変数 t は、一時的に値を記憶するものとして利用している。例えば、第 1 列において、 3 行目第 1 係数の絶対値が最大だったら、 A[1][1] と A[3][1] を交換し、 A[1][2] と A[3][2] を交換し、・・・というように交換が進んでいくのである。つまり、第 max 行が今度は第 i 行となり、そして、先程の解説にもかいたからここでは、説明を省くが、前進消去の計算へとうつるのである。

1.1.8 実行結果

ピボット選択がない場合に 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

1.1.9 考察 (ピボット選択あり)

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

そこで、ここでは、ピボット選択についての説明をしたいと思う。例えば、次のような連立方程式がある。

0x + 1y = 1
2x + 3y = 5

この場合、ガウスの消去法を行なうと 0 による割り算が発生してしまう。このようなことが起こらないようにするためにピボット選択を行なう。ここでは方程式の順序を適当に入れ換えることにとどまらず、第 i 列で絶対値が最大の要素をもつ行を使うのが最もよい。その理由は行のスケーリング (定数倍) に使われる枢軸要素が非常に小さいと大きな計算誤差が生じるからである。 A[i][i] が非常に小さいと、第 i 変数を第 j 方程式 (i+1¡=j¡=N) から消去する時に使われるスケーリング因子 A[j][i]/A[i][i] が非常に大きくなり、 A[j][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;

}

ガウスの消去法の第 i 段において、

a_{i+1,i}, a_{i+2,i}, …, a_{n,i}

の中で絶対値最大のもの

a_{imax,i}

を探し、第 i 行と第 imax 行を交換する。これを部分選択という。部分選択ならば、絶対値最大の係数を探す範囲あ狭いので、計算時間の増加は微々だるものである。数学的には、単に方程式の順序を入れ替えているだけであるから、あとで結果の順序を入れ換える必要はない。したがってプログラムは簡単ですむ。また、完全選択というものもある。完全選択では、第 i 行より下、第 i 列より右すべての係数の中で絶対値が最大のものを探す。しかし、精度の面では有利であるが、絶対値最大の係数を探す範囲が広く、そのために計算時間が約 2 倍に増加してしまう。また、変数の順序が変わってしまうので、逆変換が必要である、などの問題がある。

1.2 問題 2

1.2.1 問題

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

1.2.2 作成したプログラム

#include <stdio.h>
#define MAX  30
double A[MAX][MAX+1],B[MAX],l[MAX][MAX],u[MAX][MAX];
double x[MAX];
int N;
void equation();
void decomp();
void solve();

/* main 関数 */
int main() {
equation();
    decomp();
solve();
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]);
}
}
for(i=1; i<=N; i++){
    printf("B[%d] = ", i);
        scanf("%lf", &B[i]);
    }
    printf("\n");

/* 方程式の出力 */

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

}

/* LU 分解の計算 */
void decomp(){
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("**** under-triangle matrix 'L' ****\n");
for(i=1 ; i<=N ; i++){
    for(j=1 ; j<=N ; j++){
printf("%10f ",l[i][j]);
    }
    printf("\n");
}
printf("\n");
printf("**** upper-triangle matrix 'U' ****\n");
for(i=1 ; i<=N ; i++){
    for(j=1 ; j<=N ; j++){
printf("%10f ",u[i][j]);
    }
    printf("\n");
}
printf("\n");
}

/* 連立一次方程式の解の計算 */

void solve(){
    int i,j;
    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];
    }
    x[N]=y[N]/u[N][N];
    for(i=N-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]);
}
}

1.2.3 プログラムの解説

1. 方程式の入力連立一次方程式の入力はさきほどのガウスの消去法のときに使用した関数 equation を用いる。従って、連立一次方程式の入力と表示に関しては説明を省く。ただし、ガウスの消去法の際は、第 i 方程式の右辺値を A[i][N+1] としていたが、今度は、右辺値を保持するための配列として B[N] を作り、左辺値を入力したあとに、右辺値は配列 B[N] に入力させるようにした。それは以下の部分である。

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

そして、ガウスの消去法のときと同様に、入力された方程式を出力させるようにした。

2. LU 分解の計算次に LU 分解の計算である。 LU 分解の計算は関数 decomp で行われており、それは main 関数の中で、関数 equation が呼び出されたあとに呼び出されている。 LU 分解では、 L か U どちらかの対角線要素を全て 1 とする。ここでは、 L の対角線要素、つまり l[1][1], l[2][2], …, l[N][N] を 1 とするようにしている。まず、最初の for 文は次の通りである。

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

ここでは、 L の第 1 行目の要素をすべて、 0 にしている。これは、 L が下三角行列であるので、第 1 行は l[1][1] を除いて、すべて 0 であるからである。また、 A[1][j] = l[1][1] ∗ u[j][1] となるが、 l[1][1] は後に代入するが 1 となるので、 A[1][j] = u[j][1] となるので、 u の第 1 行を A[1][j] として、代入している。

そして、 2 番目の for 文は次のようになっている。

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;

}

次に p = A[1][1] とし、ここで、先程も触れたが、 l[1][1] に 1.0 を代入している。ここで、 A[i][1] = l[i][1] ∗ u[1][1] + l[i][2] ∗ u[1][2] + … + l[i][N] ∗ u[1][N] となるが、 U の第 1 列要素は u[1][1] を除いて全て 0 なので、 A[i][1] = l[i][1] ∗ u[1][1] とすることができる。ここで、 u[1][1] = A[1][1] なので、 l[i][1] = A[i][1]/p(p は A[1][1]) として、の第 1 列を求めることができるのである。そして、 u[i][1] = 0.0 として、 u[1][1] 以外の U の第 1 列を全て 0 にしている。ここまでの段階で、 L と U それぞれの第 1 行目と第 1 列目は求められている。

では、いよいよ、 LU 分解法の主要な計算部分である。その部分は次の通りである。

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;
}
    }

まず、 L と U の第 1 行と第 1 列はそれぞれ、求められた。従って、一番外側の for 文は i = 2 から始まっている。ここで、 A[i][i] = l[i][1] ∗ u[1][i] + l[i][2] ∗ u[2][i] + … + l[i][i] ∗ u[i][i] となっている。これは、 L の第 i 行において、 i 列以降は 0、 U の第 i 列において、第 i 行以降は 0 となっているからである。また、 l[i][i] は 1 であるので、この式より、 u[i][i] は次のように求められる。 u[i][i] = A[i][i] − (l[i][1] ∗ u[1][i] + l[i][2] ∗ u[2][i] + … + l[i][i − 1] ∗ u[i − 1][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;

また、 l[i][i]=1.0 として、 L の対角線要素を全て 1 にしている。

次は、例えば A[2][3]=l[2][1]*u[1][2]+l[2][2]*u[2][2] となる。これは、 l[2][3] 以降が全て 0 であるためである。また、 l[2][2] は 1 であるので、 u[2][2] = A[2][3] − (l[2][1] ∗ u[1][2]) となる。よってこれを以下の for 文で実現する。

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;

ここでは、 k を増やす for 文において、 A[i][j] から引く要素を増やしている。また、 j を増やす for 文で j を増やして、 U の i 行における要素を次々に求めていくことになる。また、 l[i][j] = 0.0 で、 l の i 行における i 列以降の要素を 0 としている。

その次の for 文は次のようになっている。

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;

ここでは、今度は L の第 i 列における要素を求めている。つまり、 l[i + 1][i], l[i + 2][i], …, l[N][i] を求めているのである。そして、先程とは逆に、 u[j][i] = 0.0 で、 u の i 列における i 行以降の要素を 0 としているのである。

こうして、これらを一番外側の for 文で i を増やしながら、まわすことで、 U の第 2 列を求め、 L の第 2 列を求め、 U の第 3 列を求め、 L の第 3 列を求め、というように、左上から順に U と L が求められていくのである。

そして、 L と U の表示である。これはそれぞれ、 for 文を回して、表示させていけばよい。

3. L と U から解 x を求める。

L と U を使って解を求める。 ly = B、 ux = y として、解 x を求める。 • LY=B まず、 L の第 1 要素は L[1][1] 以外は 0 なので、 B[1] = l[1][1]∗B[1] となるので、 y[1] = B[1]/l[1][1] とする。そして、次に for 文で i = 2, 3, …, N の順に

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

としている。例えば、 B[3] = l[3][1] ∗ y[1] + l[3][2] ∗ y[2] + l[3][3] ∗ y[3] となる。ここで、 l[3][4] 以降は 0 なので、その分の計算はしていない。従って、 y[3] = (B[3]−l[3][1]∗y[1]−l[3][2]∗y[2])/l[3][3] となるのである。この割られる側が j¡i のので、それで for 文を回している。この操作は外側の for 文でまわして、 y[2]〜y[N] を求めている。

• UX=B ここでは、まず、 u は上三角行列で、 N 行目は u[N][N] 以外はすべて 0 なので、まず、 x[N] を求めることができる。 x_{N} = y_{N}/u_{NN} を行い、 x[N] が求められる。次にそれを N-1 行目に代入して今度は x[N-1] を求める。このように、解は下から求められていくので、 i = N − 1 からはじまり、 − − i と、 1 ずつ i を減らしていく。ここでも、同様に、解を求めるために、次々と代入していくのである。そして、 x[1]〜x[N] が求められる。

4. 解の出力ここでは、 for 文を回して、求められた、 x[1]〜x[N] を次々に表示させていっている。

1.2.4 実行結果

3 変数の連立方程式を入力し、解を求める。

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

-1.000000  -0.500000   1.000000
**** upper-triangle matrix 'U' ****
1.000000   3.000000  -4.000000
0.000000  -2.000000   2.000000
0.000000   0.000000   2.000000
****** The Solution *******
x1 = 1.000000
x2 = 5.000000
x3 = 2.000000

これは手計算で求めた結果とも一致した。

1.2.5 考察

1. LU 分解ここでは、 LU 分解の説明をする。正方行列の一種で、主対角線の片側が全部 0 であるような行列のことを三角行列といい、非 0 要素が対角の左下半分にだけあるものを下三角行列、非 0 要素が対角の右上半分にだけあるものを上三角行列という。例えば、次のようなものである。

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

任意の正則行列 A は、下三角行列と上三角行列の積の形に分解できる。すなわち、下三角行列 L、上三角行列 U、が存在して、 A = LU と書ける。これを A の LU 分解という。ここでは、下三角行列 L、上三角行列 U、を求めて表示するために、クラウト法を用いる。

2. クラウト法基本的な計算手順は以下のとおりである。

l_{i1} = a_{i1}(i = 1, 2, …, n)
l_{1j} = 0(j = 2, 3, …, n)
u_{11} = 1
u_{1j} = a_{1j}/l_{11}(j = 2, 3, …, n)
u_{i1} = 0(i = 2, 3, …, n)

i = 2, 3, …, n の順に

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

u_{ii} = 1

j = i + 1, …, n について

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

u_{ji} = 0

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

u_{ij} = 0

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

L = ( \begin{matrix}l_{11} & 0 & … & 0 \\ l_{21} & l_{22} & … & 0 \\ … \\ … & … & … & 0 \\ l_{n1} & l_{n2} & … & l_{nn} \end{matrix} ) U = ( \begin{matrix}u_{11} & u_{12} & … & u_{1n} \\ 0 & l_{22} & … & u_{2n} \\ … \\ 0 & … \\ 0 & 0 & … & u_{nn} \end{matrix} )

3. 連立 1 次方程式の解

連立 1 次方程式

Ax = B

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

A = LU

ができていれば、

LUx = B

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

y = Ux

と置いて、まず方程式

Ly = B

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

Ux = y

を解くことにより、もとの方程式の解 x を得ることができる。それらは以下の通りである。

( \begin{matrix}l_{11} & 0 & … & 0 \\ l_{21} & l_{22} & … & 0 \\ … \\ … & … & … & 0 \\ l_{N1} & l_{N2} & … & l_{NN} \end{matrix} ) ( \begin{matrix}y_{1} \\ y_{2} \\ … \\ … \\ y_{N} \end{matrix} ) = ( \begin{matrix}B_{1} \\ B_{2} \\ … \\ … \\ B_{N} \end{matrix} )

( \begin{matrix}u_{11} & u_{12} & … & u_{1N} \\ 0 & l_{22} & … & u_{2N} \\ … \\ 0 & … \\ 0 & 0 & … & u_{NN} \end{matrix} ) ( \begin{matrix}x_{1} \\ x_{2} \\ … \\ … \\ x_{N} \end{matrix} ) = ( \begin{matrix}y_{1} \\ By_{2} \\ … \\ … \\ y_{N} \end{matrix} )

1.3 問題 3

1.3.1 問題

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

1.3.2 作成したプログラム

#include <stdio.h>
#define MAX  30
double A[MAX][MAX],B[MAX],l[MAX][MAX],u[MAX][MAX],
       x[MAX],X[MAX][MAX];
int N;
void matrix();
void decomp();
void solve();
void inverse();

/* main 関数 */
int main() {
    int i,j,k;
    matrix();
    decomp();
    for(k=1 ; k<=N; k++){
for(i=1 ; i<=N ; i++){
            if(i==k) B[i]=1;
        else     B[i]=0;
        }
        solve();
j=k;
for(i=1 ; i<=N ; i++){
    X[i][j]=x[i];
}
    }
    inverse();
    return 0;
}

/* 行列の入力と表示 */

void matrix(){
    int i,j;

/* 次元の入力 */

    while(N<MAX){
printf("How many dementions ?\n");
printf("Don't over 30 !\n");
    printf("Please input:");scanf("%d",&N);
    printf("\n");
    if(N>MAX){
        printf("Over 30 !!\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("Input Matrix  !\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("********** you inputted this matrix **********\n");
for(i=1; i<=N; i++){
    for(j=1;j<=N;j++){
       printf("%10f",A[i][j]);
}
    printf("\n");
}
printf("\n");
}

/* LU 分解の計算 */
void decomp(){
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("********** under-triangle matrix 'L' **********\n");
for(i=1 ; i<=N ; i++){
for(j=1 ; j<=N ; j++){
printf("%10f ",l[i][j]);
}
printf("\n");
}
printf("\n");
printf("********** upper-triangle matrix 'U' **********\n");
for(i=1 ; i<=N ; i++){
for(j=1 ; j<=N ; j++){
printf("%10f",u[i][j]);
}
printf("\n");
}
printf("\n");
}

/* 連立一次方程式の解の計算 */

void solve(){
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];
    }
}

/* 逆行列の出力*/

void inverse(){
int i,j;
printf("********** inverse matrix ‘ X ’ **********\n");
    for(i=1 ; i<=N ; i++){
for(j=1 ; j<=N ; j++){
    printf("%10f ",X[i][j]);
}
printf("\n");
    }
}

1.3.3 プログラムの解説

まず、行列の入力を行う。これはガウスの消去法でも、方程式の入力の際に使った関数 equation を利用して、あらたに逆行列を求めたい行列を入力し、その行列を出力する関数 matrix を設けた。ここでの内容はほとんど関数 equation と同じなので、この部分の説明は問題 1 のプログラムを参照してもらいたい。ここでは、方程式ではないので、右辺値は入力する必要がないので、入力された値は A[N][N] にのみ入る。

行列が入力されたら、関数 decomp を呼び出して、入力された行列を LU 分解し、行列 L と行列 U を求める。関数 decomp についても、問題 2 における LU 分解のプログラムのものと全く同じものを利用しているので、ここでは説明を省略する。

そして、 L と U が求められたら、ここからが、逆行列を求める部分の主要な部分となる。まず、 AX = e_{i} となるわけだが、この単位行列 e_{i} を左から 1 列ずつつくり、そして、解 x も左から 1 列ずつもとめていくのである。まず、外側の for 文で k を増やし、その内側で i を増やしている。そして、 i と k の値が一致したら、 B[i] を 1 としている。つまり X と B の第 1 列について求めるときは B[1] = 1 となり、第 2 列について求めるときは、 B[2] = 1 となり、それ以外は 0 となるのである。そして、 N 行 1 列の行列 x と同じく N 行 1 列の行列 B について、 LU 分解の解を求める計算関数 solve を呼び出している。関数 solve についても、問題 2 で説明してあるのでそちらを参考にしてもらいたい。そして、まず、左端から解 x[1]〜x[N] が求められ、それは X[i][1] に格納される。そこまでの操作を for 文で繰り返し、次に x[1]〜x[N] が B[2]=1 について求められたら、 X[i][2] に格納される。このようにして、 X[i][N] まで求められたら、それは逆行列となっているので、逆行列を表示する関数 inverse を呼び出す。

関数 inverse では、 for 文を 2 重に回すことで、 X[i][j] を順に表示させる。そして、 1 つの行について内側の for 文でまわして、表示されたら、改行して、次の行の表示を行うというようにしてある。

1.3.4 実行結果

1. 2 行 2 列の逆行列

How many dementions ?
Don't over 30 !
Please input:2
Input Matrix  !
A[11] = 1
A[12] = 2
A[21] = 3
A[22] = 4
********** you inputted this matrix **********
  1.000000  2.000000
  3.000000  4.000000
********** under-triangle matrix 'L' **********
  1.000000   0.000000
  3.000000   1.000000
********** upper-triangle matrix 'U' **********
  1.000000  2.000000
  0.000000 -2.000000
********** inverse matrix ‘ X ’ **********
 -2.000000   1.000000
  1.500000  -0.500000

この実行結果は手計算で求めたものと一致した。

2. 3 行 3 列の逆行列

How many dementions ?
Don't over 30 !

Please input:3
Input Matrix  !
A[11] = -6
A[12] = 9
A[13] = -1
A[21] = 7
A[22] = -9
A[23] = 3
A[31] = -5
A[32] = 4
A[33] = -5
********** you inputted this matrix **********
 -6.000000  9.000000 -1.000000
  7.000000 -9.000000  3.000000
 -5.000000  4.000000 -5.000000
********** under-triangle matrix 'L' **********
  1.000000   0.000000   0.000000
 -1.166667   1.000000   0.000000
  0.833333  -2.333333   1.000000
********** upper-triangle matrix 'U' **********
 -6.000000  9.000000 -1.000000
  0.000000  1.500000  1.833333
  0.000000  0.000000  0.111111
********** inverse matrix ‘ X ’ **********
-33.000000 -41.000000 -18.000000
-20.000000 -25.000000 -11.000000
 17.000000  21.000000   9.000000

この実行結果は手計算で求めたものと一致した。

1.3.5 考察

まず、 LU 分解法を利用して逆行列を求める場合、 A が逆行列を持つとして成り立つ式 AA^{−1} = E を、 x_{i} を A^{−1} の第 i 列として連立方程式が適用できる形に書き直すと、 Ax_{i} = e_{i} となる。これを LU 分解法を使って、繰り返し解いて A^{−1} を求める。例として、 A が 3 行 3 列の行列のときに、 A^{−1} が求められる様子を検証してみる。

A ( \begin{matrix}x_{11} \\ x_{21} \\ x_{31} \end{matrix} ) = ( \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ) A ( \begin{matrix}x_{12} \\ x_{22} \\ x_{32} \end{matrix} ) = ( \begin{matrix}0 \\ 1 \\ 0 \end{matrix} ) A ( \begin{matrix}x_{13} \\ x_{23} \\ x_{33} \end{matrix} ) = ( \begin{matrix}0 \\ 0 \\ 1 \end{matrix} )

これを LU 分解によって解いていくのである。そして、その結果から以下のような A^{−1} を求める。

A^{−1} = ( \begin{matrix}x_{11} x_{12} x_{13} \\ x_{21} x_{22} x_{23} \\ x_{31} x_{32} x_{33} \end{matrix} )

この時、３回連立方程式を解くが、左辺は変わっていない。よって、ガウスの消去法を使うよりも LU 分解を使う方が、左辺の結果をとっておいて計算するため、計算量が少なくてすむ。

1.4 問題 4

1.4.1 問題行列

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

の LU 分解を求めよ。

1.4.2 机上による計算結果クラウト法を用いて実際に 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} ]

となり、解を得る。

1.4.3 プログラムによる実行結果

問題 2 において作成した LU 分解を行うプログラムを改良し、方程式の右辺の入力と方程式の解の計算ならびに表示をする部分を削除したものを作った。つまり、 L と U に分解し、それを表示させるプログラムでこの行列を入力し、 L と U を求める。

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

**** 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

1.4.4 考察

机上で求めた結果とプログラムから求めた結果が一致したので、この結果は正しいものであると考えられる。

1.5 問題 5

1.5.1 問題

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.5.2 机上による計算結果

4 で求めた、 L と U を利用する。

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} ]

1. 計算結果 (1)

[ \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_{1} + y_{2} = 0
3y_{1} + 4y_{2} + y_{3} = 0

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

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

となり、

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

が得られる。

ここで、

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

とすると、

[ \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_{1} + 2x_{2} + 3x_{3} = 1
2x_{1} + 5x_{2} = −2
6x_{3} = 5

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

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

となり、

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

解を得る。

2. 計算結果 (2)

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

となり、

y_{1} = 1
2y_{1} + y2 = 2
3y_{1} + 4y_{2} + y_{3} = 3

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

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

となり、

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

が得られる。ここで、

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

とすると、

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

を得る。故に、

x_{1} + 2x_{2} + 3x_{3} = 1
2x_{1} + 5x_{2} = 0
6x_{3} = 0

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

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

となり、

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

解を得る。

3. 計算結果 (3)

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

となり、

y_{1} = −1
2y_{1} + y_{2} = 0
3y_{1} + 4y_{2} + y_{3} = 1

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

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

となり、

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

が得られる。ここで、

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

とすると、

[ \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_{1} + 2x_{2} + 3x_{3} = −1
2x_{1} + 5x_{2} = 2
6x_{3} = −5

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

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

となり、

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

解を得る。

1.5.3 プログラムによる実行結果

問題 2 で作った LU 分解をし、連立一次方程式の近似解を求めるプログラムに入力し、実行する。

1. 実行結果 (1)

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] = 25
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

2. 実行結果 (2)

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

3. 実行結果 (3)

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

1.5.4 考察

考え方としては以下の通りである。

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

として、

Ax = b

を解く。

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

A = LU

から、

Ax = L(UX) = b

を得る。故に、

UX = y

とおいて

Ly = b

を解き、

y = L^{−1}b

を求め、次に

Ux = y

から

x = U^{−1}y

を求めればよい。

手計算で求めたものは分数で解を出し、プログラムによる実行結果では、浮動小数点数による近似解を求めたが、それぞれ一致しているといえるので、それらの解は正しいといえる。

2 性質 3 と 5 と 6 を証明せよ。 (性質 1,2,4 は仮定してよい)

2.1 性質 3 の証明

2.1.1 性質 3

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

∥RA − I∥ < 1

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

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

が成り立つ。

2.1.2 証明 (性質 3)

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

バナッハの摂動定理 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∥}

が得られる。

2.2 性質 5 の証明

2.2.1 性質 5

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

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

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

2.2.2 証明 (性質 5)

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}|

を満たす。

2.3 性質 6 の証明

2.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}

2.3.2 証明 (性質 6)

|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

3 感想

1. ガウスの消去法今回は、ガウスの消去法という比較的簡単なプログラムであったためにそれほど時間はかからなかった。しかし、出力の際の見栄えなどを気にしてしまい、そのあたりで時間を要した。ガウスの消去法を考える際に、 3 変数くらいの具体的な連立方程式を自分で解いて、その手で解いた式の流れをもとにして、前進消去を考えた。また、後退代入も具体例から規則性を見つけることで解くことができた。

2. LU 分解法まず、手計算により、 3 行 3 列の行列を実行結果としてものせたが、実際に解いて見ることで、どのように L と U の行列が変化しているかということを見ることでプログラムの流れを考えた。 L と U のどちらの対角線要素を 1 としてもよかったが、教科書にならって、 L の方を 1 にすることにした。 LU 分解の計算では for 文が何重にも重なっていたので、その部分を理解することに時間を要したが、解説を読んでもらえればわかると思うが、プログラムの流れはほぼ完璧に理解した。

3. 逆行列ほとんど、 LU 分解のプログラムだったが、解の行列と単位行列をそれぞれ左から 1 列ずつ求めていくということがポイントである。 1 列ずつ求め、その際に関数 solve を呼び出していた。この部分の工夫だけで、 LU 分解を利用して逆行列を求めるプログラムを作ることができた。

4. 証明性質 3 についてはいまいち自信が持てなかったが、教科書を読んで、これ以外に思いつかなかった。また、性質 5 についても、理解があやふやな部分もあるので、今後の学習でより理解を深めていきたい。性質 6 は比較的簡単に証明することができた。

全体としては、本当に頑張ってレポートに取り組んだので、膨大な量となってしまった。また、問題量が多かったために、時間もかなりかかってしまった。そのおかげで、 C 言語、数値計算、 LaTex、それぞれの知識がかなり深まったと思う。また、プログラムもそれぞれ、自分で考え組み立てたので、ほぼ完璧に理解できている。力は確実についたと思うが、かなりきつかった。今後のレポート、授業にも真剣に取り組んで、さらに、数値計算、 C 言語に慣れ親しみ、習得したいと思う。