fork download
copy 
  1.  #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define SIZE 4 // 行列のサイズ
  5. #define MAX_ITER 1000 // 最大反復回数
  6. #define TOL 1e-6 // 収束判定の許容誤差
  7.  
  8. // QR分解を行うための関数（グラムシュミット直交化を使用）
  9. void qr_decompose(double A[SIZE][SIZE], double Q[SIZE][SIZE], double R[SIZE][SIZE]) {
  10.     int i, j, k;
  11.     double norm;
  12.  
  13.     for (i = 0; i < SIZE; i++) {
  14.         for (j = 0; j < SIZE; j++) {
  15.             Q[i][j] = 0.0;
  16.             R[i][j] = 0.0;
  17.         }
  18.     }
  19.  
  20.     for (k = 0; k < SIZE; k++) {
  21.         for (i = 0; i < SIZE; i++) {
  22.             Q[i][k] = A[i][k];
  23.         }
  24.         for (j = 0; j < k; j++) {
  25.             R[j][k] = 0.0;
  26.             for (i = 0; i < SIZE; i++) {
  27.                 R[j][k] += Q[i][j] * A[i][k];
  28.             }
  29.             for (i = 0; i < SIZE; i++) {
  30.                 Q[i][k] -= R[j][k] * Q[i][j];
  31.             }
  32.         }
  33.         norm = 0.0;
  34.         for (i = 0; i < SIZE; i++) {
  35.             norm += Q[i][k] * Q[i][k];
  36.         }
  37.         R[k][k] = sqrt(norm);
  38.         for (i = 0; i < SIZE; i++) {
  39.             Q[i][k] /= R[k][k];
  40.         }
  41.     }
  42. }
  43.  
  44. // QR法で固有値を計算する関数
  45. int qr_algorithm(double A[SIZE][SIZE], double eigenvalues[SIZE]) {
  46.     double Q[SIZE][SIZE], R[SIZE][SIZE], A_next[SIZE][SIZE];
  47.     int i, j, k;
  48.     int iteration_count = 0;
  49.  
  50.     for (k = 0; k < MAX_ITER; k++) {
  51.         iteration_count++;
  52.         // QR分解
  53.         qr_decompose(A, Q, R);
  54.         // 新しい行列 A = R * Q を計算
  55.         for (i = 0; i < SIZE; i++) {
  56.             for (j = 0; j < SIZE; j++) {
  57.                 A_next[i][j] = 0.0;
  58.                 for (int l = 0; l < SIZE; l++) {
  59.                     A_next[i][j] += R[i][l] * Q[l][j];
  60.                 }
  61.             }
  62.         }
  63.         // 収束判定
  64.         int converged = 1;
  65.         for (i = 0; i < SIZE; i++) {
  66.             if (fabs(A[i][i] - A_next[i][i]) > TOL) {
  67.                 converged = 0;
  68.                 break;
  69.             }
  70.         }
  71.         if (converged) break;
  72.         // 更新
  73.         for (i = 0; i < SIZE; i++) {
  74.             for (j = 0; j < SIZE; j++) {
  75.                 A[i][j] = A_next[i][j];
  76.             }
  77.         }
  78.     }
  79.     // 固有値を対角成分として取得
  80.     for (i = 0; i < SIZE; i++) {
  81.         eigenvalues[i] = A[i][i];
  82.     }
  83.     return iteration_count;
  84. }
  85.  
  86. // メイン関数
  87. int main() {
  88.     double A[SIZE][SIZE] = {
  89.         {16, -1, 1, 2},
  90.         {2, 12, 1, -1},
  91.         {1, 3, 24, 2},
  92.         {4, -2, 1, 20}
  93.     };
  94.     double eigenvalues[SIZE];
  95.  
  96.     printf("初期行列:\n");
  97.     for (int i = 0; i < SIZE; i++) {
  98.         for (int j = 0; j < SIZE; j++) {
  99.             printf("%8.4f ", A[i][j]);
  100.         }
  101.         printf("\n");
  102.     }
  103.  
  104.     int iteration_count = qr_algorithm(A, eigenvalues);
  105.  
  106.     printf("\n固有値:\n");
  107.     for (int i = 0; i < SIZE; i++) {
  108.         printf("lambda[%d] = %8.4f\n", i, eigenvalues[i]);
  109.     }
  110.  
  111.     printf("\n反復回数: %d\n", iteration_count);
  112.     return 0;
  113. }
  114.  
Success #stdin #stdout 0s 5284KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
初期行列:
 16.0000  -1.0000   1.0000   2.0000 
  2.0000  12.0000   1.0000  -1.0000 
  1.0000   3.0000  24.0000   2.0000 
  4.0000  -2.0000   1.0000  20.0000 

固有値:
lambda[0] =  25.0136
lambda[1] =  20.8385
lambda[2] =  14.1274
lambda[3] =  12.0204

反復回数: 80
https://ideone.com/pevkMd
反復法、固有値
language:
C (gcc 8.3)
created:
12 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!