fork(1) download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. #define SIZE 4
  6. #define EPSILON_POWER 1e-9  
  7. #define EPSILON_NEWTON 1e-7  
  8. #define MAX_ITER 1000  
  9.  
  10.  
  11. double vector_norm(double *v, int n) {
  12.     double sum = 0.0;
  13.     for (int i = 0; i < n; i++) {
  14.         sum += v[i] * v[i];
  15.     }
  16.     return sqrt(sum);
  17. }
  18.  
  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.  
  29. double f(double lambda) {
  30.     return pow(lambda, 4) - 2.0 * pow(lambda, 3) - 43.0 * pow(lambda, 2) - 36.0 * lambda + 130.0;
  31. }
  32.  
  33. double df(double lambda) {
  34.     return 4.0 * pow(lambda, 3) - 6.0 * pow(lambda, 2) - 86.0 * lambda - 36.0;
  35. }
  36.  
  37. void solve_newton(int attempt, double lambda0) {
  38.     double lambda = lambda0;
  39.     int iter = 0;
  40.  
  41.     while (iter <= MAX_ITER) {
  42.         double f_val = f(lambda);
  43.         double df_val = df(lambda);
  44.  
  45.         if (fabs(df_val) < 1e-12) return; 
  46.  
  47.         double delta = f_val / df_val;
  48.         lambda -= delta;
  49.         iter++;
  50.  
  51.         if (fabs(delta) < EPSILON_NEWTON) {
  52.             printf("試行 #%d (初期値 %5.1f) -> 収束成功 (%2d回) | 固有値解 λ = %.10f\n", 
  53.                    attempt, lambda0, iter, lambda);
  54.             return;
  55.         }
  56.     }
  57. }
  58.  
  59. // --- メイン関数 ---
  60. int main() {
  61.     // 2つ目の行列 A
  62.     double A[SIZE][SIZE] = {
  63.         {2.0,  1.0, 2.0,  5.0},
  64.         {3.0, -3.0, 1.0,  6.0},
  65.         {0.0,  1.0, 2.0,  7.0},
  66.         {1.0,  1.0, 3.0,  1.0}
  67.     };
  68.  
  69.     double x[SIZE], y[SIZE], Ax[SIZE];
  70.     double lambda_old = 0.0, lambda_new = 0.0;
  71.  
  72.     // 1. べき乗法を最初に実行
  73.     printf("==================================================\n");
  74.     printf("              【 手法 1 : べき乗法 】             \n");
  75.     printf("==================================================\n");
  76.  
  77.     for (int i = 0; i < SIZE; i++) {
  78.         x[i] = 1.0;
  79.     }
  80.  
  81.     printf("初期ベクトル x(0): [%.1f, %.1f, %.1f, %.1f]\n\n", x[0], x[1], x[2], x[3]);
  82.     printf("%-5s | %-12s | %-12s | %-30s\n", "Iter", "固有値(新)", "差分値", "固有ベクトルx(k)");
  83.  
  84.     int actual_iters = 0;
  85.     for (int iter = 1; iter <= MAX_ITER; iter++) {
  86.         actual_iters = iter;
  87.         mat_vec_mul(A, x, y, SIZE);
  88.  
  89.         double norm_y = vector_norm(y, SIZE);
  90.         if (norm_y < 1e-12) break;
  91.         for (int i = 0; i < SIZE; i++) {
  92.             y[i] /= norm_y;
  93.         }
  94.  
  95.         mat_vec_mul(A, y, Ax, SIZE);
  96.         double num = 0.0, den = 0.0;
  97.         for (int i = 0; i < SIZE; i++) {
  98.             num += y[i] * Ax[i];
  99.             den += y[i] * y[i];
  100.         }
  101.         lambda_new = num / den;
  102.         double diff = fabs(lambda_new - lambda_old);
  103.  
  104.         if (iter <= 5 || iter % 5 == 0) {
  105.             printf("%-5d | %-12.8f | %-12.4e | [%.5f, %.5f, %.5f, %.5f]\n", 
  106.                    iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
  107.         }
  108.  
  109.         if (diff < EPSILON_POWER) {
  110.             if (iter % 5 != 0 && iter > 5) {
  111.                 printf("%-5d | %-12.8f | %-12.4e | [%.5f, %.5f, %.5f, %.5f]\n", 
  112.                        iter, lambda_new, diff, y[0], y[1], y[2], y[3]);
  113.             }
  114.             break;
  115.         }
  116.  
  117.         for (int i = 0; i < SIZE; i++) {
  118.             x[i] = y[i];
  119.         }
  120.         lambda_old = lambda_new;
  121.     }
  122.  
  123.     printf("反復回数 %d 回で収束しました。\n", actual_iters);
  124.     printf("\n最大固有値 (λ1)          : %.10f\n", lambda_new);
  125.     printf("対応する固有ベクトル (x1) : [");
  126.     for (int i = 0; i < SIZE; i++) {
  127.         printf("%.10f%s", y[i], (i == SIZE - 1) ? "]\n" : ", ");
  128.     }
  129.  
  130.     printf("\n   検証: Ax と λx の比較   \n");
  131.     mat_vec_mul(A, y, Ax, SIZE);
  132.     printf("Ax    : [");
  133.     for (int i = 0; i < SIZE; i++) {
  134.         printf("%.6f%s", Ax[i], (i == SIZE - 1) ? "]\n" : ", ");
  135.     }
  136.     printf("λ * x : [");
  137.     for (int i = 0; i < SIZE; i++) {
  138.         printf("%.6f%s", lambda_new * y[i], (i == SIZE - 1) ? "]\n" : ", ");
  139.     }
  140.  
  141.  
  142.     // 【指示通りに2行空ける】
  143.     printf("\n\n");
  144.  
  145.  
  146.     // 2. ニュートン・ラフソン法を実行
  147.     printf("==================================================\n");
  148.     printf("        【 手法 2 : ニュートン・ラフソン法 】     \n");
  149.     printf("==================================================\n");
  150.     printf("対象の特性方程式: f(λ) = λ^4 - 2.0λ^3 - 43.0λ^2 - 36.0λ + 130.0 = 0\n\n");
  151.  
  152.     double initial_lambda[] = {10.0, 2.0, -2.5, -6.0};
  153.     for (int i = 0; i < 4; i++) {
  154.         solve_newton(i + 1, initial_lambda[i]);
  155.     }
  156.     printf("==================================================\n");
  157.  
  158.     return 0;
  159. }
Success #stdin #stdout 0.01s 5316KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
==================================================
              【 手法 1 : べき乗法 】             
==================================================
初期ベクトル x(0): [1.0, 1.0, 1.0, 1.0]

Iter  | 固有値(新) | 差分値    | 固有ベクトルx(k)        
1     | 7.58947368   | 7.5895e+00   | [0.59235, 0.41464, 0.59235, 0.35541]
2     | 7.97210119   | 3.8263e-01   | [0.59901, 0.42786, 0.53677, 0.41230]
3     | 7.72878414   | 2.4332e-01   | [0.59682, 0.44176, 0.55001, 0.38228]
4     | 7.86852048   | 1.3974e-01   | [0.60110, 0.42802, 0.54559, 0.39724]
5     | 7.78782087   | 8.0700e-02   | [0.59823, 0.43820, 0.54642, 0.38926]
10    | 7.81916110   | 4.7095e-03   | [0.59949, 0.43391, 0.54638, 0.39217]
15    | 7.81737213   | 2.6407e-04   | [0.59941, 0.43420, 0.54635, 0.39202]
20    | 7.81747243   | 1.4759e-05   | [0.59941, 0.43418, 0.54635, 0.39203]
25    | 7.81746682   | 8.2460e-07   | [0.59941, 0.43418, 0.54635, 0.39203]
30    | 7.81746714   | 4.6068e-08   | [0.59941, 0.43418, 0.54635, 0.39203]
35    | 7.81746712   | 2.5737e-09   | [0.59941, 0.43418, 0.54635, 0.39203]
37    | 7.81746712   | 8.1175e-10   | [0.59941, 0.43418, 0.54635, 0.39203]
反復回数 37 回で収束しました。

最大固有値 (λ1)          : 7.8174671191
対応する固有ベクトル (x1) : [0.5994100822, 0.4341839804, 0.5463547658, 0.3920309865]

   検証: Ax と λx の比較   
Ax    : [4.685869, 3.394219, 4.271110, 3.064689]
λ * x : [4.685869, 3.394219, 4.271110, 3.064689]


==================================================
        【 手法 2 : ニュートン・ラフソン法 】     
==================================================
対象の特性方程式: f(λ) = λ^4 - 2.0λ^3 - 43.0λ^2 - 36.0λ + 130.0 = 0

試行 #1 (初期値  10.0) -> 収束成功 ( 6回) | 固有値解 λ = 7.8174671194
試行 #2 (初期値   2.0) -> 収束成功 ( 5回) | 固有値解 λ = 1.3593330358
試行 #3 (初期値  -2.5) -> 収束成功 ( 4回) | 固有値解 λ = -2.7864676299
試行 #4 (初期値  -6.0) -> 収束成功 ( 7回) | 固有値解 λ = -4.3903325253
==================================================
https://ideone.com/wsXzOd
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!