数値計算（大石先生）第三回課題

{学籍番号：ｇ９９ｐ１２７５ \atop 氏名：益子 理絵}

出題日：２００１年６月０１日提出期限：２００１年６月１５日

提出日２００１年６月１５日

1 課題内容

1.1 教科書ｐ３１の問３の１から５を解け。 1.2 授業でやった定理の性質 3,5,6 を証明せよ。（性質 1,2,4 は仮定してよい。）

2 課題：教科書の問題

2.1 教科書ｐ３１問３−１

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

2.1.1 作成したプログラム

1. ピボット選択なしのプログラムファイル名は、「gauss0.c」と「matrixSize.h」である。解きたい連立一次方程式が入力できるようになっている。変数の個数は２０個まで可であるが、もっと大きくしたい場合には、 matrixSize.h でサイズを変更すればよい。 gauss0.c の概要は以下のとおりである。

• main 関数（行列入力・前進消去・後退代入）を呼ぶ

• 行列入力関数

• 前進消去関数

• 後退代入関数

• 行列表示関数

(a) matrixSize.h

/* Header File for program: Gauss.c */
#define NMAX  20           /*変数の最大数*/

(b) gauss0.c

#include "matrixSize.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
float A[NMAX + 1][NMAX + 2];
int N;
void getmatrix();
void forward();
void back();
void Printmatrix();

main(){
  getmatrix();             /* 行列入力 */
  forward();               /* 前進消去 */
  back();                  /* 後退代入 */
}

void getmatrix() {
int i, j;
float t;
while(N<NMAX){

printf("変数の数はいくつですか? （最大 20 変数まで） ");

scanf("%d", &N);
printf("\n");
if(N>NMAX){
  printf("エラー：変数が多すぎます。 \n\n");
  N = 0;
}else{
  break;
}
}
for(i=1; i<=N; i++) {
for(j=1; j<=N; j++) {
  printf("A[%d%d] = ", i, j );
  scanf("%f", &t);
  A[i][j] = t;
}
printf("B[%d] = ", i);
scanf("%f", &t);
A[i][j] = t;
printf("\n");
}
printf("入力された値\n");
Printmatrix();
}
/*前進消去*/
void forward() {
int i, j, k;
float t;
for(j=1; j<N; j++) {
Printmatrix();
for(i=j+1; i<=N; i++) {
  t = A[i][j] / A[j][j];
  for(k=j; k<=N+1; k++) {
    A[i][k] -= t * A[j][k];
  }

  }
}
printf("前進消去後\n");
Printmatrix();
}
/*後退代入*/
void back() {
float s[NMAX], t;
int i, j;
printf("答（後退代入後） \n");
s[N] = A[N][N+1] / A[N][N];
printf("X[%d] = %lf\n",N ,s[N] );
for(i=N-1; i>=1; i--) {
  t = 0;
  for(j=i+1; j<=N; j++) {
    t += A[i][j] * s[j];
  }
  s[i] = ( A[i][N+1] - t ) / A[i][i];
  printf("X[%d] = %lf\n",i ,s[i] );
}
}
/*行列表示*/
void Printmatrix() {
int i, j;
for(i=1; i<=N; i++) {
  for(j=1; j<=N; j++) {
    printf("%10.6f ", A[i][j]);
  }
  printf("| %10.6f\n", A[i][j]);
}
printf("\n");
}

2. ピボット選択ありのプログラムファイル名は、「gauss.c」と「matrixSize.h」で添付しようかと思ったが統合性がなくなるので今回はやめた。 matrixSize.h は上と同じものを使っているので省略する。 gauss.c では、 gauss0.c にピボット選択を付け加えただけである。付け加えた内容は以下のとおりである。

void pivotchoose(int j) {
int i, k, imax;
float amax, T[NMAX + 1][NMAX + 2];
for(i=j; i<=N; i++) {
  imax = i;
  amax = fabs(A[i][j]);
  for(k=i; k<=N; k++){
    if(fabs(A[k][j])>amax) {
      imax = k;
      amax = fabs(A[k][j]);
    }
  }
  for(k=1; k<=N+1; k++) {
    T[i][k] = A[i][k];
    A[i][k] = A[imax][k];
    A[imax][k] = T[i][k];
  }
}
}

あとは、前進消去のループところで、 pivotchoose(j); を呼べばよい。

全体のプログラムは以下のようになる。

(a) gauss.c

#include "matrixSize.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
float A[NMAX + 1][NMAX + 2];
int N;
void getmatrix();
void pivotchoose(int j);
void forward();
void back();
void Printmatrix();

main(){
  getmatrix();             /* 行列入力 */
  forward();               /* 前進消去 */
  back();                  /* 後退代入 */
}

while(N<NMAX){

printf("変数の数はいくつですか? （最大 20 変数まで） ");

scanf("%d", &N);
printf("\n");
if(N>NMAX){
  printf("エラー：変数が多すぎます。 \n\n");
  N = 0;
}else{
  break;
}
}
for(i=1; i<=N; i++) {
for(j=1; j<=N; j++) {
  printf("A[%d%d] = ", i, j );
  scanf("%f", &t);
  A[i][j] = t;
}
printf("B[%d] = ", i);
scanf("%f", &t);
A[i][j] = t;
printf("\n");
}
printf("入力された値\n");
Printmatrix();
}

/*ピボット選択*/

void pivotchoose(int j) {
int i, k, imax;
float amax, T[NMAX + 1][NMAX + 2];
for(i=j; i<=N; i++) {
imax = i;
amax = fabs(A[i][j]);
for(k=i; k<=N; k++){
  if(fabs(A[k][j])>amax) {

      imax = k;
      amax = fabs(A[k][j]);
    }
  }
  for(k=1; k<=N+1; k++) {
    T[i][k] = A[i][k];
    A[i][k] = A[imax][k];
    A[imax][k] = T[i][k];
  }
}
}
/*前進消去*/
void forward() {
int i, j, k;
float t;
for(j=1; j<N; j++) {
  pivotchoose(j);
  Printmatrix();
  for(i=j+1; i<=N; i++) {
    t = A[i][j] / A[j][j];
    for(k=j; k<=N+1; k++) {
      A[i][k] -= t * A[j][k];
    }
  }
}
printf("前進消去後\n");
Printmatrix();
}
/*後退代入*/
void back() {
float s[NMAX], t;
int i, j;
printf("答（後退代入後） \n");
s[N] = A[N][N+1] / A[N][N];
printf("X[%d] = %lf\n",N ,s[N] );
for(i=N-1; i>=1; i--) {
  t = 0;
  for(j=i+1; j<=N; j++) {

    t += A[i][j] * s[j];
  }
  s[i] = ( A[i][N+1] - t ) / A[i][i];
  printf("X[%d] = %lf\n",i ,s[i] );
}
}
/*行列表示*/
void Printmatrix() {
int i, j;
for(i=1; i<=N; i++) {
  for(j=1; j<=N; j++) {
    printf("%10.6f ", A[i][j]);
  }
  printf("| %10.6f\n", A[i][j]);
}
printf("\n");
}

2.1.2 実行結果

1. ピボット選択なしのプログラム結果は以下のようになった。

(a) 実行結果：その１変数の数はいくつですか? （最大 20 変数まで） 3 入力された値

2.000000   3.000000  -3.000000 |  23.000000
-1.000000  -5.000000   2.000000 | -21.000000
3.000000  -9.000000   1.000000 | -20.000000
2.000000   3.000000  -3.000000 |  23.000000
-1.000000  -5.000000   2.000000 | -21.000000
3.000000  -9.000000   1.000000 | -20.000000
2.000000   3.000000  -3.000000 |  23.000000
0.000000  -3.500000   0.500000 |  -9.500000
0.000000 -13.500000   5.500000 | -54.500000

前進消去後

2.000000   3.000000  -3.000000 |  23.000000
0.000000  -3.500000   0.500000 |  -9.500000
0.000000   0.000000   3.571429 | -17.857141

答（後退代入後）

X[3] = -5.000000
X[2] = 2.000000
X[1] = 1.000001

(b) 実行結果：その２（ピボット選択しないと解けない連立方程式）変数の数はいくつですか? （最大 20 変数まで） 3 入力された値

 -1.000000   2.000000  -1.000000 |  16.000000
 -1.000000   2.000000  -3.000000 |  18.000000
  3.000000  -1.000000   4.000000 |  -4.000000
 -1.000000   2.000000  -1.000000 |  16.000000
 -1.000000   2.000000  -3.000000 |  18.000000
  3.000000  -1.000000   4.000000 |  -4.000000
 -1.000000   2.000000  -1.000000 |  16.000000
  0.000000   0.000000  -2.000000 |   2.000000
  0.000000   5.000000   1.000000 |  44.000000

前進消去後

 -1.000000   2.000000  -1.000000 |  16.000000
  0.000000   0.000000  -2.000000 |   2.000000
  0.000000  -1.#IND00   1.#INF00 |  -1.#INF00

答（後退代入後）

X[3] = -1.#IND00
X[2] = -1.#IND00
X[1] = -1.#IND00

やはり０除算が行われるので、答えがエラーになっている。

2. ピボット選択ありのプログラム色々入れる値を変えて試してみたが、プログラム自体は正しく動いている。２つの実行結果の例を以下に示す。計算過程を示すようになっているのでわかりやすいと思う。

• その１

変数の数はいくつですか? （最大 20 変数まで） 3 入力された値

    2.000000   4.000000   6.000000 |   6.000000
    3.000000   8.000000   7.000000 |  15.000000
    5.000000   7.000000  21.000000 |  24.000000
    5.000000   7.000000  21.000000 |  24.000000
    3.000000   8.000000   7.000000 |  15.000000

2.000000   4.000000   6.000000 |   6.000000
5.000000   7.000000  21.000000 |  24.000000
-0.000000   3.800000  -5.600000 |   0.599999
-0.000000   1.200000  -2.400000 |  -3.600000

前進消去後

5.000000   7.000000  21.000000 |  24.000000
-0.000000   3.800000  -5.600000 |   0.599999
-0.000000  -0.000000  -0.631579 |  -3.789474

答（後退代入後）

X[3] = 6.000001
X[2] = 9.000003
X[1] = -33.000008

• その２

変数の数はいくつですか? （最大 20 変数まで） 5 入力された値

1.000000  -2.000000   3.000000   1.000000   2.000000 |  13.000000
2.000000   1.000000  -1.000000   2.000000   4.000000 |   7.000000
3.000000  -1.000000  -2.000000   1.000000  -1.000000 |  -5.000000
1.000000   3.000000   1.000000  -4.000000  -2.000000 |   2.000000
4.000000  -2.000000   1.000000   3.000000  -3.000000 |  -7.000000
4.000000  -2.000000   1.000000   3.000000  -3.000000 |  -7.000000
3.000000  -1.000000  -2.000000   1.000000  -1.000000 |  -5.000000
2.000000   1.000000  -1.000000   2.000000   4.000000 |   7.000000
1.000000   3.000000   1.000000  -4.000000  -2.000000 |   2.000000
1.000000  -2.000000   3.000000   1.000000   2.000000 |  13.000000
4.000000  -2.000000   1.000000   3.000000  -3.000000 |  -7.000000
0.000000   3.500000   0.750000  -4.750000  -1.250000 |   3.750000
0.000000   2.000000  -1.500000   0.500000   5.500000 |  10.500000
0.000000  -1.500000   2.750000   0.250000   2.750000 |  14.750000
0.000000   0.500000  -2.750000  -1.250000   1.250000 |   0.250000
4.000000  -2.000000   1.000000   3.000000  -3.000000 |  -7.000000
0.000000   3.500000   0.750000  -4.750000  -1.250000 |   3.750000
0.000000   0.000000   3.071429  -1.785714   2.214286 |  16.357143
0.000000  -0.000000  -2.857143  -0.571429   1.428571 |  -0.285714
0.000000  -0.000000  -1.928571   3.214286   6.214286 |   8.357142
4.000000  -2.000000   1.000000   3.000000  -3.000000 |  -7.000000

     0.000000   3.500000   0.750000  -4.750000  -1.250000 |   3.750000
     0.000000   0.000000   3.071429  -1.785714   2.214286 |  16.357143
     0.000000  -0.000000  -0.000000  -2.232558   3.488372 |  14.930234
     0.000000  -0.000000  -0.000000   2.093023   7.604651 |  18.627907

前進消去後

     4.000000  -2.000000   1.000000   3.000000  -3.000000 |  -7.000000
     0.000000   3.500000   0.750000  -4.750000  -1.250000 |   3.750000
     0.000000   0.000000   3.071429  -1.785714   2.214286 |  16.357143
     0.000000  -0.000000  -0.000000  -2.232558   3.488372 |  14.930234
     0.000000  -0.000000  -0.000000  -0.000000  10.875000 |  32.625000

答（後退代入後）

   X[5] = 3.000000
   X[4] = -2.000000
   X[3] = 2.000000
   X[2] = -1.000001
   X[1] = 1.000000

2.1.3 考察（ガウスの消去法）

• 前進消去と後退代入方程式 Ax = b の係数行列 A を、対角要素の下側を順次消去することで上三角行列に変換し、 (前進消去) 次に一番下側の方程式から解を求め、順次それを代入し次々と解を求めていく (後退代入) 方法である。（この過程については、詳しく柏木先生のレポートでやったので省略します。）

• ピボット選択係数の割る値のことをピボット (pivot,、軸) という。これは 0 であってはならず、かつ丸め誤差の影響を小さくするためには、なるべく絶対値の大きい方がよい。このため、ピボットのところに絶対値の大きな係数をもってくるという操作を行なう。これがピボット選択である。ピボット選択には、単に方程式（すなわち、行）の順序を入れ換える部分選択法 (partial pivoting) と変数（すなわち、列）の順序も入れ換える完全選択法 (complete pivoting) がある。部分選択法に要する処理時間は、方程式を解くための全体の時間に比較してごくわずかである。完全選択法は部分選択法に比べ、精度の面では有利であるが、処理時間は消去法の正味の計算時間と同程度であり、全体の計算時間は部分選択法に比べて約２倍になる。また、解の順番を元通りに並べ換える後処理も必要となる。このピボット選択をすると、ガウスの消去法は浮動小数点計算などの丸め誤差に対して安定となる。それでも、入力値によっては正しい値が得られない（わずかにずれた値になる）ことがある。

• 計算量ガウスの消去法のうち、前進消去のプロセス全体においては、ｎ３／３＋ O （ｎ２）回の乗除算が行われている。また、後退代入のプロセス全体では、ｎ２／２＋ O （ｎ）回の乗除算が行われている。よって、前進消去と合わせると、ガウスの消去アルゴリズムでは、ｎ３／３＋ O （ｎ２）回の乗除算で解が求められる。ガウス・ジョルダン（Gauss-Jordan）法というアルゴリズムの必要とする演算回数 (同様に各約ｎ３／２回) より約３割少ないのである。

2.2 教科書ｐ３１問３−２

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

2.2.1 問題の解法・アルゴリズム

前進消去が終わると、左辺は上三角行列になっている。ここまで変形ができれば、最終行から順に x を求めることができる。これが後退代入であった。前進消去の際の計算は、行列の演算としては行基本変形である。この変換を行列 L^{−1}、出来上がった上三角行列を U で表せば、 L^{−1}A = U となる。ここで L は対角要素が全て 1 である下三角行列になる。すなわち、前進消去とは行列 A を LU に分解することと等価であり、 Ax = b は LUx = b と書き直せる。すると、 Ly = b、 Ux = y の 2 本の方程式を解くことと等価である。まず、係数行列 A を上三角行列 U に変換する。その変換仕様が L^{−1} である。次に 1 本目の方程式を y について解く。この式は右辺にも同じ変換を掛けることを示している (∵ y=L^{−1}b)。次に 2 本目の式を x について解く。これらの方程式はどちらも三角方程式であるから、簡単に解を求めることができる。これが LU 分解法であり、以下のプログラムではこのアルゴリズムを用いて方程式を解く。

2.2.2 作成したプログラム

• lu.c

#include<stdio.h>
#include<stdlib.h>
#include<float.h>
#include<math.h>
#define N 4
void decomp(double a[N][N], double m[N][N], int e[N], int n);
void solve(double m[N][N], double b[N], double x[N], int e[N], int n);
void print_matrix(double a[N][N], int n);
void print_vector(double a[N], int n);
void exchange_vector(double a[N], int e[N], int n);
int init_matrix(double a[N][N], double b[N]);
int compare(const void *big, const void *small);
int main(){
 int i,j,k;
 int n;
 double a[N][N];
 double b[N];
 double m[N][N];

double x[N];
int e[N];
double big, small;
while(1){
n=init_matrix(a,b);
if(n==0){
   exit(0);
}
for(i=0;i<n;i++){
   e[i]=i;
}
printf("initial matrix:\n");
printf("a:\n");
print_matrix(a,n);
printf("b:\n");
print_vector(b,n);
printf("\n");
decomp(a,m,e,n);
printf("LU decomposed matrix:\n");
print_matrix(m,n);
printf("\n");
solve(m,b,x,e,n);
printf("solution x:\n");
print_vector(x,n);
printf("\n");
}
}
void decomp(double a[N][N], double m[N][N], int e[N], int n){
int i,j,k,p,max;
double temp[N];
for(k=0;k<n-1;k++){
/* cheque the maximum coefficient of rows */
   max=k;
   for(p=k+1;p<n;p++){
       if(fabs(a[k][k]) < fabs(a[p][k]))
       max=p;
   }
/* cheque coefficient */
/* if the maximum number is smaller than DBL_EPSILON, */
/* the coefficient matrix is singlular, cannot solve. */
if(fabs(a[max][k])<=DBL_EPSILON){
   printf("an error occurred: coefficient matrix is singular.\n");
   exit(1);
}

/* exchange rows */
for(p=0;p<n;p++){
    temp[p]=a[max][p];
}
for(p=0;p<n;p++){
    a[max][p]=a[k][p];
}
for(p=0;p<n;p++){
    a[k][p]=temp[p];
}
/* exchange index vector */
p=e[max];
e[max]=e[k];
e[k]=p;
for(i=k+1;i<n;i++){
/* calculate multiplier */
    m[i][k]=-a[i][k]/a[k][k];
/* forward elimination */
    for(j=k+1;j<n;j++){
        a[i][j]+=m[i][k]*a[k][j];
    }
}
}
/* coefficient of a[n-1][n-1] must be grater than DBL_EPSILON */
/* if it is smaller than DBL_EPSILON, coefficient matrix is */
/* singular. */
if(fabs(a[n-1][n-1])<=DBL_EPSILON){
    printf("an error occurred: coefficient matrix is singular.\n");
    exit(1);
}
/* a[][]:=U */
/* m[][]:=-L */
/* marge them into m[][] */
for(i=0;i<n;i++){
    for(j=0;j<n;j++){
        m[i][j]=-m[i][j];
    }
}
for(i=0;i<n;i++){
    for(j=i;j<n;j++){
        m[i][j]=a[i][j];
    }
}
}

void solve(double m[N][N], double b[N], double x[N], int e[N], int n){
int i,j,k;
double sum=0;
double y[N];
/* exchange b[] by e[] */
exchange_vector(b,e,n);
/* solve Ly=b */
/* L is lower triangular matrix */
y[0]=b[0];
for(i=1;i<n;i++){
    double temp[N];
    sum=0;
    for(j=0;j<i;j++){
        temp[j]=m[i][j]*y[j];
    }
qsort(temp,i,sizeof(double),compare);
for(j=i-1;j>=0;j--){
    sum+=temp[j];
}
y[i]=b[i]-sum;
}
/* solve Ux=y */
/* U is upper triangular matrix */
for(i=n-1;i>=0;i--){
    double temp[N];
    sum=0;
for(j=i+1;j<n;j++){
    temp[j]=m[i][j]*x[j];
}
qsort(temp,i,sizeof(double),compare);
    for(j=n-1;j>=i+1;j--){
        sum+=temp[j];
    }
x[i]=(y[i]-sum)/m[i][i];
}
}
void exchange_vector(double a[N], int e[N], int n){
int i,j;
double temp[N];
for(i=0;i<n;i++){
temp[i]=a[i];
}
for(i=0;i<n;i++){
    a[i]=temp[e[i]];

}
}
void print_matrix(double a[N][N], int n){
int i,j;
for(i=0;i<n;i++){
    printf("|");
    printf(" ");
for(j=0;j<n;j++){
    printf("% 10.5f",a[i][j]);
    printf(" ");
}
printf("|");
printf("\n");
}
printf("\n");
}
void print_vector(double a[N], int n){
int j;
printf("|");
printf(" ");
for(j=0;j<n;j++){
    printf("% 10.5f",a[j]);
    printf(" ");
}
printf("|");
printf("\n");
printf("\n");
}
int init_matrix(double a[N][N], double b[N]){
int i,j,k,in,n;
printf(" 0:self input\n");
printf("1-6:pre-installed matrix\n? ");
scanf("%d",&in);
if(in==0){
    printf("enter degree of the matrix:\n? ");
    scanf("%d",&n);
    printf("enter coefficient matrix:\n? ");
    for(i=0;i<n;i++){
        for(j=0;j<n;j++,k++){
            scanf("%lf",&a[i][j]);
        }
    }
    printf("enter right hand side matrix:\n? ");
    for(i=0;i<n;i++){

    scanf("%lf",&b[i]);
}
return n;
}
if(in==1){
double temp[N*N+N]={ 33.0, 16.0, 72.0,\
    -24.0, -10.0, -57.0,\
    -8.0, -4.0, -17.0,\
    -359.0, 281.0, 85.0};
n=3;
k=0;
for(i=0;i<n;i++){
    for(j=0;j<n;j++,k++){
        a[i][j]=temp[k];
    }
}
for(i=0;i<n;i++,k++){
    b[i]=temp[k];
}
return n;
}
if(in==2){
double temp[N*N+N]={ 1.0, -2.0, 3.0, 1.0,\
    -2.0, 1.0, -2.0, -1.0,\
    3.0, -2.0, 1.0, 5.0,\
1.0, -1.0, 5.0, 3.0,\
3.0, -4.0, 7.0, 8.0};
n=4;
k=0;
for(i=0;i<n;i++){
for(j=0;j<n;j++,k++){
    a[i][j]=temp[k];
}
}
for(i=0;i<n;i++,k++){
    b[i]=temp[k];
}
return n;
}
if(in==3){
double temp[N*N+N]={ 5.0, 7.0, 6.0, 5.0,\
7.0, 10.0, 8.0, 7.0,\
6.0, 8.0, 10.0, 9.0,\
5.0, 7.0, 9.0, 10.0,\

23.0, 32.0, 33.0, 31.0};
n=4;
k=0;
for(i=0;i<n;i++){
for(j=0;j<n;j++,k++){
a[i][j]=temp[k];
}
}
for(i=0;i<n;i++,k++){
    b[i]=temp[k];
}
return n;
}
(in==4){
double temp[N*N+N]={ 0.4096, 0.1234, 0.3678, 0.2943,\
0.2246, 0.3872, 0.4015, 0.1129,\
0.3645, 0.1920, 0.3781, 0.0643,\
0.1784, 0.4002, 0.2786, 0.3927,\
0.4043, 0.1550, 0.4240,-0.2557};
n=4;
k=0;
for(i=0;i<n;i++){
for(j=0;j<n;j++,k++){
a[i][j]=temp[k];
}
}
for(i=0;i<n;i++,k++){
b[i]=temp[k];
}
return n;
}
if(in==5){
double temp[N*N+N]={ 2.4, 1.5, 10.1,\
3.6, 4.3, 5.6,\
4.8, 4.6, 7.7,\
51.02, 62.99, 78.31};
n=3;
k=0;
for(i=0;i<n;i++){
    for(j=0;j<n;j++,k++){
        a[i][j]=temp[k];
    }
}
for(i=0;i<n;i++,k++){

b[i]=temp[k];
}
return n;
}
if(in==5){
double temp[N*N+N]={ 2.4, 1.5, 10.1,\
3.6, 4.3, 5.6,\
4.8, 4.6, 7.7,\
51.02, 62.99, 78.31};
n=3;
k=0;
for(i=0;i<n;i++){
    for(j=0;j<n;j++,k++){
        a[i][j]=temp[k];
    }
}
for(i=0;i<n;i++,k++){
    b[i]=temp[k];
}
return n;
}
if(in==6){
double temp[N*N+N]={ 12.0, 38.0, 56.0,\
18.0, 57.0, 67.0,\
24.0, 26.0, 35.0,\
412.2, 575.8, 396.9};
n=3;
k=0;
for(i=0;i<n;i++){
    for(j=0;j<n;j++,k++){
        a[i][j]=temp[k];
    }
}
for(i=0;i<n;i++,k++){
b[i]=temp[k];
}
return n;
}
if(in==7){
double temp[N*N+N]={ 5.0, 10.0, 13.0,\
10.0, 20.0, 26.0,\
24.0, 26.0, 35.0,\
3.0, 6.0, 22.0};
n=3;

k=0;
for(i=0;i<n;i++){
    for(j=0;j<n;j++,k++){
        a[i][j]=temp[k];
    }
}
for(i=0;i<n;i++,k++){
    b[i]=temp[k];
}
return n;
}
if(in==8){
double temp[N*N+N]={ 0.01, 0.15, 0.21,\
1000, 2000, 2600,\
12500, 33000, 25000,\
3.0, 66.5, 1020};
n=3;
k=0;
for(i=0;i<n;i++){
    for(j=0;j<n;j++,k++){
        a[i][j]=temp[k];
    }
}
for(i=0;i<n;i++,k++){
    b[i]=temp[k];
}
return n;
}
printf("an error occurred: inputed number was not defined.\n");
exit(2);
 }
int compare(const void *big, const void *small){
double b=*(double*)big;
double s=*(double*)small;
if(fabs(b)>fabs(s)){
return -1;
}
if(fabs(b)==fabs(s)){
return 0;
}
if(fabs(b)<fabs(s)){
return 1;}
}

• プログラムの説明プログラムが長くなったので簡単に説明する。

1. 主な変数 int n 行列の次数 double a[N] [N] 係数行列 double b[N] 右辺の行列 double m[N] [N] L- E+U を格納する配列（E:単位行列） double x[N] 解を格納する配列 e[ N] 行の入れ換えを行った結果を格納する配列

2. プリプロセッサヘッダ各 include 命令と、行列の最大次数を書き込んでおく。課題として与えられている行列の次元は最大で 4 なので、この値で各配列を確保するものとする。

3. 関数 main()

使用する変数をここで全て宣言しておく。今回各関数間で受け渡す変数は、すべて配列すなわちポインタであるから、あとは関数を呼びだす際に引数にこれらを与えればよい。次に係数行列を init_matrix() を呼び出して初期化する。それから初期状態を表示したあと、 decomp() を呼び出して LU 分解を行い、その結果を表示する。そして最後に solve() で解を求めて表示する。

4. 関数 decomp()

引数： double a[ N] [ N] ,double mn] [ N] ,int e[ N] ,int n 戻り値：なしこの関数で LU 分解を行う。まず、部分枢軸法に基づいてピボットを選ぶ。数値計算を行うにあたり、丸め誤差を減らすため、係数の絶対値がもっとも大きい行を探して入れ換えを行う（結果を e[ ] に格納）。ここでもし 0 でないか、マシンイプシロンよりも大きな係数を持つ行が存在しない場合、係数行列 A は正則ではなく、解は存在しないので、エラーメッセージを表示したあとプログラムを終了する。そのあと、乗数を計算して m[ ] [ ] に代入し、前進消去を行う。このような処理を行うと、 a[ ] [ ] は上三角行列になり、 m[ ] [ ] には下三角行列が格納される。そこで最後にこれらを m[ ] [ ] にまとめ上げる。

5. 関数 solve()

引数： double m[ N] [ N] ,double b[ N] ,double x[ N] ,int e[ N] ,int n 戻り値：なし分解された行列から、解を出す。まず、 decomp() で行われた入れ換えを右辺の係数列ベクトル b[ ] に反映させる。次に Ly =b を解く。 L は m[ ] [ ] の下三角部分に格納されている下三角行列であるから、第 1 行から順に解を求めていく。次に Ux =y を解く。 U は m[ ] [ ] の上三角部分であるから、こちらもやはり簡単に求めることができる。

6. 関数 exchange_vect or()

引数： double a[ N] ,int e[ N] ,int n 戻り値：なし引数として受取った a[ ] を、インデックスベクトル e[ ] を使って入れ換える。

7. 関数 print _matrix()

引数： double a[N] [N] ,int n 戻り値：なし引数 a[ ] [ ] を n × n 行列として表示する。

      8. 関数 print_vector()

引数： double a[N] ,int n

戻り値：なし

引数 a[ ] を n 行の列ベクトルとして表示する。

      9. 関数 init_matrix()

引数： double a[N] [N] ,double b[N]

戻り値：行列の次数 int n

予め与えておく 6 つの係数行列を格納し、かつユーザからの入力を受け付けて初期化を行う。

2.2.3 実行結果

３通りについての結果を示す。 LU の行列式を途中で出力している。

その１

a:
|    1.00000   -2.00000    3.00000    1.00000 |
|   -2.00000    1.00000   -2.00000   -1.00000 |
|    3.00000   -2.00000    1.00000    5.00000 |
|    1.00000   -1.00000    5.00000    3.00000 |
b:
|    3.00000   -4.00000    7.00000    8.00000 |
LU decomposed matrix:
|    3.00000   -2.00000    1.00000    5.00000 |
|   -0.66667   -0.33333    4.66667    1.33333 |
|    0.33333    4.00000  -16.00000   -6.00000 |
|    0.33333    1.00000    0.37500    3.25000 |
solution x:
|    5.30769    5.84615    3.15385   -0.07692 |

その２

a:
|    5.00000    7.00000    6.00000    5.00000 |
|    7.00000   10.00000    8.00000    7.00000 |
|    6.00000    8.00000   10.00000    9.00000 |
|    5.00000    7.00000    9.00000   10.00000 |
b:
|   23.00000   32.00000   33.00000   31.00000 |
LU decomposed matrix:
|    6.00000    8.00000   10.00000    9.00000 |
|    1.16667    0.66667   -3.66667   -3.50000 |
|    0.83333    0.50000    2.50000    4.25000 |
|    0.83333    0.50000   -0.20000    0.10000 |
solution x:
|    1.00000    1.00000    1.00000    1.00000 |

その３

a:
|    2.40000    1.50000   10.10000 |
|    3.60000    4.30000    5.60000 |
|    4.80000    4.60000    7.70000 |
b:
|   51.02000   62.99000   78.31000 |
LU decomposed matrix:
|    4.80000    4.60000    7.70000 |
|    0.75000    0.85000   -0.17500 |
|    0.50000   -0.94118    6.08529 |
solution x:
|    6.81575    5.54582    2.60826 |

プログラムは期待通りに動いている。手計算での期待した値とは違う答え、すなわち誤差があるのは承知だが、近似値としてはまずまずだろう。

2.2.4 考察

LU 分解法を用いるに当っての問題は、 (i) 連立方程式そのものが持つ問題 (ii) 数値計算を行うことで起こる問題の二つが存在する。前者には更に、 (i) 係数行列が正則ではない (ii) 係数行列が悪い線形系であるという二つの問題が存在する。 (i) の場合は一意な解（非自明解）は存在しない。 (ii) の場合、残差を小さく抑えることができても誤差が大きくなってしまう。これらはプログラミングによって改良するのは無理である。一方後者には (i) 丸め誤差が出る (ii) 三角方程式を出す際の加算で桁落ちが起こるという問題がある。 (i) については、それを避けるべくピボット選択（部分枢軸法）を導入した。また、 (ii) については、絶対値の小さい順に加えていくことで回避することとした。それと、ガウスの消去法で作成したプログラムと同じ式を入力して、その結果の値を比較しようと思ったが、既に遅かったのでレポート出した後に時間ができたらやってみようと思います。

2.3 教科書ｐ３１問３−３

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

2.3.1 問題の解法・アルゴリズム

行列 A の逆行列 X = A^{−1} を計算する。ここで、 X = [x_{1}, x_{2}, x_{3}, …, x_{n}] とする。ただし、 x_{i} は逆行列 X の第 i 列を表す n 次元ベクトルとする。 e_{i} を第 i 成分が１でその他の成分が０の n 次元ベクトルとするとき、 Ax_{i} = e_{i} であるから、 A の LU 分解を求め、上式を、 i = 1, 2, …, n と n 回とけば、逆行列が求められることになる。その乗除算の回数はｎ^{3} + Ｏ（n^{2}）である。このアルゴリズムを使ってプログラムを作成する。

2.3.2 作成したプログラム

#include <stdio.h>

#define MAX  10
/*decomp*/
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;

  }
}
          printf("**** triangular 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("**** triangular matrix U****\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];
}

}
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("------Inverse matrix A------ \n");
            printf("Input size of A(max:10) \n size:");scanf("%d",&n);
            printf("\n");
            if(n>MAX){
                printf(":too many parameters:\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("Input:\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("**** inverse matrix A****\n");
     for(i=1; i<=n; i++){
          for(j=1; j<=n; j++){
                 printf("%10f ",A[i][j]);
}
          printf("\n");
}
      return 0;
}

2.3.3 実行結果

• その１

------Inverse matrix A------
Input size of A(max:10)
 size:2
Input:
a[11] = 1
a[12] = 0
a[21] = -1
a[22] = 1
**** triangular matrix L****
  1.000000   0.000000
 -1.000000   1.000000
**** triangular matrix U****
  1.000000   0.000000
  0.000000   1.000000

**** inverse matrix A****
  1.000000   0.000000
  1.000000   1.000000

• その２

------Inverse matrix A------
Input size of A(max:10)
 size:3
Input:
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
**** triangular matrix L****
  1.000000   0.000000   0.000000
  2.000000   1.000000   0.000000
  3.000000   4.000000   1.000000
**** triangular matrix U****
  1.000000   2.000000   3.000000
  0.000000   4.000000   5.000000
  0.000000   0.000000   6.000000
**** inverse matrix A****
  1.583333  -0.166667  -0.083333
 -1.541667   1.083333  -0.208333
  0.833333  -0.666667   0.166667

• その３

------Inverse matrix A------
Input size of A(max:10)
 size:4
Input:
a[11] = 3

a[12] = 2
a[13] = 4
a[14] = 5
a[21] = 7
a[22] = -2
a[23] = 3
a[24] = 3
a[31] = 1
a[32] = 9
a[33] = 5
a[34] = 4
a[41] = 3
a[42] = 3
a[43] = 8
a[44] = 8
**** triangular matrix L****
  1.000000   0.000000   0.000000   0.000000
  2.333333   1.000000   0.000000   0.000000
  0.333333  -1.250000   1.000000   0.000000
  1.000000  -0.150000  -0.717647   1.000000
**** triangular matrix U****
  3.000000   2.000000   4.000000   5.000000
  0.000000  -6.666667  -6.333333  -8.666667
  0.000000   0.000000  -4.250000  -8.500000
  0.000000   0.000000   0.000000  -4.400000
**** inverse matrix A****
  0.066845   0.163102   0.066845  -0.136364
  0.125668  -0.013369   0.125668  -0.136364
 -0.909091   0.181818   0.090909   0.454545
  0.836898  -0.237968  -0.163102  -0.227273

2.3.4 考察

すべてを手計算でやることはめんどくさくってやらなかったのだが、その１だけはやってみた。結果は合っていたのでプログラムは正しく動いている。使用したアルゴリズムについては上で述べたので省略する。

2.4 教科書ｐ３１問３−４

次の行列 A の LU 分解を求めよ。

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

2.4.1 解 1. ２行目−２×（１行目）と３行目−３×（１行目）

( \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 16 & 26 \end{matrix} )

2. 上の結果の３行目−４×２行目

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

よって、

U = ( \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{matrix} ) ・・・・・・（答）

となる。また、 1 での変形過程より、行列 L の要素 l_{21} = \frac{2}{1}, l_{31} = \frac{3}{1} となる。 2 の過程より、行列 L の要素 l_{32} = \frac{16}{4} となる。よって、行列 L は、

L = ( \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & 1 \end{matrix} ) ・・・・・・（答）

となる。

2.4.2 たしかめ

LU = ( \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & 1 \end{matrix} ) ( \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{matrix} ) = ( \begin{matrix}1 & 2 & 3 \\ 2 & 8 & 11 \\ 3 & 22 & 35 \end{matrix} ) = A

ゆえに、 A = LU となっている。

2.5 教科書ｐ３１問３−５

問３−４の結果を用いて、

(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 を解け。 2.5.1 問題の解法・アルゴリズム

問題３−４より、 Ax = LUx = b なので、 y = ux としてまず Ly = b を解いて y を求める。続いて、 Ux = y から x を求める。ここで、問３−４の結果から、

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

であるのでこれを使用して解く。 Ly = b と Ux = y のそれぞれの計算は L が下三角行列、 U が上三角行列なので簡単に解ける。

2.5.2 解答

1. （１）

b = ( \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ) のとき

• Ly = b

( \begin{matrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 3 & 4 & 1 \end{matrix} ) ( \begin{matrix}y_{0} \\ y_{1} \\ y_{2} \end{matrix} ) = ( \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ) から、( \begin{matrix}y_{0} = 1 \\ 2y_{0} + y_{1} = 0 \\ 3y_{0} + 4y_{1} + y_{2} = 0 \end{matrix} ) となるので、

順に下式を上式から係数倍引いて計算すると、 y_{0} = 1, y_{1} = −2, y_{2} = 5 すなわち、

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

となる。

• Ux = y

( \begin{matrix}1 & 2 & 3 \\ 0 & 4 & 5 \\ 0 & 0 & 6 \end{matrix} ) ( \begin{matrix}x_{0} \\ x_{1} \\ x_{2} \end{matrix} ) = ( \begin{matrix}1 \\ −2 \\ 5 \end{matrix} ) より、( \begin{matrix}x_{0} + 2x_{1} + 3x_{2} = 1 \\ 4x_{1} + 5x_{2} = −2 \\ 6x_{2} = 5 \end{matrix} )

となるので、今度は順に下式から代入して解いていく。すると、

x_{2} = \frac{5}{6}, x_{1} = −\frac{37}{24}, x_{0} = \frac{19}{12}すなわち、 x = ( −\frac{37}{24}_{\frac{5}{6}}^{\frac{19}{12}} ) ・・・・・（答） である。

2. （２）

b = ( \begin{matrix}1 \\ 2 \\ 3 \end{matrix} ) のとき

（１）と同様にして計算する。連立方程式の左辺は以下すべて等しくなるので省略する。もともと、 LU 分解法は Ax = b を解くにあたって、右辺が変わった場合にもう一度最初からときなおさずに b の部分だけ計算しなおせばいいというアルゴリズムなのである。問題から左辺が同じなのは当たり前である。 • Ly = b

同様にして、 y = ( \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ) となる。

• Ux = y

同様にして、 x = ( \begin{matrix}1 \\ 0 \\ 0 \end{matrix} ) ・・・・・・（答え） となる。

3. （３）

b = ( \begin{matrix}−1 \\ 0 \\ 1 \end{matrix} ) のとき

これも上の二つと同様にして計算する。 • Ly = b

同様にして、 y = ( \begin{matrix}−1 \\ 2 \\ 4 \end{matrix} ) となる。

• Ux = y

同様にして、 x = ( −\frac{5}{3} ; \frac{4}{3} ; −\frac{2}{3} ) ・・・・・・（答え） となる。

3 課題：授業でやった定理の性質の証明

3.0.3 性質３

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

∥ RA − I ∥< 1

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

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

が成立する。

2. 証明 ∥ RA − I ∥< 1 のとき、 I + (RA − I) の逆行列が存在することが示せる。 G = RA − I とおく。 ∥ G ∥< 1 のとき、 I + G の逆行列が存在するということは以下のように証明される。今、 I+G が特異であるとすると、あるベクトル x ≠ 0 が存在して (I + G)x = 0 を満たす。 ∥ x ∥=1 であるように規格化されているとしても一般性を失わない。このとき次の式が成立する。

Gx = −x, ∥ x ∥= 1

これは、 ∥ Gx ∥=∥ x ∥= 1 から、 ∥ G ∥≥ 1 を意味している。しかし、これは仮定と矛盾しているので、 I + G は正則である。次にこのことをふまえて考えると、次の式が成り立つ。

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

この式を書き換えると (I + G)^{−1} = I − (I + G)^{−1}G となる。よって、以下の式を得る。

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

よって、

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

となる。このとき、 ∥ G ∥< 1 なので、以下のようになる。

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

このことを元に次の評価を行うと定義が証明されることになる。

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

3.0.4 性質５

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

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

の中に存在する。

2. 証明 x が行列 A の固有値 λ に対する固有ベクトルであるとする。今、 n × n 行列 B を任意に選び、 λ は B の個有値でないとすると、次の式が成立する。

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

λ が B の固有値でないことから (λI − B)^{−1} が存在することがわかるので次のようになる。

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

これより、 ∥ (λI − B)^{−1}(A − B) ∥≥ 1 が得られる。ここで、特に行列 A の対角成分以外をすべて零とした行列を A_{D} として、 B = A_{D} の場合を考える。ノルムとして、最大ノルムを考えると、以下の式を得る事ができる。

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

これより、あるｊが存在して次の式が成立し、ゲルシュゴリンの定理が証明されたことになる。

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

3.0.5 性質６

1. 性質６ IEEE754 の下での浮動小数点演算の下で議論すると次のことが成り立つ。 x,y∈F,· ∈ {+,-,*,/} とすると、 x ∘ y = (x · y)(1 + δ) で | δ |≤ u となる。このとき u = 2^{−53} とする。

2. 証明 x ∈ R が F の範囲内にあれば、

∘x = x(1 + δ), | δ |< u

となる。ただし、 u は単位の相対丸め（unit roundoff）といい、 u = 2^{−N} で与えられる。倍制度浮動小数点数においては、 N ＝５３であるから、 u = 2^{−53} ≤ 1.12 × 10^{−16} となる。実際、 x ∈ R が、 F の範囲内であれば、

x_{1} ≤ x ≤ x_{1} + 2^{e−N} = x_{2}

と書ける。但し、

x_{1} = (\frac{1}{2} + \frac{d_{2}}{2^{2}} + \frac{d_{3}}{2^{3}} + … + \frac{d_{N}^{1}}{2^{N}}) × 2^{e}

とする。ここで、一般性を失うことなく、 x ≤ 0 と仮定した。よって、 |\frac{∘x−x}{x}| ≤ \frac{x_{2}−x_{1}}{2} × \frac{1}{2^{e}−1} ≤ \frac{2^{e−N−1}}{2^{e−1}} = 2^{−N} = u を得る。よって、 IEEE754 における不動小数点数演算の定義から、明らかに x, y ∈ F, · ∈ {+,-,*,/} のとき、

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

となる。

4 感想

今週は発表と他のレポートが８つから１０ほど重なっていたので本当に疲れました。しかも、家で証明のところは tex で打って学校にいってコンパイルしてみたら、エラーが１００個くらいでたので、あきらめて打ち直したので悲しかったです。毎日がどれかの課題の締め切り日というのが多くて、へとへとだった。電車でも立ったまま寝れるくらいでふらふらします。必死でやるから、わかんなくてもむりやりわからせるというかとりあえずわかれって感じではっぱかけているので時間ができたときにゆっくり見直したいです。たかが連立一次方程式でも５つくらいからになると手計算でやるのはしんどいっていうことがよくわかりました。計算機にやらせた結果と比較するために自分で解いてみたのだが、えらくめんどくさかったです。ガウスの消去法でのプログラムは、柏木先生の課題で作成したものを使いました。

参考文献

[1] C によるアルゴリズム辞典（奥村晴彦著）

[2] 大石先生の書いた本と教科書（数値計算）

[3] 理工系ユーザーの為の UNIX 入門（武藤覚著）

[4] 基礎 C 言語（土屋範久著）