第二回課題 6/1 出題、大石教授担当分

g99p134-9 宮之前謙一 2001 年 6 月 14 日

1 課題

1.1 教科書「数値計算」 P.31 の問題３の問題をすべて解け

１．ガウスの消去法のプログラムを作り、連立一次方程式の近似解を求めよ。最初はピボット選択なしで、次にピボット選択をするプログラムを作れ。２．ＬＵ分解法のプログラムを作り、連立一次方程式の近似解を求めてみよ。３．ＬＵ分解法を利用して逆行列を求めるプログラムを作れ。４．行列

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

のＬＵ分解を求めよ。５．４の結果を用いて、

(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 解答

2.1 問一・ピボット選択なしのプログラム

まずは、ピボット選択なしのプログラムを下に載せる。プログラムは全部で３つ。ヘッダーファイルを使って main 関数を持った gauss.c というプログラムと、ヘッダーファイル gauss.h と、そこで宣言された関数の中身を載せた、 gaussmath.c である。今回作ったプログラムは、 main 関数を呼び出すと、結果の表示まで全てを行うプログラムになっている。

/****************** gauss.c *******************************/
#include<stdio.h>
#include"gauss.h"
int main(){
  gauss();
  return 0;
}
_____________________________________________________________________
/******************* gauss.h ******************************/
void gauss();

＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿

/****************** gaussmath.c *******************************/
#include<stdio.h>
#include<stdlib.h>
double a[30][30], b[30], result[30];
int ma;
void init_gauss(){
  int i,j;

  for(i=0;i<ma;i++){
     printf("Please type %d equation's %d number!\n",i+1,ma);
     for(j=0;j<ma;j++){
        scanf("%lf", &a[i][j]);
     }
 printf("Please type its vector's 1 number!\n");
 scanf("%lf", &b[i]);
  }
}
void print_gauss(){
  int i, j;
  printf("Your equations are ...\n");
  for(i=0; i<ma; i++) {
    for(j=0; j<ma-1; j++) {
       printf(" %lf*X%d +",a[i][j],j+1);
}
printf(" %lf*X%d = %lf\n",a[i][ma-1],ma,b[i]);
    }
    printf("\n");
}

void print_answer(){
  int i;
  printf("***The answer*** \n");
  for(i=0; i<ma; i++){
      printf(" X%d = %lf \n",i+1,result[i]);
  }
}
void g() {

int i, j, k;
double tmp, keika, pp;
/**** forward elimination ***/
for(i=0; i<ma; i++) {
  tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is not 1 ***/
  if(tmp!=1) {
    for(j=i; j<ma; j++) {
   a[i][j] /= tmp;
    }
    b[i] /= tmp;
  }
  tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is 1 ***/
  for(j=i+1; j<ma;j++) {
     pp = a[j][i] / tmp;
     for(k=i;k<ma;k++) {
    a[j][k] -= pp * a[i][k];
     }
     b[j] -= pp * b[i];
  }
}
/**** backword elimination ****/
for(i=ma-1;i>=0;i--) {
   keika=0.0;
   for(j=i+1;j<ma;j++) {
      keika += a[i][j]*result[j];
   }
   result[i] = (b[i] - keika) / a[i][i];
}
}
void gauss() {
int i, j, k, re;
double time1, time2;
printf("I can answer the equation by using Gauss's elimination!\n");
printf("Please select the equation size under 30.\n");
scanf("%d",&ma);
printf("%dquations!\n",ma);
init_gauss();
printf("The x's demention is %d\n",ma);

for(i=0;i<ma;i++) {
  result[i]=0.0;
}
print_gauss();
for(i=0;i<ma;i++) {
  if(a[0][i]==0) {
    for(j=1;j<ma;j++) {
    if(a[j][i]!=0)  break;
   if(a[ma-1][i]==0) {
     printf("I can't answer.Because a group of a vertical numbers are 0!\n");
     exit(-1);
}
    }
  }
}
for(k=0;k<ma;k++) {
 if(a[k][0]!=0) {
   for(i=k+1;i<ma;i++) {
      re = (int)a[i][0] % (int)a[k][0];
      if(re==0) {
   time1 = a[i][0] / a[k][0];
for(j=0;j<ma;j++) {
if(a[k][j]!=0) {
           re = (int)a[i][j] % (int)a[k][j];
           if(re!=0) break;
time2 = a[i][j] / a[k][j];
       if(time1!=time2) break;
} else {
if(a[i][j]!=a[k][j]) break;
}
       if(j==ma-1) {
          printf("I can't answer.Because you type the same pair of matrixs!\n");
         exit(-1);
}
 }
      }
   }
  } else {
    for(i=k+1;i<ma;i++) {
       if(a[i][0]==0) {
         for(j=0;j<ma;j++) {

              if(a[i][j]!=a[k][j]) break;
     if(j==ma-1) {
            printf("I can't answer.Because you type the same pair of matrixs!\n");
       exit(-1);
     }
           }
         }
      }
    }
  }
  g();
  print_answer();
}

2.2 問一のピボット選択無しのプログラムの解説

まず、 gauss.c についての説明です。これには、 main 関数がおいてあって、ガウスの消去法を行う場合には、 1 番を選択してすれば、 gauss という関数が呼び出される。しかしこれは、このプログラムには無く、ヘッダーファイルの宣言より、リンクされている gaussmath.c の中身を探す。すると、そこにあるので、それを読み出す。 gauss.h では、必要な関数を宣言している。一番に必要な gauss という関数と、二番に必要な current という関数である。 gaussmath.c では、具体的な、消去法などの処理が書かれている。まず始めに、方程式をいくつ連立させるのか聞かれる。それに答えると、自動的に、 init-gauss という関数が呼び出される。そこでは、値の入力が行われる。それが終了すると、次に、結果を表す配列の初期化が行われて、続いて入力した方程式を表示する関数 print-gauss を呼び出して、方程式を表示する。次に、作った方程式が、解けない形でないか、いくつかのパターンについて確認する。このうちうどれかで引っかかった場合は、解が求められなくなるので、プログラムを強制終了する。具体的には、一つ目のパターンとして、同じ列の値が全て０である組み合わせが存在するかを見る。もしそれが存在するならば、解は不定となり、特定できないので、解けない、というメッセージを表示して終わる。調べ方は、まず 1 行目の値を調べていき、 0 であるものが見つかったら、その列全ての数を調べて、どれかひとつでも 0 でなかったら、探索を終了する。というものだ。二つ目のパターンとして、見た目は違う形の方程式でも、ある一つの方程式の全ての数を定数倍した数である場合もある。この場合は、その方程式は、ひとつの方程式とは考えられない無意味なものとなる。よって、解が不定となり、解くことができない。具体的な発見方法は、まず、各行の同じ列について、ある列の値が、調べている元の行列の値の倍数であるかどうかを調べて、もしも同じ値であったら、実際に何倍かを調べて、他の列についても同様に検査して、全ての列が、同じ倍数でないか調べる。もしそうであったら、貝が求められない旨を表示して終了するようにする。この処理をするときに気をつけないといけないことは、もしも、注目している列が 0 であった時は、もう片方の項も 0 であるという点だった。とりあえず今の段階では解ける形であると分かったら、処理を行う関数 g を呼び出す。その後、正常にガウスの消去法が終了すれば、答えを出力する関数、 print-answer を呼び出して処理を終わる。関数 g の説明に移ると、そこでは、まず、前進消去が行われる。まず、 i 行 i 列の値が 1 であるか見る。違った場合は、値を 1 にするように、方程式を操作する。 1 であったら、他の各行の、同じ列の値を 0 にするように計算を行う。この作業を、全ての行について行う。この作業が全て終了すると、続いて、後退代入に移る。前進消去を行った下の行から答えを求めていって、全ての答えを求める。

2.3 実行結果

2.3.1 実行結果１

まずは、探索が成功する例を示す。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
4
4quations!
Please type 1 equation's 4 number!
2 2 4 2
Please type its vector's 1 number!
26
Please type 2 equation's 4 number!
2 3 4 3
Please type its vector's 1 number!
32
Please type 3 equation's 4 number!
1 2 1 2
Please type its vector's 1 number!
16
Please type 4 equation's 4 number!
1 2 2 1
Please type its vector's 1 number!
15
The x's demention is 4
Your equations are ...
 2.000000*X1 + 2.000000*X2 + 4.000000*X3 + 2.000000*X4 = 26.000000
 2.000000*X1 + 3.000000*X2 + 4.000000*X3 + 3.000000*X4 = 32.000000
 1.000000*X1 + 2.000000*X2 + 1.000000*X3 + 2.000000*X4 = 16.000000
 1.000000*X1 + 2.000000*X2 + 2.000000*X3 + 1.000000*X4 = 15.000000
***The answer***
 X1 = 1.000000
 X2 = 2.000000
 X3 = 3.000000
 X4 = 4.000000
Press any key to continue

これにより、整数入力だと探索が成功することが分かる。 2.3.2 実行結果２

探索が成功する例。

I can answer the equation by using Gauss's elimination
Please select the equation size under 30.
3
3quations!
Please type 1 equation's 3 number!
2 1 3
Please type its vector's 1 number!
16.05
Please type 2 equation's 3 number!
5.2 2.2 1.7
Please type its vector's 1 number!
18.14
Please type 3 equation's 3 number!
3.9 4.2 0.7
Please type its vector's 1 number!
17.98
The x's demention is 3
Your equations are ...
 2.000000*X1 + 1.000000*X2 + 3.000000*X3 = 16.050000
 5.200000*X1 + 2.200000*X2 + 1.700000*X3 = 18.140000
 3.900000*X1 + 4.200000*X2 + 0.700000*X3 = 17.980000
***The answer***
 X1 = 1.200000
 X2 = 2.550000
 X3 = 3.700000
Press any key to continue

少数を入れた場合も成功している。

2.3.3 実行結果３

もう一つ、探索が成功する例。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
3
3quations!
Please type 1 equation's 3 number!
2 1 3
Please type its vector's 1 number!
31

Please type 2 equation's 3 number!
1 1 1
Please type its vector's 1 number!
99
Please type 3 equation's 3 number!
5 6 98
Please type its vector's 1 number!
21
The x's demention is 3
Your equations are ...
 2.000000*X1 + 1.000000*X2 + 3.000000*X3 = 31.000000
 1.000000*X1 + 1.000000*X2 + 1.000000*X3 = 99.000000
 5.000000*X1 + 6.000000*X2 + 98.000000*X3 = 21.000000
***The answer***
 X1 = -54.361702
 X2 = 160.180851
 X3 = -6.819149
Press any key to continue

この場合もまた、探索が成功する。

2.3.4 実行結果４

次に、 0 を含む形だが、ピボット選択を必要としない例を入力してみる。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
3
3quations!
Please type 1 equation's 3 number!
2 0 1
Please type its vector's 1 number!
12
Please type 2 equation's 3 number!
1 2 3
Please type its vector's 1 number!
24
Please type 3 equation's 3 number!
2 2 1
Please type its vector's 1 number!
15
The x's demention is 3
Your equations are ...
 2.000000*X1 + 0.000000*X2 + 1.000000*X3 = 12.000000

 1.000000*X1 + 2.000000*X2 + 3.000000*X3 = 24.000000
 2.000000*X1 + 2.000000*X2 + 1.000000*X3 = 15.000000
***The answer***
 X1 = 3.000000
 X2 = 1.500000
 X3 = 6.000000
Press any key to continue

これも、しっかりと処理を行ってくれる。

2.3.5 実行結果５

ピボット選択についての話をする前に、入力から、解けないパターンの方程式がいくつかあるので、それについて、下に実行例を挙げてみる。まずは、ある行列の全ての数の倍数の形の方程式を含むものを入れてみる。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
4
4quations!
Please type 1 equation's 4 number!
1 2 3 4
Please type its vector's 1 number!
30
Please type 2 equation's 4 number!
2 4 6 8
Please type its vector's 1 number!
60
Please type 3 equation's 4 number!
2 1 2 3
Please type its vector's 1 number!
20
Please type 4 equation's 4 number!
4 5 10 2
Please type its vector's 1 number!
121
The x's demention is 4
Your equations are ...
 1.000000*X1 + 2.000000*X2 + 3.000000*X3 + 4.000000*X4 = 30.000000
 2.000000*X1 + 4.000000*X2 + 6.000000*X3 + 8.000000*X4 = 60.000000
 2.000000*X1 + 1.000000*X2 + 2.000000*X3 + 3.000000*X4 = 20.000000
 4.000000*X1 + 5.000000*X2 + 10.000000*X3 + 2.000000*X4 = 121.000000
I can't answer.Because you type the same pair of matrixs!

Press any key to continue

これも、失敗する。

2.3.6 実行結果６

今度は、ある列の全ての数が 0 である方程式を入れてみる。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
4
4quations!
Please type 1 equation's 4 number!
2 0 2 2
Please type its vector's 1 number!
10
Please type 2 equation's 4 number!
3 0 5 1
Please type its vector's 1 number!
17
Please type 3 equation's 4 number!
9 0 3 4
Please type its vector's 1 number!
23
Please type 4 equation's 4 number!
7 0 0 8
Please type its vector's 1 number!
20
The x's demention is 4
Your equations are ...
 2.000000*X1 + 0.000000*X2 + 2.000000*X3 + 2.000000*X4 = 10.000000
 3.000000*X1 + 0.000000*X2 + 5.000000*X3 + 1.000000*X4 = 17.000000
 9.000000*X1 + 0.000000*X2 + 3.000000*X3 + 4.000000*X4 = 23.000000
 7.000000*X1 + 0.000000*X2 + 0.000000*X3 + 8.000000*X4 = 20.000000
I can't answer.Because a group of a vertical numbers are 0!
Press any key to continue

これも、失敗する。

2.3.7 実行結果７

次に、ピボット選択を行っていないので、解けない方程式を入れてみる。ピボット選択は、注目している列の数が、 0 である場合を入れてやればよい。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
3
3quations!
Please type 1 equation's 3 number!
0 2 1
Please type its vector's 1 number!
12
Please type 2 equation's 3 number!
2 0 3
Please type its vector's 1 number!
-20
Please type 3 equation's 3 number!
3 2 0
Please type its vector's 1 number!
10
The x's demention is 3
Your equations are ...
 0.000000*X1 + 2.000000*X2 + 1.000000*X3 = 12.000000
 2.000000*X1 + 0.000000*X2 + 3.000000*X3 = -20.000000
 3.000000*X1 + 2.000000*X2 + 0.000000*X3 = 10.000000
***The answer***
 X1 = -1.#IND00
 X2 = -1.#IND00
 X3 = -1.#IND00
Press any key to continue

このように、解が正しく出てこない。これは、一番最初の列が１ではない処理行っているが、 0 である場合は特別で、プログラムで書かれている方法では対処しきれないためである。次で、これを解消する。

2.4 部分ピボット選択を入れたプログラム

部分ピボット選択を入れた gaussmath.c の関数 g のプログラムを下に示す。

void g() {
  int i, j, k, maxabv;
  double tmp, abvalue, change, keika, pp;
/**** forward elimination ***/
  for(i=0; i<ma; i++) {
  tmp = a[i][i];
 /***** the case of the number of matrix[i][i] is 0 ***/
    if(tmp==0) {
      maxabv = 0;        /* first definition */

for(j=i+1; j<ma; j++) {
if(a[j][i]<0) {
abvalue = -a[j][i];
} else {
     abvalue = a[j][i];
}
if(abvalue>tmp) {
tmp = a[j][i];
     maxabv = j;
}
}
if(maxabv!=0) {
  for(k=0; k<ma; k++) {
     change = a[i][k];
     a[i][k] = a[maxabv][k];
     a[maxabv][k] = change;
  }
  change = b[i];
b[i] = b[maxabv];
b[maxabv] = change;
} else {
printf("I can't answer these equations!\n");
exit(-1);
}
}
tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is not 1 ***/
if(tmp!=1) {
for(j=i; j<ma; j++) {
a[i][j] /= tmp;
}
b[i] /= tmp;
}
tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is 1 ***/
for(j=i+1; j<ma;j++) {
 pp = a[j][i] / tmp;
 for(k=i;k<ma;k++) {
a[j][k] -= pp * a[i][k];
 }
 b[j] -= pp * b[i];
}
}

/**** backword elimination ****/
  for(i=ma-1;i>=0;i--) {
     keika=0.0;
     for(j=i+1;j<ma;j++) {
        keika += a[i][j]*result[j];
     }
     result[i] = (b[i] - keika) / a[i][i];
  }
}

2.5 プログラムの解説

ピボット選択は、関数 g の一番最初に置いてある。まず、入力された配列 a の、 i 行 i 列の値が 0 であるか見ている。 0 であった場合は、部分ピボット選択を行って、一番絶対値が大きいものの行全ての数と、今見ている行の数を、全て交換する。もしもここで、絶対値が大きいものがない（つまり、全ての列が 0 である）という結果になった場合は、解けない旨を表示してプログラムを強制終了する。この場合は、一体どういう場合なのかというと、結果を参照して欲しい。今までとはまた違った形の、解を求められない方程式を弾くためのものだ。

2.6 実行結果

2.6.1 実行結果８

先ほど失敗した探索が失敗している、ピボット選択を行う必要がある例を入力してみる。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
3
3quations!
Please type 1 equation's 3 number!
0 2 1
Please type its vector's 1 number!
12
Please type 2 equation's 3 number!
2 0 3
Please type its vector's 1 number!
-20
Please type 3 equation's 3 number!
3 2 0
Please type its vector's 1 number!
10
The x's demention is 3
Your equations are ...
 0.000000*X1 + 2.000000*X2 + 1.000000*X3 = 12.000000
 2.000000*X1 + 0.000000*X2 + 3.000000*X3 = -20.000000

 3.000000*X1 + 2.000000*X2 + 0.000000*X3 = 10.000000
***The answer***
 X1 = -2.363636
 X2 = 8.545455
 X3 = -5.090909
Press any key to continue

このように、探索が成功している。

2.6.2 実行結果９

最後に、パターンがつかみづらいが、失敗するものを入れてみる。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
5
5quations!
Please type 1 equation's 5 number!
1 2 1 1 1
Please type its vector's 1 number!
6
Please type 2 equation's 5 number!
2 4 2 1 2
Please type its vector's 1 number!
11
Please type 3 equation's 5 number!
3 6 3 1 3
Please type its vector's 1 number!
16
Please type 4 equation's 5 number!
1 1 1 1 1
Please type its vector's 1 number!
5
Please type 5 equation's 5 number!
5 6 2 1 3
Please type its vector's 1 number!
17
The x's demention is 5
Your equations are ...
 1.000000*X1 + 2.000000*X2 + 1.000000*X3 + 1.000000*X4 + 1.000000*X5 = 6.000000
 2.000000*X1 + 4.000000*X2 + 2.000000*X3 + 1.000000*X4 + 2.000000*X5 = 11.000000
 3.000000*X1 + 6.000000*X2 + 3.000000*X3 + 1.000000*X4 + 3.000000*X5 = 16.000000

 1.000000*X1 + 1.000000*X2 + 1.000000*X3 + 1.000000*X4 + 1.000000*X5 = 5.000000
 5.000000*X1 + 6.000000*X2 + 2.000000*X3 + 1.000000*X4 + 3.000000*X5 = 17.000000

I can't answer these equations!
Press any key to continue

このように、失敗する。この失敗が、先ほど解説の所で、見落としせずに、失敗すると忘れずに表示をしなくてはいけない、と言った所だ。この場合は、係数の値が似かより過ぎているため、必要以上に値が消えてしまう。手計算でも、行き詰まってしまう。答えは、全ての X が 1 になるようにしていて、式の計算も合っているが、それに関係なく失敗する。実際の処理では、順調に最後まで前進消去で消していって、最後の 4 行目のループの所で、 4 行 4 列目の値が 0 であるので、下の行の値を探そうとするが、そんな行は無いため、計算ができない、とエラーを出して終わる。同じ理由で、次のような例も失敗する。

2.6.3 実行結果１０

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
5
5quations!
Please type 1 equation's 5 number!
1 2 1 1 1
Please type its vector's 1 number!
6
Please type 2 equation's 5 number!
2 4 2 1 2
Please type its vector's 1 number!
11
Please type 3 equation's 5 number!
3 6 3 1 3
Please type its vector's 1 number!
16
Please type 4 equation's 5 number!
4 8 4 1 4
Please type its vector's 1 number!
21
Please type 5 equation's 5 number!
5 10 5 1 5
Please type its vector's 1 number!
26
The x's demention is 5
Your equations are ...
 1.000000*X1 + 2.000000*X2 + 1.000000*X3 + 1.000000*X4 + 1.000000*X5 = 6.000000
 2.000000*X1 + 4.000000*X2 + 2.000000*X3 + 1.000000*X4 + 2.000000*X5 = 11.000000

 3.000000*X1 + 6.000000*X2 + 3.000000*X3 + 1.000000*X4 + 3.000000*X5 = 16.000000
 4.000000*X1 + 8.000000*X2 + 4.000000*X3 + 1.000000*X4 + 4.000000*X5 = 21.000000
 5.000000*X1 + 10.000000*X2 + 5.000000*X3 + 1.000000*X4 + 5.000000*X5 = 26.000000
I can't answer these equations!
Press any key to continue

今までの実行結果のように、結構、失敗する形はある。適当に方程式を入れても、それが失敗に終わることは、考えていたより結構あり得るのだ。

2.7 問二・ＬＵ分解法のプログラム

先ほど一番で説明したプログラムを少し変える。まずはプログラムを載せる。上のプログラムの gauss.c を lu.c、 gauss.h を lu.c、 gaussmath.c を lumath.c に変える。ここでは、 lumath.c しか変えていないので、それだけ載せる。

/************* lumath.c ***************/
#include<stdio.h>
#include<stdlib.h>
double a[30][30], b[30], result[30];
int change[30];
double notone[30], divide[30][30];
int ma;
void init_lu(){
  int i,j;
  for(i=0;i<ma;i++){
     printf("Please type %d equation's %d number!\n",i+1,ma);
     for(j=0;j<ma;j++){
        scanf("%lf", &a[i][j]);
     }
 printf("Please type its vector's 1 number!\n");
 scanf("%lf", &b[i]);
  }
}
void print_lu(){
  int i, j;
  printf("Your equations are ...\n");
  for(i=0; i<ma; i++) {

  for(j=0; j<ma-1; j++) {
     printf(" %lf*X%d +",a[i][j],j+1);
}
printf(" %lf*X%d = %lf\n",a[i][ma-1],ma,b[i]);
  }
  printf("\n");
}

void print_answer(){
int i;
printf("***The answer*** \n");
for(i=0; i<ma; i++){
    printf(" X%d = %lf \n",i+1,result[i]);
}
}
void lu_line(double a[30][30]) {
int i, j, k;
double tmp, abvalue, ch, pp;
/**** forward elimination ***/
for(i=0; i<ma; i++) {
tmp = a[i][i];
/***** the case of the number of matrix[i][i] is 0 ***/
  if(tmp==0) {
    change[i] = 0;        /* first definition */
    for(j=i+1; j<ma; j++) {
if(a[j][i]<0) {
 abvalue = -a[j][i];
} else {
                 abvalue = a[j][i];
}
   if(abvalue>tmp) {
 tmp = a[j][i];
                 change[i] = j;
   }
    }
    if(change[i]!=0) {
      for(k=0; k<ma; k++) {
         ch = a[i][k];
         a[i][k] = a[change[i]][k];
         a[change[i]][k] = ch;
      }

} else {
printf("I can't answer these equations!\n");
exit(-1);
}
}
tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is not 1 ***/
if(tmp!=1) {
notone[i] = tmp;
for(j=i; j<ma; j++) {
a[i][j] /= tmp;
}
}
tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is 1 ***/
for(j=i+1; j<ma;j++) {
pp = a[j][i] / tmp;
divide[j][i] = pp;
for(k=i;k<ma;k++) {
a[j][k] -= pp * a[i][k];
}
}
}
}
void lu_vector(double a[30][30], double b[30]) {
int i, j;
double ch, keika;
/**** forward elimination ***/
for(i=0; i<ma; i++) {
/***** the case of the number of matrix[i][i] is 0 ***/
if(change[i]!=0) {
ch = b[i];
b[i] = b[change[i]];
b[change[i]] = ch;
}
/*** the case of the  number of matrix[i][i] is not 1 ***/
if(notone[i]!=0) {
b[i] /= notone[i];
}

/*** the case of the  number of matrix[i][i] is 1 ***/
   for(j=i+1; j<ma;j++) {
      b[j] -= divide[j][i] * b[i];
   }
}
print_lu();
/**** backword elimination ****/
for(i=ma-1;i>=0;i--) {
   keika=0.0;
   for(j=i+1;j<ma;j++) {
      keika += a[i][j]*result[j];
   }
   result[i] = (b[i] - keika) / a[i][i];
}
}

void lu() {
int i, j, k, re;
double time1, time2;
printf("I can answer the equation by using Gauss's elimination!\n");
printf("Please select the equation size under 30.\n");
scanf("%d",&ma);
printf("%dquations!\n",ma);
init_lu();
printf("The x's dimention is %d\n",ma);
for(i=0;i<ma;i++) {
  result[i]=0.0;
}
print_lu();
for(i=0;i<ma;i++) {
  if(a[0][i]==0) {
    for(j=1;j<ma;j++) {
    if(a[j][i]!=0)  break;
   if(a[ma-1][i]==0) {
     printf("I can't answer.Because a group of a vertical numbers are 0!\n");
     exit(-1);
}
    }

  }
}
for(k=0;k<ma;k++) {
 if(a[k][0]!=0) {
   for(i=k+1;i<ma;i++) {
      re = (int)a[i][0] % (int)a[k][0];
      if(re==0) {
   time1 = a[i][0] / a[k][0];
for(j=0;j<ma;j++) {
if(a[k][j]!=0) {
           re = (int)a[i][j] % (int)a[k][j];
           if(re!=0) break;
time2 = a[i][j] / a[k][j];
       if(time1!=time2) break;
} else {
if(a[i][j]!=a[k][j]) break;
}
       if(j==ma-1) {
          printf("I can't answer.Because you type the same pair of matrixs!\n");
         exit(-1);
}
 }
      }
   }
  } else {
    for(i=k+1;i<ma;i++) {
       if(a[i][0]==0) {
         for(j=0;j<ma;j++) {
            if(a[i][j]!=a[k][j]) break;
   if(j==ma-1) {
          printf("I can't answer.Because you type the same pair of matrixs!\n");
     exit(-1);
   }
         }
       }
    }
  }
}
for(i=0;i<ma;i++) {
   change[i] = i;
   notone[i] = 0;
}
lu_line(a);
lu_vector(a,b);

  print_answer();
  for(i=0;i<ma;i++) {
  b[i] = i*i;
  }
  lu_vector(a,b);
  print_answer();
}

2.8 プログラムの解説

行列についての処理と、ベクトルの処理を分けるために、新しく配列を定義してある。ピボット選択を行った際に、入れ替えた行列を記憶しておく int 型の配列 change と、ある行の１列目が１でなかった時に、その行を、１列目が１となるように調整に使った値を残しておく配列 notone と、前進消去に使った値を残しておく配列 divide である。これ以外にも、行列の処理とベクトルの処理を分けたときに、この後、逆行列の問題を解く為に必要な、引数を設けた。

2.9 実行結果

下に、実行結果を示す。ここでは、プログラムの中に b の値を二つ入力させている。プログラムでは、Ａの値を同じにして、計算することが可能になったことを示す。

I can answer the equation by using Gauss's elimination!
Please select the equation size under 30.
4
4quations!
Please type 1 equation's 4 number!
2 2 4 2
Please type its vector's 1 number!
26
Please type 2 equation's 4 number!
2 3 4 3
Please type its vector's 1 number!
32
Please type 3 equation's 4 number!
1 2 1 2
Please type its vector's 1 number!
16
Please type 4 equation's 4 number!
1 2 2 1
Please type its vector's 1 number!
15
The x's dimention is 4
Your equations are ...

2.000000*X1 + 2.000000*X2 + 4.000000*X3 + 2.000000*X4 = 26.000000
2.000000*X1 + 3.000000*X2 + 4.000000*X3 + 3.000000*X4 = 32.000000
1.000000*X1 + 2.000000*X2 + 1.000000*X3 + 2.000000*X4 = 16.000000
1.000000*X1 + 2.000000*X2 + 2.000000*X3 + 1.000000*X4 = 15.000000
Your equations are ...
1.000000*X1 + 1.000000*X2 + 2.000000*X3 + 1.000000*X4 = 13.000000
0.000000*X1 + 1.000000*X2 + 0.000000*X3 + 1.000000*X4 = 6.000000
0.000000*X1 + 0.000000*X2 + 1.000000*X3 + 0.000000*X4 = 3.000000
0.000000*X1 + 0.000000*X2 + 0.000000*X3 + 1.000000*X4 = 4.000000
***The answer***
X1 = 1.000000
X2 = 2.000000
X3 = 3.000000
X4 = 4.000000
Your equations are ...
1.000000*X1 + 1.000000*X2 + 2.000000*X3 + 1.000000*X4 = 0.000000
0.000000*X1 + 1.000000*X2 + 0.000000*X3 + 1.000000*X4 = 1.000000
0.000000*X1 + 0.000000*X2 + 1.000000*X3 + 0.000000*X4 = -3.000000
0.000000*X1 + 0.000000*X2 + 0.000000*X3 + 1.000000*X4 = -8.000000
***The answer***
X1 = 5.000000
X2 = 9.000000
X3 = -3.000000
X4 = -8.000000
Press any key to continue

二つ目の式の表示は、行列の値はすでに操作されているので、入力した値とは違う。だが、これで、行列の処理とベクトルの処理を分かれやっていることが分かるだろう。

2.10 問三・逆行列のプログラム

下に、逆行列のプログラムを載せる。

/**** lugyaku.c ****/
#include<stdio.h>
#include<stdlib.h>
double a[30][30], a_back[30][30], result[30];
double notone[30], divide[30][30];
int change[30];
int ma;

void init_lu(){
int i,j;

for(i=0;i<ma;i++){
   printf("Please type %d line's %d number!\n",i+1,ma);
   for(j=0;j<ma;j++){
      scanf("%lf", &a[i][j]);
   }
}
}
void print_lu(double c[30][30]){
int i, j;
printf("Your equations are ...\n");
for(i=0; i<ma; i++) {
  for(j=0; j<ma; j++) {
     printf(" %lf ",c[i][j]);
  }
  printf("\n");
}
}
void print_answer(double c[30][30]){
int i, j;
printf("***The answer*** \n");
for(i=0; i<ma; i++) {
  for(j=0; j<ma; j++) {
     printf(" %lf ",c[i][j]);
  }
  printf("\n");
}
}
void lu_line(double a[30][30]) {
int i, j, k;
double tmp, abvalue, ch, pp;
/**** forward elimination ***/
for(i=0; i<ma; i++) {
tmp = a[i][i];
/***** the case of the number of matrix[i][i] is 0 ***/
  if(tmp==0) {
    change[i] = 0;        /* first definition */

for(j=i+1; j<ma; j++) {
if(a[j][i]<0) {
abvalue = -a[j][i];
} else {
             abvalue = a[j][i];
}
if(abvalue>tmp) {
tmp = a[j][i];
             change[i] = j;
}
}
if(change[i]!=0) {
  for(k=0; k<ma; k++) {
     ch = a[i][k];
     a[i][k] = a[change[i]][k];
     a[change[i]][k] = ch;
  }
} else {
printf("I can't answer these equations!\n");
exit(-1);
}
}
tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is not 1 ***/
if(tmp!=1) {
notone[i] = tmp;
for(j=i; j<ma; j++) {
a[i][j] /= tmp;
}
}
tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is 1 ***/
for(j=i+1; j<ma;j++) {
 pp = a[j][i] / tmp;
 divide[j][i] = pp;
 for(k=i;k<ma;k++) {
a[j][k] -= pp * a[i][k];
 }
}
}
}
void lu_vector(double a[30][30], double b[30]) {
int i, j;

double ch, keika;
/**** forward elimination ***/
for(i=0; i<ma; i++) {
/***** the case of the number of matrix[i][i] is 0 ***/
if(change[i]!=0) {
  ch = b[i];
  b[i] = b[change[i]];
  b[change[i]] = ch;
}
/*** the case of the  number of matrix[i][i] is not 1 ***/
if(notone[i]!=0) {
  b[i] /= notone[i];
}
/*** the case of the  number of matrix[i][i] is 1 ***/
for(j=i+1; j<ma;j++) {
   b[j] -= divide[j][i] * b[i];
}
}
/**** backword elimination ****/
for(i=ma-1;i>=0;i--) {
keika=0.0;
for(j=i+1;j<ma;j++) {
   keika += a[i][j]*result[j];
}
result[i] = (b[i] - keika) / a[i][i];
}
}
void make_back(double c[30][30], double c_back[30][30]) {
int i, j;
double back[30];
lu_line(c);
for(k=0;k<ma;k++) {
for(i=0;i<ma;i++) {
   back[i] = 0.0;
}
back[k] = 1.0;
lu_vector(c,back);

   for(i=0;i<ma;i++) {
  c_back[i][k] = result[i];
}
}
}
void lu() {
int i, j;
for(i=0;i<ma;i++) {
for(j=0;j<ma;j++) {
a[i][j] = 0;
a_back[i][j] = 0.0;
}
}
printf("Now begin LUgyaku!\n");
printf("Please select normal line size under 30.\n");
scanf("%d",&ma);
printf("%dlines!\n",ma);
for(i=0;i<ma;i++) {
result[i] = 0.0;
}
init_lu();
print_lu(a);
make_back(a,a_back);
print_answer(a_back);
}
int main() {
 lu();
 return 0;
}

2.11 プログラムの解説

今までのプログラムから残っている部分も、ほとんどが内容を少し変えている。おおきく変わったところは、それぞれの関数に、引数を用意し始めたことだ。これは、次の問５で、今まで作成したすべてのプログラムをくっつけるための準備である。ここでの詳細な関数の説明は、関数 make-back についてのみとする。この関数では、入力された行列の逆行列を作る。内容は簡単で、単位行列のある列を配列 back に入れて、それについてＬＵ分解法を行う、というものだ。これをループで回して、逆行列を得る。

2.12 実行結果

以下に実行結果を示す。

Now begin LUgyaku!
Please select normal line size under 30.
3
3lines!
Please type 1 line's 3 number!
1 3 3
Please type 2 line's 3 number!
1 1 4
Please type 3 line's 3 number!
-1 2 -6
Your equations are ...
 1.000000  3.000000  3.000000
 1.000000  1.000000  4.000000
 -1.000000  2.000000  -6.000000
***The answer***
 -14.000000  24.000000  9.000000
 2.000000  -3.000000  -1.000000
 3.000000  -5.000000  -2.000000
Press any key to continue

このように、逆行列の処理がしっかりと行われていることがわかる。手計算でも、今回の結果と同じ答えが出た。

2.13 問四・ＬＵ分解のプログラム

行列Ａはここでは載せない。下に、プログラムを載せる。

/**** lu4.c ****/
#include<stdio.h>
#include<stdlib.h>
double a[30][30], result[30];
double l[30][30], u[30][30];
int ma;
void init_lu(){
  int i,j;

  for(i=0;i<ma;i++){
     printf("Please type %d line's %d number!\n",i+1,ma);
     for(j=0;j<ma;j++){

      scanf("%lf", &a[i][j]);
   }
}
}
void print_lu(){
int i, j;
printf("Your lines are ...\n");
for(i=0;i<ma;i++) {
  for(j=0;j<ma;j++) {
     printf(" %lf ",a[i][j]);
  }
  printf("\n");
}
printf("L------  \n");
for(i=0;i<ma;i++) {
  for(j=0;j<ma;j++) {
     printf(" %lf ",l[i][j]);
  }
  printf("\n");
}
printf("U------  \n");
for(i=0;i<ma;i++) {
  for(j=0;j<ma;j++) {
     printf(" %lf ",u[i][j]);
  }
  printf("\n");
}
}
void lu_apart() {
int i, j, k;
double total;
/***** first definition *****/
for(i=0;i<ma;i++) {
   l[i][i] = 1.0;
   for(j=i+1;j<ma;j++) {
      l[i][j] = 0.0;
   }
   for(j=0;j<i;j++) {
      u[i][j] = 0.0;
   }

}

for(k=0;k<ma;k++) {      /***** line number *****/
   for(i=k;i<ma;i++) {      /***** retu number *****/
      total = 0.0;
      for(j=0;j<k;j++) {      /***** teisuukou no genzan *****/
         total += l[k][j] * u[j][i];
      }
      u[k][i] = a[k][i] - total;
   }
   for(i=0;i<k+1;i++) {      /***** retu number *****/
      total = 0.0;
      for(j=0;j<i;j++) {
         total += l[k+1][j] * u[j][i];
      }
      l[k+1][i] = (a[k+1][i] - total) / u[i][i];
   }
}
}
void lu() {
printf("Now begin LUbunkai!\n");
printf("Please select normal line size under 30.\n");
scanf("%d",&ma);
printf("%dlines!\n",ma);
init_lu();
lu_apart();
print_lu();
}
int main() {
 lu();
 return 0;
}

2.14 プログラムの解説

このプログラムでは、問３のプログラムほど複雑ではないので、各関数から引数を取っ払って、簡単なプログラムにするよう努めた。ＬＵ分解を行う関数は、 lu-apart であるので、これについて説明する。まず始めに、配列 L の 1 行目は全て値が決まっているので、それを代入する。同時に、 U の 0 の部分も入れてしまう。これを初期値として、ループに入る。ここでは、 U の行を上から見て、値が埋まっていない最初の一行の値を決定して、 L の、同様な状況の一行の値を決定する、という事を繰り返す。行列 U の値の求め方は、 L,U の行列について、 0 が存在することを利用して、行列 L の、値が１である列を、行列 U の、求めたい行と同じになるように行列 L を選べば良い。同様な考えで、行列 L の値については、今、ある行を見ているとして、その行のある列の値を求めたい場合、行列 U から、その列と対応する行番号の一つ下以降の値が 0 であるようなものを選んで計算していけば良い。このようなプログラムを使って、下の実行結果を得る。

2.15 実行結果

以下に、教科書の問題の実行結果を示す。

Now begin LUbunkai!
Please select normal line size under 30.
3
3lines!
Please type 1 line's 3 number!
1 2 3
Please type 2 line's 3 number!
2 8 11
Please type 3 line's 3 number!
3 22 35
Your lines are ...
 1.000000  2.000000  3.000000
 2.000000  8.000000  11.000000
 3.000000  22.000000  35.000000
L------
 1.000000  0.000000  0.000000
 2.000000  1.000000  0.000000
 3.000000  4.000000  1.000000
U------
 1.000000  2.000000  3.000000
 0.000000  4.000000  5.000000
 0.000000  0.000000  6.000000
Press any key to continue

これより、処理ができていることがわかる。

2.16 問五・プログラム

今までのを利用して作ったプログラムを下に載せる。

/**** lu5.c ****/
#include<stdio.h>
#include<stdlib.h>

double a[30][30], l[30][30], u[30][30], b[30], result[30];
double l_back[30][30], u_back[30][30], m[30], n[30];
int change[30];
double notone[30], divide[30][30];
int ma;
void init_lu(){
int i,j;
for(i=0;i<ma;i++){
   printf("Please type %d equation's %d number!\n",i+1,ma);
   for(j=0;j<ma;j++){
      scanf("%lf", &a[i][j]);
   }
printf("Please type its vector's 1 number!\n");
scanf("%lf", &b[i]);
}
}

void print_alu(double c[30][30], double d[30][30], double e[30][30]){
int i, j;
printf("Your lines are ...\n");
for(i=0;i<ma;i++) {
  for(j=0;j<ma;j++) {
     printf(" %lf ",c[i][j]);
  }
  printf("\n");
}
printf("L------  \n");
for(i=0;i<ma;i++) {
  for(j=0;j<ma;j++) {
     printf(" %lf ",d[i][j]);
  }
  printf("\n");
}
printf("U------  \n");
for(i=0;i<ma;i++) {
  for(j=0;j<ma;j++) {
     printf(" %lf ",e[i][j]);
  }
  printf("\n");
}
}

void print_back(double c[30][30], double d[30][30]){
int i, j;
printf("normal ----\n");
for(i=0;i<ma;i++) {
  for(j=0;j<ma;j++) {
     printf(" %lf ",c[i][j]);
  }
  printf("\n");
}
printf("back  ----\n");
for(i=0;i<ma;i++) {
  for(j=0;j<ma;j++) {
     printf(" %lf ",d[i][j]);
  }
  printf("\n");
}
}
void print_lu(double c[30][30],char X, double d[30]){
int i, j;
printf("Equations are ...\n");
for(i=0; i<ma; i++) {
  for(j=0; j<ma-1; j++) {
     printf(" %lf*%c%d +",a[i][j],X,j+1);
}
printf(" %lf*%c%d = %lf\n",c[i][ma-1],X,ma,d[i]);
  }
  printf("\n");
}
void print_ans(char X, double m[30]) {
int i;
for(i=0; i<ma; i++){
    printf(" %c%d = %lf \n",X,i+1,m[i]);
}
}
void print_answer(){
int i;
printf("***The answer*** \n");

for(i=0; i<ma; i++){
printf(" X%d = %lf \n",i+1,result[i]);
}
}
void lu_line(double a[30][30]) {
int i, j, k;
double tmp, abvalue, ch, pp;
/**** forward elimination ***/
for(i=0; i<ma; i++) {
tmp = a[i][i];
/***** the case of the number of matrix[i][i] is 0 ***/
if(tmp==0) {
change[i] = 0;        /* first definition */
for(j=i+1; j<ma; j++) {
if(a[j][i]<0) {
abvalue = -a[j][i];
} else {
             abvalue = a[j][i];
}
if(abvalue>tmp) {
tmp = a[j][i];
             change[i] = j;
}
}
if(change[i]!=0) {
  for(k=0; k<ma; k++) {
     ch = a[i][k];
     a[i][k] = a[change[i]][k];
     a[change[i]][k] = ch;
  }
} else {
printf("I can't answer these equations!\n");
exit(-1);
}
}
tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is not 1 ***/
if(tmp!=1) {
notone[i] = tmp;
for(j=i; j<ma; j++) {
a[i][j] /= tmp;
}

}
tmp = a[i][i];
/*** the case of the  number of matrix[i][i] is 1 ***/
for(j=i+1; j<ma;j++) {
  pp = a[j][i] / tmp;
  divide[j][i] = pp;
  for(k=i;k<ma;k++) {
 a[j][k] -= pp * a[i][k];
  }
}
}
}
void lu_vector(double a[30][30], double b[30], double ret[30]) {
int i, j;
double ch, keika;
/**** forward elimination ***/
for(i=0; i<ma; i++) {
/***** the case of the number of matrix[i][i] is 0 ***/
if(change[i]!=0) {
  ch = b[i];
  b[i] = b[change[i]];
  b[change[i]] = ch;
}
/*** the case of the  number of matrix[i][i] is not 1 ***/
if(notone[i]!=0) {
  b[i] /= notone[i];
}
/*** the case of the  number of matrix[i][i] is 1 ***/
for(j=i+1; j<ma;j++) {
   b[j] -= divide[j][i] * b[i];
}
}
/**** backword elimination ****/
for(i=ma-1;i>=0;i--) {
keika=0.0;
for(j=i+1;j<ma;j++) {
   keika += a[i][j]*ret[j];
}

ret[i] = (b[i] - keika) / a[i][i];
}
}
void make_back(double c[30][30], double c_back[30][30]) {
int i, j;
double back[30];
lu_line(c);
for(j=0;j<ma;j++) {
for(i=0;i<ma;i++) {
   back[i] = 0.0;
}
back[j] = 1.0;
lu_vector(c,back,n);
for(i=0;i<ma;i++) {
c_back[i][j] = n[i];
}
}
}
void lu_apart(double line[30][30], double line_l[30][30], double line_u[30][30]) {
int i, j, k;
double total;
/***** first definition *****/
for(i=0;i<ma;i++) {
line_l[i][i] = 1.0;
for(j=i+1;j<ma;j++) {
   line_l[i][j] = 0.0;
}
for(j=0;j<i;j++) {
   line_u[i][j] = 0.0;
}
}

for(k=0;k<ma;k++) {      /***** line number *****/
for(i=k;i<ma;i++) {      /***** retu number *****/
   total = 0.0;
   for(j=0;j<k;j++) {      /***** teisuukou no genzan *****/
      total += line_l[k][j] * line_u[j][i];
   }

      line_u[k][i] = line[k][i] - total;
   }
   for(i=0;i<k+1;i++) {      /***** retu number *****/
      total = 0.0;
      for(j=0;j<i;j++) {
         total += line_l[k+1][j] * line_u[j][i];
      }
      line_l[k+1][i] = (line[k+1][i] - total) / line_u[i][i];
   }
}
}
void multi(double p[30][30], double q[30], double r[30]) {
int i, j;
for(i=0;i<ma;i++) {
   r[i] = 0.0;
}

for(i=0;i<ma;i++) {
   for(j=0;j<ma;j++) {
      r[i] += p[i][j] * q[j];
   }
}
}

int main() {
int i, j, k, re;
double time1, time2;
char X, Y;
printf("I can answer the equation by using LUbunkai!\n");
printf("Please select the equation size under 30.\n");
scanf("%d",&ma);
printf("%dequations!\n",ma);
init_lu();
printf("The x's dimention is %d\n",ma);
for(i=0;i<ma;i++) {
  result[i]=0.0;
}

X='X'; Y='Y';
print_lu(a,X,b);
for(i=0;i<ma;i++) {
  if(a[0][i]==0) {
    for(j=1;j<ma;j++) {
    if(a[j][i]!=0)  break;
   if(a[ma-1][i]==0) {
     printf("I can't answer.Because a group of a vertical numbers are 0!\n");
     exit(-1);
}
    }
  }
}
for(k=0;k<ma;k++) {
 if(a[k][0]!=0) {
   for(i=k+1;i<ma;i++) {
      re = (int)a[i][0] % (int)a[k][0];
      if(re==0) {
   time1 = a[i][0] / a[k][0];
for(j=0;j<ma;j++) {
if(a[k][j]!=0) {
           re = (int)a[i][j] % (int)a[k][j];
           if(re!=0) break;
time2 = a[i][j] / a[k][j];
       if(time1!=time2) break;
} else {
if(a[i][j]!=a[k][j]) break;
}
       if(j==ma-1) {
          printf("I can't answer.Because you type the same pair of matrixs!\n");
         exit(-1);
}
 }
      }
   }
  } else {
    for(i=k+1;i<ma;i++) {
       if(a[i][0]==0) {
         for(j=0;j<ma;j++) {
            if(a[i][j]!=a[k][j]) break;
   if(j==ma-1) {
          printf("I can't answer.Because you type the same pair of matrixs!\n");
     exit(-1);

   }
         }
       }
    }
  }
}
for(i=0;i<ma;i++) {
   change[i] = i;
   notone[i] = 0;
}
lu_apart(a,l,u);
printf("Next print is LUapart\n");
print_alu(a,l,u);
print_lu(l,Y,b);
make_back(l,l_back);
printf("Next print is L and L_back\n");
print_back(l,l_back);
multi(l_back,b,m);
printf("Y's answer is ----\n");
print_ans(Y,m);
print_lu(u,X,m);
make_back(u,u_back);
printf("Next print is U and U_back\n");
print_back(u,u_back);
multi(u_back,m,result);
printf("X's answer is ----\n");
print_answer();
return 0;
}

今回追加した重要な関数は特にない。ただ、行列の掛け算を行う関数 multi を用意したが、これについては、ただ単に行列の掛け算を行っているだけである。また、 LU 分解を行っている様子がわかりやすいように、いくつか出力を行う関数を追加してある。これについての説明は省略しておく。これを使って問題を解いた結果を下に示す。

2.17 実行結果

2.17.1 （１）

I can answer the equation by using LUbunkai!
Please select the equation size under 30.
3
3equations!
Please type 1 equation's 3 number!
1 2 3
Please type its vector's 1 number!
1
Please type 2 equation's 3 number!
2 8 11
Please type its vector's 1 number!
0
Please type 3 equation's 3 number!
3 22 35
Please type its vector's 1 number!
0
The x's dimention is 3
Equations are ...
 1.000000*X1 + 2.000000*X2 + 3.000000*X3 = 1.000000
 2.000000*X1 + 8.000000*X2 + 11.000000*X3 = 0.000000
 3.000000*X1 + 22.000000*X2 + 35.000000*X3 = 0.000000
Next print is LUapart
Your lines are ...
 1.000000  2.000000  3.000000
 2.000000  8.000000  11.000000
 3.000000  22.000000  35.000000
L------
 1.000000  0.000000  0.000000
 2.000000  1.000000  0.000000
 3.000000  4.000000  1.000000
U------
 1.000000  2.000000  3.000000
 0.000000  4.000000  5.000000
 0.000000  0.000000  6.000000
Equations are ...
 1.000000*Y1 + 2.000000*Y2 + 0.000000*Y3 = 1.000000
 2.000000*Y1 + 8.000000*Y2 + 0.000000*Y3 = 0.000000
 3.000000*Y1 + 22.000000*Y2 + 1.000000*Y3 = 0.000000
Next print is L and L_back

normal ----
1.000000  0.000000  0.000000
0.000000  1.000000  0.000000
0.000000  0.000000  1.000000
back  ----
1.000000  0.000000  0.000000
-2.000000  1.000000  0.000000
5.000000  -4.000000  1.000000
Y's answer is ----
Y1 = 1.000000
Y2 = -2.000000
Y3 = 5.000000
Equations are ...
1.000000*X1 + 2.000000*X2 + 3.000000*X3 = 1.000000
2.000000*X1 + 8.000000*X2 + 5.000000*X3 = -2.000000
3.000000*X1 + 22.000000*X2 + 6.000000*X3 = 5.000000
Next print is U and U_back
normal ----
1.000000  2.000000  3.000000
0.000000  1.000000  1.250000
0.000000  0.000000  1.000000
back  ----
1.000000  -0.500000  -0.083333
0.000000  0.250000  -0.208333
0.000000  0.000000  0.166667
X's answer is ----
***The answer***
X1 = 1.583333
X2 = -1.541667
X3 = 0.833333
Press any key to continue

、簡単に出力について説明しておくと、まず最初に、入力された方程式を出力する。ここでの X が、求める解となる。次に、行列の LU 分解を示している。その次に、 Ly=b の組み合わせの方程式を表示して、その後 L の逆行列を表示して、ｙの答えを表示する。そのあと、 Ux=y の組み合わせの方程式を表示して、 U の逆行列を表示して、最後に X の答えを表示している。ここで手計算も行ったので、その結果を下に示す。

X_{1} = \frac{19}{12}, X_{2} = −\frac{37}{24}, X_{3} = \frac{5}{6}

これより、実行結果の答えが正しい近似値であることがわかる。

2.17.2 （２）

I can answer the equation by using LUbunkai!

Please select the equation size under 30.
3
3equations!
Please type 1 equation's 3 number!
1 2 3
Please type its vector's 1 number!
1
Please type 2 equation's 3 number!
2 8 11
Please type its vector's 1 number!
2
Please type 3 equation's 3 number!
3 22 35
Please type its vector's 1 number!
3
The x's dimention is 3
Equations are ...
 1.000000*X1 + 2.000000*X2 + 3.000000*X3 = 1.000000
 2.000000*X1 + 8.000000*X2 + 11.000000*X3 = 2.000000
 3.000000*X1 + 22.000000*X2 + 35.000000*X3 = 3.000000
Next print is LUapart
Your lines are ...
 1.000000  2.000000  3.000000
 2.000000  8.000000  11.000000
 3.000000  22.000000  35.000000
L------
 1.000000  0.000000  0.000000
 2.000000  1.000000  0.000000
 3.000000  4.000000  1.000000
U------
 1.000000  2.000000  3.000000
 0.000000  4.000000  5.000000
 0.000000  0.000000  6.000000
Equations are ...
 1.000000*Y1 + 2.000000*Y2 + 0.000000*Y3 = 1.000000
 2.000000*Y1 + 8.000000*Y2 + 0.000000*Y3 = 2.000000
 3.000000*Y1 + 22.000000*Y2 + 1.000000*Y3 = 3.000000
Next print is L and L_back
normal ----
 1.000000  0.000000  0.000000
 0.000000  1.000000  0.000000
 0.000000  0.000000  1.000000

back  ----
 1.000000  0.000000  0.000000
 -2.000000  1.000000  0.000000
 5.000000  -4.000000  1.000000
Y's answer is ----
 Y1 = 1.000000
 Y2 = 0.000000
 Y3 = 0.000000
Equations are ...
 1.000000*X1 + 2.000000*X2 + 3.000000*X3 = 1.000000
 2.000000*X1 + 8.000000*X2 + 5.000000*X3 = 0.000000
 3.000000*X1 + 22.000000*X2 + 6.000000*X3 = 0.000000
Next print is U and U_back
normal ----
 1.000000  2.000000  3.000000
 0.000000  1.000000  1.250000
 0.000000  0.000000  1.000000
back  ----
 1.000000  -0.500000  -0.083333
 0.000000  0.250000  -0.208333
 0.000000  0.000000  0.166667
X's answer is ----
***The answer***
 X1 = 1.000000
 X2 = 0.000000
 X3 = 0.000000
Press any key to continue

手計算も行ったので、その結果を下に示す。

X_{1} = 1, X_{2} = 0, X_{3} = 0

これより、この答えは正しいことがわかる。

2.17.3 （３）

I can answer the equation by using LUbunkai!
Please select the equation size under 30.
3
3equations!
Please type 1 equation's 3 number!
1 2 3
Please type its vector's 1 number!
-1
Please type 2 equation's 3 number!

2 8 11
Please type its vector's 1 number!
0
Please type 3 equation's 3 number!
3 22 35
Please type its vector's 1 number!
1
The x's dimention is 3
Equations are ...
1.000000*X1 + 2.000000*X2 + 3.000000*X3 = -1.000000
2.000000*X1 + 8.000000*X2 + 11.000000*X3 = 0.000000
3.000000*X1 + 22.000000*X2 + 35.000000*X3 = 1.000000
Next print is LUapart
Your lines are ...
1.000000  2.000000  3.000000
2.000000  8.000000  11.000000
3.000000  22.000000  35.000000
L------
1.000000  0.000000  0.000000
2.000000  1.000000  0.000000
3.000000  4.000000  1.000000
U------
1.000000  2.000000  3.000000
0.000000  4.000000  5.000000
0.000000  0.000000  6.000000
Equations are ...
1.000000*Y1 + 2.000000*Y2 + 0.000000*Y3 = -1.000000
2.000000*Y1 + 8.000000*Y2 + 0.000000*Y3 = 0.000000
3.000000*Y1 + 22.000000*Y2 + 1.000000*Y3 = 1.000000
Next print is L and L_back
normal ----
1.000000  0.000000  0.000000
0.000000  1.000000  0.000000
0.000000  0.000000  1.000000
back  ----
1.000000  0.000000  0.000000
-2.000000  1.000000  0.000000
5.000000  -4.000000  1.000000
Y's answer is ----
Y1 = -1.000000
Y2 = 2.000000
Y3 = -4.000000

Equations are ...
 1.000000*X1 + 2.000000*X2 + 3.000000*X3 = -1.000000
 2.000000*X1 + 8.000000*X2 + 5.000000*X3 = 2.000000
 3.000000*X1 + 22.000000*X2 + 6.000000*X3 = -4.000000
Next print is U and U_back
normal ----
 1.000000  2.000000  3.000000
 0.000000  1.000000  1.250000
 0.000000  0.000000  1.000000
back  ----
 1.000000  -0.500000  -0.083333
 0.000000  0.250000  -0.208333
 0.000000  0.000000  0.166667
X's answer is ----
***The answer***
 X1 = -1.666667
 X2 = 1.333333
 X3 = -0.666667
Press any key to continue

手計算も行ったので、その結果を下に示す。

X_{1} = −\frac{5}{3}, X_{2} = \frac{4}{3}, X_{3} = −\frac{2}{3}

これより、答えは近似値であることが分かる。

2.18 性質３の証明

A を n × n 行列、 ||A|| を誘導されたノルム、 R を n × n 行列、 I を n × n 単位行列として、この時、 ||RA − I|| < 1 という条件である。 G = RA − I とすると、 ||G|| < 1 であり、この時、単位行列 I について、 I + G すなわち RA が正則行列である。よって、 A も正則である。ここで、 I + G が正則であるから、 I が単位行列であることより、次の式が成り立つ。

I = (I + G)^{−1}(I + G)

これより、

(I + G)^{−1} = I − (I + G)^{−1}G

性質 2 より、

||(I + G)^{−1}|| ≤ ||I|| + ||(I + G)^{−1}||||G||

これを整理して、

(1 − ||G||)(I + G)^{−1} ≤ 1

||G|| < 1 より、

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

となる。これを利用して、 ||A^{−1}|| について調べる。

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

となるので、 (1) より、

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

が成り立つ。

2.19 性質５の証明

行列 A の固有値 λ に対する固有ベクトル x を考える。 n × n 行列 B を適当に選び λ が B の固有値ではないときに、下の式が成り立つ。

(A − B)x = (λI − B)x

λ は B の固有値ではないから、

x = (λI − B)^{−1}(A − B)x

と書きかえられる。さらに、 ||Ax|| ≤ ||A||||x|| より、

||(λI − B)^{−1}(A − B)|| ≥ 1

が成り立つ。ここで、行列 A の体格セイブン以外を全て 0 とした行列を C として考えたとき、これを、 B = C とすると、最大値ノルムは、

||(λI − A_{D})^{−1}(A − A_{D}||_{∞} = max_{1≤n} \frac{1}{|λ − a_{ii}|} ∑_{k=1,k≠i}^{n} |a_{ik}| ≥ 1

これより、ある j が存在して、

|λ − a_{jj}| ≤ ∑_{k=1,(k≠j)}^{n} |a_{jk}|

が成り立つ。これより、ゲルシュゴリンの定理は証明された。

2.20 性質６の証明

IEEE754 の浮動小数点演算の定義から、 x, y ∈ F, · ∈ +, −, ×, / のとき、

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

となる。

x = (x_{1}, x_{2}, …, x_{n})^{t} ∈ R^{n}, x_{i} ∈ F,
y = (y_{1}, y_{2}, …, y_{n})^{t} ∈ R^{n}, y_{i} ∈ F

として話を進める。内積 s において、 s = x^{t}y = x_{1}y_{1} + … + x_{n}y_{n}、 T_{i} = x_{1}y_{1} + … + x_{i}y_{i}とすると、 T_{1} = ⊙(x_{1}y_{1}) = x_{1}y_{1}(1 + δ_{1}) T_{2} = ⊙(T_{1} + (x_{2}y_{2}) = (T_{1} + (x_{2}y_{2}(1 + δ_{2}))(1 + δ_{3}) = (x_{1}y_{1}(1 + δ_{1}) + x_{2}y_{2}(1 + δ_{2}))(1 + δ_{3})

= x_{1}y_{1}(1 + δ_{1})(1 + δ_{3}) + x_{2}y_{2}(1 + δ_{2})(1 + δ_{3}) これを一般化すると、 |δ_{i}| ≤ u で nu < 1 とする。このとき、

∏_{i=1}^{n} 1 + δ_{i} = 1 + θ_{n}

が成り立つ。ただし、条件は、

θ_{n} ≤ \frac{nu}{1 − nu} = γ_{n}

である。また適当な δ_{i} によってかける 1 + δ_{i} の n 個の積を θ_{n} とかくと、

T = x_{1}y_{1}(1 + θ_{n}) + x_{2}y_{2}(1 + θ′_{n}) + x_{3}y_{3}(1 + θ_{n−1}) + … + x_{n}y_{n}(1 + θ_{2})

と書けるので

|s − T | ≤ γ_{n} ∑_{i=1}^{n} |x_{i}y_{i}|

が成り立つ。以上により、性質６が証明された。

3 感想

今回の課題は、かなりきつかったです。教科書の問題については、プログラムを、部分に分けて少しずつ作っていったので、最後の５番の問題を解くときは、かなり気分よくできました。ループの山なので、エラーが出ずにコンパイルできても、出力される値がすごく変なものになるときもありましたが、終わってみると、それぞれ単純なことしかやっていなかったな、と思えます。教科書のように、順を追ってプログラムを作っていくと、結構面白かったです。次回水曜日までの課題もあるので、今回の課題はこの辺りで提出します。