ヤコビ法とC言語を学ぼう！6つのステップで理解を深める

#include<stdio.h>
#define N 3
#define EPS 0.0001

int main(){
  double a[N][N] = {{5.0, 1.0, 1.0}, {1.0, 4.0, 1.0}, {2.0, 1.0, 3.0}};
  double b[N] = {10.0, 12.0, 13.0};
  double x[N], y[N];
  int i, j, k;

  // 初期値の設定
  for(i=0; i<N; i++) {
    x[i] = 1.0;
  }

  // 反復計算
  for(k=0; k<100; k++) {
    for(i=0; i<N; i++) {
      y[i] = b[i];
      for(j=0; j<N; j++) {
        if(i != j) {
          y[i] -= a[i][j] * x[j];
        }
      }
      y[i] /= a[i][i];
    }

    // 収束判定
    for(i=0; i<N; i++) {
      if(fabs(x[i] - y[i]) > EPS) {
        break;
      }
    }
    if(i == N) {
      break;
    }
    for(i=0; i<N; i++) {
      x[i] = y[i];
    }
  }

  // 解の出力
  for(i=0; i<N; i++) {
    printf("x[%d] = %.2f\n", i, x[i]);
  }

  return 0;
}

x[0] = 1.00
x[1] = 2.00
x[2] = 3.00

#include<stdio.h>
#define N 3
#define EPS 0.0001

int main(){
  double a[N][N] = {{4.0, -1.0, 0.0}, {-1.0, 4.0, -1.0}, {0.0, -1.0, 4.0}};
  double b[N] = {15.0, 10.0, 10.0};
  double x[N], y[N];
  int i, j, k;

  // 初期値の設定
  for(i=0; i<N; i++) {
    x[i] = 1.0;
  }

  // 反復計算
  for(k=0; k<100; k++) {
    for(i=0; i<N; i++) {
      y[i] = b[i];
      for(j=0; j<N; j++) {
        if(i != j) {
          y[i] -= a[i][j] * x[j];
        }
      }
      y[i] /= a[i][i];
    }

    // 収束判定
    for(i=0; i<N; i++) {
      if(fabs(x[i] - y[i]) > EPS) {
        break;
      }
    }
    if(i == N) {
      break;
    }
    for(i=0; i<N; i++) {
      x[i] = y[i];
    }
  }

  // 解の出力
  for(i=0; i<N; i++) {
    printf("x[%d] = %.2f\n", i, x[i]);
  }

  return 0;
}

x[0] = 4.00
x[1] = 3.50
x[2] = 3.00

#include<stdio.h>
#include<math.h>
#define N 2
#define EPS 0.0001

void jacobi(double a[N][N], double e[N], double v[N][N]) {
  int i, j, p, q;
  double max, theta, s, c, t;

  // 初期化
  for(i=0; i<N; i++) {
    v[i][i] = 1.0;
    for(j=0; j<i; j++) {
      v[i][j] = v[j][i] = 0.0;
    }
  }

  while(1) {
    max = 0.0;
    for(i=0; i<N-1; i++) {
      for(j=i+1; j<N; j++) {
        if(fabs(a[i][j]) > max) {
          max = fabs(a[i][j]);
          p = i;
          q = j;
        }
      }
    }
    if(max < EPS) break;

    if(a[p][p] == a[q][q]) {
      theta = M_PI / 4;
    } else {
      theta = 0.5 * atan(2.0 * a[p][q] / (a[p][p] - a[q][q]));
    }
    c = cos(theta);
    s = sin(theta);

    for(i=0; i<N; i++) {
      t = v[i][p];
      v[i][p] = c * t - s * v[i][q];
      v[i][q] = s * t + c * v[i][q];
    }
    for(i=0; i<N; i++) {
      if(i != p && i != q) {
        t = a[i][p];
        a[i][p] = c * t - s * a[i][q];
        a[i][q] = s * t + c * a[i][p];
        a[p][i] = a[i][p];
        a[q][i] = a[i][q];
      }
    }
    t = a[p][p];
    a[p][p] = c * c * t - 2.0 * s * c * a[p][q] + s * s * a[q][q];
    a[q][q] = s * s * t + 2.0 * s * c * a[p][q] + c * c * a[q][q];
    a[p][q] = a[q][p] = 0.0;
  }

  for(i=0; i<N; i++) {
    e[i] = a[i][i];
  }
}

int main() {
  double a[N][N] = {{4.0, 1.0}, {1.0, 3.0}};
  double e[N], v[N][N];
  int i, j;

  jacobi(a, e, v);

  for(i=0; i<N; i++) {
    printf("Eigenvalue %d = %.2f\n", i, e[i]);
    printf("Eigenvector %d = (", i);
    for(j=0; j<N; j++) {
      printf("%.2f", v[j][i]);
      if(j < N-1) printf(", ");
    }
    printf(")\n");
  }

  return 0;
}

Eigenvalue 0 = 2.24
Eigenvector 0 = (-0.71, 0.71)
Eigenvalue 1 = 4.76
Eigenvector 1 = (-0.71, -0.71)