fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. #define SIZE 4
  6. #define EPSILON 1e-9   // 収束判定条件
  7. #define MAX_ITER 1000  // 最大反復回数
  8.  
  9. // ベクトルのノルム（ユークリッド長）を計算
  10. double vector_norm(double *v, int n) {
  11.     double sum = 0.0;
  12.     for (int i = 0; i < n; i++) {
  13.         sum += v[i] * v[i];
  14.     }
  15.     return sqrt(sum);
  16. }
  17.  
  18. // 行列とベクトルの積 y = A * x
  19. void mat_vec_mul(double A[SIZE][SIZE], double *x, double *y, int n) {
  20.     for (int i = 0; i < n; i++) {
  21.         y[i] = 0.0;
  22.         for (int j = 0; j < n; j++) {
  23.             y[i] += A[i][j] * x[j];
  24.         }
  25.     }
  26. }
  27.  
  28. int main() {
  29.     // 自身で設定した行列A
  30.     double A[SIZE][SIZE] = {
  31.         {1.0,  2.0, 3.0,  5.0},
  32.         {3.0, -5.0, 1.0,  4.0},
  33.         {5.0,  9.0, 2.0, -6.0},
  34.         {1.0,  7.0, 4.0,  1.0}
  35.     };
  36.  
  37.     double x[SIZE], y[SIZE], Ax[SIZE];
  38.     double lambda_old = 0.0, lambda_new = 0.0;
  39.  
  40.     // 1. 初期値の設定 [1, 1, 1, 1]^T
  41.     for (int i = 0; i < SIZE; i++) {
  42.         x[i] = 1.0;
  43.     }
  44.  
  45.     printf("=== べき乗法による収束状況の確認 ===\n");
  46.     printf("初期ベクトル x(0): [%.1f, %.1f, %.1f, %.1f]\n\n", x[0], x[1], x[2], x[3]);
  47.     printf("%-5s | %-12s | %-12s | %-30s\n", "Iter", "固有値(新)", "差分値", "固有ベクトルx(k)");
  48.     printf("----------------------------------------------------------------------------\n");
  49.  
  50.     // べき乗法の反復ループ
  51.     int actual_iters = 0;
  52.     for (int iter = 1; iter <= MAX_ITER; iter++) {
  53.         actual_iters = iter;
  54.  
  55.         // y = A * x
  56.         mat_vec_mul(A, x, y, SIZE);
  57.  
  58.         // ノルムで正規化
  59.         double norm_y = vector_norm(y, SIZE);
  60.         if (norm_y < 1e-12) break;
  61.         for (int i = 0; i < SIZE; i++) {
  62.             y[i] /= norm_y;
  63.         }
  64.  
  65.         // レイリー商を用いて固有値を高精度に推定
  66.         mat_vec_mul(A, y, Ax, SIZE);
  67.         double num = 0.0, den = 0.0;
  68.         for (int i = 0; i < SIZE; i++) {
  69.             num += y[i] * Ax[i];
  70.             den += y[i] * y[i];
  71.         }
  72.         lambda_new = num / den;
  73.  
  74.         // 前のステップとの変化量を追跡
  75.         double diff = fabs(lambda_new - lambda_old);
  76.  
  77.         // 収束状況の出力表示 (10回ごと、および最終ステップ)
  78.         if (iter <= 5 || iter % 5 == 0) {
  79.             printf("%-5d | %-12.8f | %-12.4e | [%.5f, %.5f, %.5f, %.5f]\n", 
  80.                    iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
  81.         }
  82.  
  83.         // 収束判定条件のチェック
  84.         if (diff < EPSILON) {
  85.             if (iter % 5 != 0 && iter > 5) { // 最終ステップを確実に表示させる
  86.                 printf("%-5d | %-12.8f | %-12.4e | [%.5f, %.5f, %.5f, %.5f]\n", 
  87.                        iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
  88.             }
  89.             break;
  90.         }
  91.  
  92.         // 次のステップに向けて x を更新
  93.         for (int i = 0; i < SIZE; i++) {
  94.             x[i] = y[i];
  95.         }
  96.         lambda_old = lambda_new;
  97.     }
  98.  
  99.     printf("----------------------------------------------------------------------------\n");
  100.     printf("反復回数 %d 回で収束しました。\n", actual_iters);
  101.  
  102.     // === 計算結果の表示 ===
  103.     printf("\n=== 計算結果 ===\n");
  104.     printf("最大固有値 (λ1)          : %.10f\n", lambda_new);
  105.     printf("対応する固有ベクトル (x1) : [");
  106.     for (int i = 0; i < SIZE; i++) {
  107.         printf("%.10f%s", y[i], (i == SIZE - 1) ? "]\n" : ", ");
  108.     }
  109.  
  110.     // === 確かに固有値・固有ベクトルであることの検証 (Ax = λx) ===
  111.     printf("\n=== 検証: Ax と λx の比較 ===\n");
  112.     mat_vec_mul(A, y, Ax, SIZE); // Axを再計算
  113.     printf("Ax    : [");
  114.     for (int i = 0; i < SIZE; i++) {
  115.         printf("%.6f%s", Ax[i], (i == SIZE - 1) ? "]\n" : ", ");
  116.     }
  117.     printf("λ * x : [");
  118.     for (int i = 0; i < SIZE; i++) {
  119.         printf("%.6f%s", lambda_new * y[i], (i == SIZE - 1) ? "]\n" : ", ");
  120.     }
  121.  
  122.     return 0;
  123. }
  124.  
Success #stdin #stdout 0.01s 5316KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=== べき乗法による収束状況の確認 ===
初期ベクトル x(0): [1.0, 1.0, 1.0, 1.0]

Iter  | 固有値(新) | 差分値    | 固有ベクトルx(k)        
----------------------------------------------------------------------------
1     | 7.06015038   | 7.0602e+00   | [0.55069, 0.15019, 0.50063, 0.65081]
2     | 7.56272200   | 5.0257e-01   | [0.68485, 0.48918, 0.14675, 0.51975]
3     | 7.79665790   | 2.3394e-01   | [0.53345, 0.20811, 0.56744, 0.59172]
4     | 8.36088120   | 5.6422e-01   | [0.66241, 0.41253, 0.25086, 0.57281]
5     | 8.52180634   | 1.6093e-01   | [0.58882, 0.28456, 0.47179, 0.59139]
10    | 8.69758623   | 9.4221e-03   | [0.62524, 0.34192, 0.38099, 0.58908]
15    | 8.69591999   | 7.0277e-04   | [0.62205, 0.33644, 0.39006, 0.58967]
20    | 8.69624282   | 7.2845e-05   | [0.62237, 0.33699, 0.38916, 0.58962]
25    | 8.69621204   | 7.2757e-06   | [0.62234, 0.33693, 0.38925, 0.58962]
30    | 8.69621514   | 7.2934e-07   | [0.62235, 0.33694, 0.38924, 0.58962]
35    | 8.69621483   | 7.3086e-08   | [0.62234, 0.33694, 0.38924, 0.58962]
40    | 8.69621486   | 7.3240e-09   | [0.62234, 0.33694, 0.38924, 0.58962]
45    | 8.69621486   | 7.3395e-10   | [0.62234, 0.33694, 0.38924, 0.58962]
----------------------------------------------------------------------------
反復回数 45 回で収束しました。

=== 計算結果 ===
最大固有値 (λ1)          : 8.6962148610
対応する固有ベクトル (x1) : [0.6223447229, 0.3369369049, 0.3892384556, 0.5896219066]

=== 検証: Ax と λx の比較 ===
Ax    : [5.412043, 2.930076, 3.384901, 5.127479]
λ * x : [5.412043, 2.930076, 3.384901, 5.127479]
https://ideone.com/R6CKgM
language:
C (gcc 8.3)
created:
1 day 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!