fork download
copy 
  1. #include <stdio.h>
  2.  #include <math.h>
  3.  #include <stdlib.h>
  4.  #define TOL 1e-6
  5.  // 許容誤差
  6. #define MAX_ITER 1000 // 最大反復回数
  7. // 行列とベクトルの積を計算する関数
  8. void mat_vec_mult(int n, double A[n][n], double x[n], double result[n]) {
  9.  for (int i = 0; i < n; i++) {
  10.  result[i] = 0.0;
  11.  for (int j = 0; j < n; j++) {
  12.  result[i] += A[i][j] * x[j];
  13.  }
  14.  }
  15.  }
  16.  // ベクトルのノルムを計算する関数
  17. double norm(int n, double x[n]) {
  18.  double sum = 0.0;
  19.  for (int i = 0; i < n; i++) {
  20.  sum += x[i] * x[i];
  21.  }
  22.  return sqrt(sum);
  23.  }
  24.  // ベクトルを正規化する関数
  25. void normalize(int n, double x[n]) {
  26.  double x_norm = norm(n, x);
  27.  for (int i = 0; i < n; i++) {
  28.  x[i] /= x_norm;
  29.  }
  30.  }
  31.  // 累乗法による最大固有値と対応する固有ベクトルを計算する関数
  32. void power_method(int n, double A[n][n], double x0[n], double *eigenvalue, double eigenvector[n]) {
  33.  double x[n], y[n];
  34.  double lambda_old = 0.0, lambda_new = 0.0;
  35.  // 初期ベクトルをコピー
  36. for (int i = 0; i < n; i++) {
  37.  x[i] = x0[i];
  38.  }
  39.  // 正規化
  40. normalize(n, x);
  41.  for (int iter = 0; iter < MAX_ITER; iter++) {
  42.  // 行列とベクトルの積
  43. mat_vec_mult(n, A, x, y);
  44.  // 新しい固有値の近似値を計算（レイリー商）
  45. lambda_new = 0.0;
  46.  for (int i = 0; i < n; i++) {
  47.  lambda_new += x[i] * y[i];
  48.  }
  49.  // 収束判定
  50. if (fabs(lambda_new- lambda_old) < TOL) {
  51.  *eigenvalue = lambda_new;
  52.  for (int i = 0; i < n; i++) {
  53.  eigenvector[i] = y[i];
  54.  }
  55.  return;
  56.  }
  57.  // 正規化
  58. normalize(n, y);
  59. // 次の反復のために更新
  60.     lambda_old=lambda_new;
  61.  for(int i =0; i<n; i++) {
  62.  x[i] =y[i];
  63.  }
  64.  }
  65.  // 収束しない場合
  66.   fprintf(stderr, "累乗法が収束しませんでした。¥n");
  67.  exit(EXIT_FAILURE);
  68.  }
  69.  // メイン関数
  70. int main() {
  71.  // 行列サイズ
  72.   int n =4;
  73.  // 対象行列
  74.   double A[4][4] = {
  75.  {16, -1, 1, 2},
  76.  {2, 12, 1, -1},
  77.  {1, 3, 24, 2},
  78.  {4, -2, 1, 20},
  79.  };
  80. // 初期ベクトル
  81. double x0[4] = {1, 1, 1, 1};
  82.  // 結果を格納する変数
  83. double eigenvalue;
  84.  double eigenvector[2];
  85.  // 累乗法を実行
  86. power_method(n, A, x0, &eigenvalue, eigenvector);
  87.  // 正規化
  88. normalize(n, eigenvector);
  89.  // 結果を出力
  90. printf("最大固有値: %lf¥n", eigenvalue);
  91.  printf("対応する固有ベクトル: [");
  92.  for (int i = 0; i < n; i++) {
  93.  printf("%lf", eigenvector[i]);
  94.  if (i < n- 1) printf(", ");
  95.  }
  96.  printf("]¥n");
  97.  return 0;
  98.  }
Success #stdin #stdout 0s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
最大固有値: 25.013627¥n対応する固有ベクトル: [0.159495, 0.074987, 0.942110, 0.285249]¥n
https://ideone.com/YDXUVw
language:
C (gcc 8.3)
created:
3 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!