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 detailed iteration step 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 step = 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.         // Check for convergence
  80.         if (max_off_diag < EPSILON) {
  81.             printf("\n[Convergence Achieved] All remaining off-diagonal elements are below the threshold.\n");
  82.             break;
  83.         }
  84.  
  85.         step++;
  86.         // Print detailed iteration step showing exactly which elements are shrinking
  87.         printf("Step %2d | Eliminating A[%d][%d] = %9.6f | Current Diagonals: [%.4f, %.4f, %.4f, %.4f]\n", 
  88.                step, p, q, A[p][q], A[0][0], A[1][1], A[2][2], A[3][3]);
  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 outputs vector equation components and tests vector validity
  141. void verify_results(double orig_A[N][N], double eigenvalues[N], double eigenvectors[N][N]) {
  142.     printf("--- Verification Check (Ax = lambda * x & Vector Validity) ---\n");
  143.  
  144.     for (int k = 0; k < N; k++) {
  145.         double x[N];
  146.         double Ax[N];
  147.         double lambda_x[N];
  148.         double vector_norm = 0.0;
  149.  
  150.         // Extract the k-th eigenvector
  151.         for (int i = 0; i < N; i++) {
  152.             x[i] = eigenvectors[i][k];
  153.             vector_norm += x[i] * x[i]; // sum of squares
  154.         }
  155.         vector_norm = sqrt(vector_norm); // magnitude of the vector
  156.  
  157.         // Calculate elements of vector Ax and vector lambda*x
  158.         for (int i = 0; i < N; i++) {
  159.             Ax[i] = 0.0;
  160.             for (int j = 0; j < N; j++) {
  161.                 Ax[i] += orig_A[i][j] * x[j];
  162.             }
  163.             lambda_x[i] = eigenvalues[k] * x[i];
  164.         }
  165.  
  166.         // Print results
  167.         printf("Pair %d Verification:\n", k + 1);
  168.         printf("  Eigenvector Norm (Vector Length) = %9.6f (Should be 1.0)\n", vector_norm);
  169.         printf("  LHS Vector (Ax)                  = [ %9.6f  %9.6f  %9.6f  %9.6f ]\n", Ax[0], Ax[1], Ax[2], Ax[3]);
  170.         printf("  RHS Vector (lambda * x)          = [ %9.6f  %9.6f  %9.6f  %9.6f ]\n\n", lambda_x[0], lambda_x[1], lambda_x[2], lambda_x[3]);
  171.     }
  172. }
  173.  
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
--- Iteration Process (収束状況) ---
Step  1 | Eliminating A[0][1] =  4.000000 | Current Diagonals: [5.0000, 5.0000, 4.0000, 4.0000]
Step  2 | Eliminating A[2][3] =  2.000000 | Current Diagonals: [9.0000, 1.0000, 4.0000, 4.0000]
Step  3 | Eliminating A[0][2] =  2.000000 | Current Diagonals: [9.0000, 1.0000, 6.0000, 2.0000]

[Convergence Achieved] All remaining off-diagonal elements are below the threshold.

--- 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 (Ax = lambda * x & Vector Validity) ---
Pair 1 Verification:
  Eigenvector Norm (Vector Length) =  1.000000 (Should be 1.0)
  LHS Vector (Ax)                  = [  6.324555   6.324555   3.162278   3.162278 ]
  RHS Vector (lambda * x)          = [  6.324555   6.324555   3.162278   3.162278 ]

Pair 2 Verification:
  Eigenvector Norm (Vector Length) =  1.000000 (Should be 1.0)
  LHS Vector (Ax)                  = [ -0.707107   0.707107   0.000000   0.000000 ]
  RHS Vector (lambda * x)          = [ -0.707107   0.707107   0.000000   0.000000 ]

Pair 3 Verification:
  Eigenvector Norm (Vector Length) =  1.000000 (Should be 1.0)
  LHS Vector (Ax)                  = [ -1.581139  -1.581139   3.162278   3.162278 ]
  RHS Vector (lambda * x)          = [ -1.581139  -1.581139   3.162278   3.162278 ]

Pair 4 Verification:
  Eigenvector Norm (Vector Length) =  1.000000 (Should be 1.0)
  LHS Vector (Ax)                  = [  0.000000   0.000000  -1.414214   1.414214 ]
  RHS Vector (lambda * x)          = [  0.000000   0.000000  -1.414214   1.414214 ]
https://ideone.com/UnmZ5R
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!