fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. #define SIZE 4
  6. #define EPSILON 1e-7   // 収束判定条件
  7. #define MAX_ITER 1000  // 最大反復回数
  8.  
  9. // 行列とベクトルの積 y = A * x
  10. void mat_vec_mul(double A[SIZE][SIZE], double *x, double *y, int n) {
  11.     for (int i = 0; i < n; i++) {
  12.         y[i] = 0.0;
  13.         for (int j = 0; j < n; j++) {
  14.             y[i] += A[i][j] * x[j];
  15.         }
  16.     }
  17. }
  18.  
  19. int main() {
  20.     // 自身で設定した行列A
  21.     double A[SIZE][SIZE] = {
  22.         {1.0,  2.0, 3.0,  5.0},
  23.         {3.0, -5.0, 1.0,  4.0},
  24.         {5.0,  9.0, 2.0, -6.0},
  25.         {1.0,  7.0, 4.0,  1.0}
  26.     };
  27.  
  28.     double x[SIZE], y[SIZE], Ax[SIZE];
  29.     double lambda_old = 0.0, lambda_new = 0.0;
  30.  
  31.     // 初期値の設定 [1, 1, 1, 1]^T
  32.     for (int i = 0; i < SIZE; i++) {
  33.         x[i] = 1.0;
  34.     }
  35.  
  36.     printf("=== べき乗法の計算および収束状況 ===\n");
  37.     printf("初期ベクトル x(0): [%.1f, %.1f, %.1f, %.1f]\n\n", x[0], x[1], x[2], x[3]);
  38.     printf("%-5s | %-14s | %-14s | %-30s\n", "Iter", "固有値の近似値", "変化量(差分)", "固有ベクトル x(k)");
  39.     printf("----------------------------------------------------------------------------\n");
  40.  
  41.     int actual_iters = 0;
  42.     for (int iter = 1; iter <= MAX_ITER; iter++) {
  43.         actual_iters = iter;
  44.  
  45.         // y = A * x の計算
  46.         mat_vec_mul(A, x, y, SIZE);
  47.  
  48.         // スライドのアルゴリズムに従い、第1成分を基準（係数）として抽出
  49.         lambda_new = y[0]; 
  50.  
  51.         // 第1成分が 1.0 になるように全体を正規化
  52.         if (fabs(lambda_new) < 1e-12) {
  53.             break;
  54.         }
  55.         for (int i = 0; i < SIZE; i++) {
  56.             y[i] /= lambda_new;
  57.         }
  58.  
  59.         // 前のステップからの変化量を算出
  60.         double diff = fabs(lambda_new - lambda_old);
  61.  
  62.         // 収束過程の表示 (最初の5回、5回ごと、および最終ステップ)
  63.         if (iter <= 5 || iter % 5 == 0) {
  64.             printf("%-5d | %-14.10f | %-14.4e | [%.6f, %.6f, %.6f, %.6f]\n", 
  65.                    iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
  66.         }
  67.  
  68.         // 収束判定条件のチェック
  69.         if (diff < EPSILON) {
  70.             if (iter % 5 != 0 && iter > 5) {
  71.                 printf("%-5d | %-14.10f | %-14.4e | [%.6f, %.6f, %.6f, %.6f]\n", 
  72.                        iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
  73.             }
  74.             break;
  75.         }
  76.  
  77.         // 次の反復のためにベクトルと値を更新
  78.         for (int i = 0; i < SIZE; i++) {
  79.             x[i] = y[i];
  80.         }
  81.         lambda_old = lambda_new;
  82.     }
  83.  
  84.     printf("----------------------------------------------------------------------------\n");
  85.     printf("反復回数 %d 回で収束判定条件を満たしました。\n", actual_iters);
  86.  
  87.     // === 最終計算結果の出力 ===
  88.     printf("\n=== 最終結果 ===\n");
  89.     printf("最大固有値 (λ1)          : %.10f\n", lambda_new);
  90.     printf("対応する固有ベクトル (x1) : [");
  91.     for (int i = 0; i < SIZE; i++) {
  92.         printf("%.10f%s", y[i], (i == SIZE - 1) ? "]\n" : ", ");
  93.     }
  94.  
  95.     // === Ax = λx の検証出力 ===
  96.     printf("\n=== 固有値・固有ベクトルであることの検証 (Ax = λx) ===\n");
  97.     mat_vec_mul(A, y, Ax, SIZE); 
  98.     printf("Ax    : [");
  99.     for (int i = 0; i < SIZE; i++) {
  100.         printf("%.6f%s", Ax[i], (i == SIZE - 1) ? "]\n" : ", ");
  101.     }
  102.     printf("λ * x : [");
  103.     for (int i = 0; i < SIZE; i++) {
  104.         printf("%.6f%s", lambda_new * y[i], (i == SIZE - 1) ? "]\n" : ", ");
  105.     }
  106.  
  107.     return 0;
  108. }
Success #stdin #stdout 0.01s 5332KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
=== べき乗法の計算および収束状況 ===
初期ベクトル x(0): [1.0, 1.0, 1.0, 1.0]

Iter  | 固有値の近似値 | 変化量(差分) | 固有ベクトル x(k)       
----------------------------------------------------------------------------
1     | 11.0000000000  | 1.1000e+01     | [1.000000, 0.272727, 0.909091, 1.181818]
2     | 10.1818181818  | 8.1818e-01     | [1.000000, 0.714286, 0.214286, 0.758929]
3     | 6.8660714286   | 3.3157e+00     | [1.000000, 0.390117, 1.063719, 1.109233]
4     | 10.5175552666  | 3.6515e+00     | [1.000000, 0.622774, 0.378709, 0.864738]
5     | 7.7053659743   | 2.8122e+00     | [1.000000, 0.483272, 0.801255, 1.004364]
10    | 8.7981535734   | 2.6057e-01     | [1.000000, 0.546856, 0.609357, 0.942157]
15    | 8.6860605584   | 2.6271e-02     | [1.000000, 0.540849, 0.627060, 0.947951]
20    | 8.6972331750   | 2.6313e-03     | [1.000000, 0.541454, 0.625276, 0.947367]
25    | 8.6961128229   | 2.6369e-04     | [1.000000, 0.541394, 0.625455, 0.947425]
30    | 8.6962250867   | 2.6425e-05     | [1.000000, 0.541400, 0.625437, 0.947420]
35    | 8.6962138365   | 2.6481e-06     | [1.000000, 0.541399, 0.625439, 0.947420]
40    | 8.6962149639   | 2.6536e-07     | [1.000000, 0.541399, 0.625439, 0.947420]
43    | 8.6962148354   | 6.6741e-08     | [1.000000, 0.541399, 0.625439, 0.947420]
----------------------------------------------------------------------------
反復回数 43 回で収束判定条件を満たしました。

=== 最終結果 ===
最大固有値 (λ1)          : 8.6962148354
対応する固有ベクトル (x1) : [1.0000000000, 0.5413991505, 0.6254386721, 0.9474201120]

=== 固有値・固有ベクトルであることの検証 (Ax = λx) ===
Ax    : [8.696215, 4.708123, 5.438949, 8.238969]
λ * x : [8.696215, 4.708123, 5.438949, 8.238969]
https://ideone.com/oDXDuk
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!