fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. #define N 4             // Matrix dimension
  6. #define MAX_ITER 100    // Maximum number of sweeps
  7. #define EPSILON 1e-6    // 収束判定条件 (Convergence threshold)
  8.  
  9. // Function prototypes
  10. void jacobi_method(double A[N][N], double eigenvalues[N], double eigenvectors[N][N]);
  11. void verify_results(double orig_A[N][N], double eigenvalues[N], double eigenvectors[N][N]);
  12.  
  13. int main() {
  14.     // 1. Define the matrix from image_4.png
  15.     double A[N][N] = {
  16.         {5.0, 4.0, 1.0, 1.0},
  17.         {4.0, 5.0, 1.0, 1.0},
  18.         {1.0, 1.0, 4.0, 2.0},
  19.         {1.0, 1.0, 2.0, 4.0}
  20.     };
  21.  
  22.     // Keep a copy of the original matrix for validation later
  23.     double orig_A[N][N];
  24.     for (int i = 0; i < N; i++) {
  25.         for (int j = 0; j < N; j++) {
  26.             orig_A[i][j] = A[i][j];
  27.         }
  28.     }
  29.  
  30.     double eigenvalues[N];
  31.     double eigenvectors[N][N];
  32.  
  33.     // 2. Run Jacobi Method (Prints iteration history)
  34.     jacobi_method(A, eigenvalues, eigenvectors);
  35.  
  36.     // 3. Print Final Results
  37.     printf("\n--- Final Results ---\n");
  38.     for (int i = 0; i < N; i++) {
  39.         printf("Eigenvalue %d (固有値): %9.6f\n", i + 1, eigenvalues[i]);
  40.         printf("Eigenvector %d (固有ベクトル): [ ", i + 1);
  41.         for (int j = 0; j < N; j++) {
  42.             printf("%9.6f ", eigenvectors[j][i]); // Columns are eigenvectors
  43.         }
  44.         printf("]\n\n");
  45.     }
  46.  
  47.     // 4. Verification Check (固有値・固有ベクトルの確認)
  48.     verify_results(orig_A, eigenvalues, eigenvectors);
  49.  
  50.     return 0;
  51. }
  52.  
  53. void jacobi_method(double A[N][N], double eigenvalues[N], double eigenvectors[N][N]) {
  54.     // Initialize eigenvector matrix as identity matrix
  55.     for (int i = 0; i < N; i++) {
  56.         for (int j = 0; j < N; j++) {
  57.             eigenvectors[i][j] = (i == j) ? 1.0 : 0.0;
  58.         }
  59.     }
  60.  
  61.     printf("--- Iteration Process (収束状況) ---\n");
  62.  
  63.     int iter_count = 0;
  64.     for (int iter = 1; iter <= MAX_ITER; iter++) {
  65.         double max_off_diag = 0.0;
  66.         int p = 0, q = 0;
  67.  
  68.         // Find the largest off-diagonal element |A[p][q]|
  69.         for (int i = 0; i < N; i++) {
  70.             for (int j = i + 1; j < N; j++) {
  71.                 if (fabs(A[i][j]) > max_off_diag) {
  72.                     max_off_diag = fabs(A[i][j]);
  73.                     p = i;
  74.                     q = j;
  75.                 }
  76.             }
  77.         }
  78.  
  79.         iter_count++;
  80.  
  81.         // Print the iteration count and the maximum off-diagonal value remaining
  82.         printf("%d: Max off-diagonal element = %e\n", iter_count, max_off_diag);
  83.  
  84.         // Check for convergence
  85.         if (max_off_diag < EPSILON) {
  86.             printf("\n[Convergence Achieved]\n");
  87.             break;
  88.         }
  89.  
  90.         // --- Pure Algebraic Calculation of cos(theta) and sin(theta) ---
  91.         double app = A[p][p];
  92.         double aqq = A[q][q];
  93.         double apq = A[p][q];
  94.  
  95.         double c, s; 
  96.  
  97.         if (fabs(apq) < 1e-12) {
  98.             c = 1.0;
  99.             s = 0.0;
  100.         } else {
  101.             double phi = 0.5 * (app - aqq) / apq;
  102.             double t;
  103.             if (phi >= 0.0) {
  104.                 t = 1.0 / (phi + sqrt(phi * phi + 1.0));
  105.             } else {
  106.                 t = -1.0 / (-phi + sqrt(phi * phi + 1.0));
  107.             }
  108.             c = 1.0 / sqrt(1.0 + t * t);
  109.             s = t * c;
  110.         }
  111.  
  112.         // Update elements of matrix A
  113.         for (int i = 0; i < N; i++) {
  114.             if (i != p && i != q) {
  115.                 double a_ip = A[i][p];
  116.                 double a_iq = A[i][q];
  117.                 A[i][p] = A[p][i] = c * a_ip + s * a_iq;
  118.                 A[i][q] = A[q][i] = -s * a_ip + c * a_iq;
  119.             }
  120.         }
  121.         A[p][p] = c * c * app + 2.0 * s * c * apq + s * s * aqq;
  122.         A[q][q] = s * s * app - 2.0 * s * c * apq + c * c * aqq;
  123.         A[p][q] = A[q][p] = 0.0; 
  124.  
  125.         // Update accumulated eigenvector matrix P
  126.         for (int i = 0; i < N; i++) {
  127.             double v_ip = eigenvectors[i][p];
  128.             double v_iq = eigenvectors[i][q];
  129.             eigenvectors[i][p] = c * v_ip + s * v_iq;
  130.             eigenvectors[i][q] = -s * v_ip + c * v_iq;
  131.         }
  132.     }
  133.  
  134.     // Extract eigenvalues from the diagonal elements
  135.     for (int i = 0; i < N; i++) {
  136.         eigenvalues[i] = A[i][i];
  137.     }
  138. }
  139.  
  140. // Verification function: Explicitly checks the defining equation: Ax - lambda * x = 0
  141. void verify_results(double orig_A[N][N], double eigenvalues[N], double eigenvectors[N][N]) {
  142.     printf("--- Verification Check (コード内での固有値・固有ベクトル関係性の確認) ---\n");
  143.     printf("Equation verified: Ax - lambda * x = 0\n\n");
  144.  
  145.     for (int k = 0; k < N; k++) {
  146.         double x[N];
  147.         double Ax[N];
  148.         double lambda_x[N];
  149.         double error_vector[N];
  150.  
  151.         // Extract the k-th eigenvector from the matrix columns
  152.         for (int i = 0; i < N; i++) {
  153.             x[i] = eigenvectors[i][k];
  154.         }
  155.  
  156.         // Calculate elements of vector Ax, vector lambda*x, and the resulting error vector
  157.         for (int i = 0; i < N; i++) {
  158.             Ax[i] = 0.0;
  159.             for (int j = 0; j < N; j++) {
  160.                 Ax[i] += orig_A[i][j] * x[j];
  161.             }
  162.             lambda_x[i] = eigenvalues[k] * x[i];
  163.             error_vector[i] = Ax[i] - lambda_x[i];
  164.         }
  165.  
  166.         // Print the element breakdown showing the equation results
  167.         printf("Pair %d Verification:\n", k + 1);
  168.         printf("  Ax vector         = [ %9.6f  %9.6f  %9.6f  %9.6f ]\n", Ax[0], Ax[1], Ax[2], Ax[3]);
  169.         printf("  lambda * x vector = [ %9.6f  %9.6f  %9.6f  %9.6f ]\n", lambda_x[0], lambda_x[1], lambda_x[2], lambda_x[3]);
  170.         printf("  Ax - lambda * x   = [ %9.6f  %9.6f  %9.6f  %9.6f ]\n\n", error_vector[0], error_vector[1], error_vector[2], error_vector[3]);
  171.     }
  172. }
  173.  
  174.  
Success #stdin #stdout 0s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
--- Iteration Process (収束状況) ---
1: Max off-diagonal element = 4.000000e+00
2: Max off-diagonal element = 2.000000e+00
3: Max off-diagonal element = 2.000000e+00
4: Max off-diagonal element = 0.000000e+00

[Convergence Achieved]

--- Final Results ---
Eigenvalue 1 (固有値): 10.000000
Eigenvector 1 (固有ベクトル): [  0.632456  0.632456  0.316228  0.316228 ]

Eigenvalue 2 (固有値):  1.000000
Eigenvector 2 (固有ベクトル): [ -0.707107  0.707107  0.000000  0.000000 ]

Eigenvalue 3 (固有値):  5.000000
Eigenvector 3 (固有ベクトル): [ -0.316228 -0.316228  0.632456  0.632456 ]

Eigenvalue 4 (固有値):  2.000000
Eigenvector 4 (固有ベクトル): [  0.000000  0.000000 -0.707107  0.707107 ]

--- Verification Check (コード内での固有値・固有ベクトル関係性の確認) ---
Equation verified: Ax - lambda * x = 0

Pair 1 Verification:
  Ax vector         = [  6.324555   6.324555   3.162278   3.162278 ]
  lambda * x vector = [  6.324555   6.324555   3.162278   3.162278 ]
  Ax - lambda * x   = [  0.000000   0.000000   0.000000   0.000000 ]

Pair 2 Verification:
  Ax vector         = [ -0.707107   0.707107   0.000000   0.000000 ]
  lambda * x vector = [ -0.707107   0.707107   0.000000   0.000000 ]
  Ax - lambda * x   = [  0.000000   0.000000   0.000000   0.000000 ]

Pair 3 Verification:
  Ax vector         = [ -1.581139  -1.581139   3.162278   3.162278 ]
  lambda * x vector = [ -1.581139  -1.581139   3.162278   3.162278 ]
  Ax - lambda * x   = [  0.000000   0.000000   0.000000   0.000000 ]

Pair 4 Verification:
  Ax vector         = [  0.000000   0.000000  -1.414214   1.414214 ]
  lambda * x vector = [  0.000000   0.000000  -1.414214   1.414214 ]
  Ax - lambda * x   = [  0.000000   0.000000  -0.000000   0.000000 ]
https://ideone.com/CHeAnv
language:
C (gcc 8.3)
created:
4 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!