fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3. #define N 4
  4.  
  5. void beki(double l[][N], double u[][N], double x[N], double lam, double lamh);
  6. void con(double a[][N], double x[N]);
  7. void pivot(double A[][N], int i);
  8. void lu(double A[][N], double l[][N], double u[][N]);
  9. void Alam(double A[][N], double lamh);
  10.  
  11. void Alam(double A[][N], double lamh){
  12.     for(int i = 0; i < N; i++){
  13.         A[i][i] -= lamh; 
  14.     }
  15. }
  16.  
  17. void pivot(double A[][N], int i) {
  18.     double c[N] = {0};
  19.     for (int j = i; j < N; j++) {
  20.         if (fabs(A[i][i]) < fabs(A[j][i])) {
  21.             for (int n = 0; n < N; n++) {
  22.                 c[n] = A[j][n];
  23.                 A[j][n] = A[i][n];
  24.                 A[i][n] = c[n];
  25.             }
  26.         }
  27.     }
  28. }
  29.  
  30. void lu(double A[][N], double l[][N], double u[][N]) {
  31.     for (int i = 0; i < N; i++) {
  32.         for(int n=i;n<N;n++){
  33.         	l[n][i] = A[n][i]; 
  34.         }
  35.         double det = A[i][i];
  36.         for(int m=0;m<N;m++){
  37.         	A[i][m]/=det;
  38.         }
  39.         for (int a = i + 1; a < N; a++) {
  40.         	double det = A[a][i];
  41.             for (int b = i; b < N; b++) {
  42.                 A[a][b] -= A[i][b]*det;
  43.             }
  44.         }
  45.         for (int j = i; j < N; j++) {
  46.             u[i][j] = A[i][j];
  47.         }
  48.     }
  49. }
  50.  
  51. void beki(double l[][N], double u[][N], double x[N], double lam, double lamh) {
  52.     int n = 1;
  53.     while(1){
  54.         double yi[N] = {0};
  55.         double xi[N] = {0};
  56.  
  57.         // yi
  58.         for(int i = 0; i < N; i++){
  59.             yi[i] = x[i];
  60.             for(int j = 0; j < N; j++){
  61.                 if(j != i){
  62.                     yi[i] -= l[i][j] * yi[j];
  63.                 }
  64.             }
  65.             yi[i] /= l[i][i];
  66.         }
  67.  
  68.         // xi
  69.         for(int i = N - 1; i >= 0; i--){
  70.             xi[i] = yi[i];
  71.             for(int j = 0; j < N; j++){
  72.                 if(j != i){
  73.                     xi[i] -= u[i][j] * xi[j];
  74.                 }
  75.             }
  76.             xi[i] /= u[i][i];
  77.         }
  78.  
  79.         lam = xi[0];
  80.         for(int i = 0; i < N; i++){
  81.             xi[i] /= lam;
  82.         }
  83.  
  84.         printf("1/(λ-λh) = %f\n", lam);
  85.         printf("x%d = ", n);
  86.         for(int i = 0; i < N; i++){
  87.             printf("%f, ", xi[i]);
  88.         }
  89.         printf("\n");
  90.         printf("\n");
  91.  
  92.         if(fabs(x[1] - xi[1]) < 1e-6 && fabs(x[2] - xi[2]) < 1e-6 && fabs(x[3] - xi[3]) < 1e-6) {
  93.             for(int i = 0; i < N; i++){
  94.                 x[i] = xi[i];
  95.             }
  96.             lam = lamh + 1 / lam;
  97.             printf("λ = %f\n", lam);
  98.             printf("\n");
  99.             break;
  100.         }
  101.  
  102.         for(int i = 0; i < N; i++){
  103.             x[i] = xi[i];
  104.         }
  105.         n++;
  106.     }
  107. }
  108.  
  109. void con(double a[][N], double x[N]) {
  110.     double xi[N] = {0};
  111.     double l = 0;
  112.  
  113.     for(int i = 0; i < N; i++){
  114.         for(int j = 0; j < N; j++){
  115.             xi[i] += a[i][j] * x[j];
  116.         }
  117.     }
  118.  
  119.     l = xi[0];
  120.         for(int i = 0; i < N; i++){
  121.         xi[i] /= l;
  122.     }
  123.  
  124.     printf("固有値・固有ベクトル確認\n");
  125.     printf("λ = %f\n", l);
  126.     printf("x = ");
  127.     for(int i = 0; i < N; i++){
  128.         printf("%f, ", xi[i]);
  129.     }
  130.     printf("\n");
  131. }
  132.  
  133. int main(void) {
  134.     double A[N][N] = {
  135.         {2, 1, 5, 3},
  136.         {3, 0, 1, 6},
  137.         {1, 4, 3, 3},
  138.         {8, 2, 0, 1}
  139.     };
  140.     double x[N] = {1, 1, 1, 1};
  141.     double lamh = -4.678+100;
  142.     double lam = 0;
  143.  
  144.     double l[N][N] = {0};
  145.     double u[N][N] = {0};
  146.  
  147.     double a[N][N];
  148.     for (int i = 0; i < N; i++) {
  149.         for (int j = 0; j < N; j++) {
  150.             a[i][j] = A[i][j];
  151.         }
  152.     }
  153.  
  154.     Alam(A, lamh);  // λhを反映させる
  155.     lu(A, l, u);  // LU分解
  156.     beki(l, u, x, lam, lamh);  // 逆反復法
  157.     con(a, x);  // 固有値・固有ベクトル確認
  158.  
  159.     return 0;
  160. }
Success #stdin #stdout 0.01s 5276KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1/(λ-λh) = -0.011858
x1 = 1.000000, 0.989629, 0.999680, 0.999910, 

1/(λ-λh) = -0.011856
x2 = 1.000000, 0.980565, 0.999124, 0.999767, 

1/(λ-λh) = -0.011854
x3 = 1.000000, 0.972651, 0.998398, 0.999600, 

1/(λ-λh) = -0.011853
x4 = 1.000000, 0.965747, 0.997555, 0.999426, 

1/(λ-λh) = -0.011851
x5 = 1.000000, 0.959728, 0.996639, 0.999262, 

1/(λ-λh) = -0.011849
x6 = 1.000000, 0.954486, 0.995682, 0.999114, 

1/(λ-λh) = -0.011848
x7 = 1.000000, 0.949925, 0.994712, 0.998989, 

1/(λ-λh) = -0.011847
x8 = 1.000000, 0.945960, 0.993748, 0.998889, 

1/(λ-λh) = -0.011846
x9 = 1.000000, 0.942518, 0.992805, 0.998815, 

1/(λ-λh) = -0.011845
x10 = 1.000000, 0.939533, 0.991894, 0.998767, 

1/(λ-λh) = -0.011844
x11 = 1.000000, 0.936948, 0.991024, 0.998744, 

1/(λ-λh) = -0.011843
x12 = 1.000000, 0.934712, 0.990200, 0.998743, 

1/(λ-λh) = -0.011842
x13 = 1.000000, 0.932780, 0.989424, 0.998762, 

1/(λ-λh) = -0.011841
x14 = 1.000000, 0.931115, 0.988699, 0.998798, 

1/(λ-λh) = -0.011841
x15 = 1.000000, 0.929681, 0.988026, 0.998850, 

1/(λ-λh) = -0.011840
x16 = 1.000000, 0.928448, 0.987403, 0.998914, 

1/(λ-λh) = -0.011839
x17 = 1.000000, 0.927391, 0.986829, 0.998988, 

1/(λ-λh) = -0.011839
x18 = 1.000000, 0.926485, 0.986302, 0.999070, 

1/(λ-λh) = -0.011839
x19 = 1.000000, 0.925712, 0.985821, 0.999158, 

1/(λ-λh) = -0.011838
x20 = 1.000000, 0.925053, 0.985383, 0.999250, 

1/(λ-λh) = -0.011838
x21 = 1.000000, 0.924492, 0.984985, 0.999345, 

1/(λ-λh) = -0.011838
x22 = 1.000000, 0.924017, 0.984625, 0.999440, 

1/(λ-λh) = -0.011837
x23 = 1.000000, 0.923616, 0.984299, 0.999536, 

1/(λ-λh) = -0.011837
x24 = 1.000000, 0.923278, 0.984006, 0.999630, 

1/(λ-λh) = -0.011837
x25 = 1.000000, 0.922994, 0.983743, 0.999723, 

1/(λ-λh) = -0.011837
x26 = 1.000000, 0.922757, 0.983506, 0.999813, 

1/(λ-λh) = -0.011837
x27 = 1.000000, 0.922559, 0.983295, 0.999901, 

1/(λ-λh) = -0.011837
x28 = 1.000000, 0.922396, 0.983106, 0.999984, 

1/(λ-λh) = -0.011836
x29 = 1.000000, 0.922261, 0.982938, 1.000064, 

1/(λ-λh) = -0.011836
x30 = 1.000000, 0.922152, 0.982789, 1.000140, 

1/(λ-λh) = -0.011836
x31 = 1.000000, 0.922062, 0.982657, 1.000212, 

1/(λ-λh) = -0.011836
x32 = 1.000000, 0.921991, 0.982539, 1.000279, 

1/(λ-λh) = -0.011836
x33 = 1.000000, 0.921934, 0.982436, 1.000343, 

1/(λ-λh) = -0.011836
x34 = 1.000000, 0.921890, 0.982344, 1.000402, 

1/(λ-λh) = -0.011836
x35 = 1.000000, 0.921856, 0.982264, 1.000457, 

1/(λ-λh) = -0.011836
x36 = 1.000000, 0.921830, 0.982193, 1.000509, 

1/(λ-λh) = -0.011836
x37 = 1.000000, 0.921812, 0.982131, 1.000556, 

1/(λ-λh) = -0.011836
x38 = 1.000000, 0.921799, 0.982077, 1.000600, 

1/(λ-λh) = -0.011836
x39 = 1.000000, 0.921791, 0.982030, 1.000641, 

1/(λ-λh) = -0.011836
x40 = 1.000000, 0.921787, 0.981988, 1.000678, 

1/(λ-λh) = -0.011836
x41 = 1.000000, 0.921785, 0.981953, 1.000712, 

1/(λ-λh) = -0.011836
x42 = 1.000000, 0.921786, 0.981921, 1.000743, 

1/(λ-λh) = -0.011836
x43 = 1.000000, 0.921789, 0.981895, 1.000772, 

1/(λ-λh) = -0.011836
x44 = 1.000000, 0.921794, 0.981871, 1.000798, 

1/(λ-λh) = -0.011836
x45 = 1.000000, 0.921799, 0.981851, 1.000821, 

1/(λ-λh) = -0.011836
x46 = 1.000000, 0.921805, 0.981834, 1.000843, 

1/(λ-λh) = -0.011836
x47 = 1.000000, 0.921812, 0.981819, 1.000862, 

1/(λ-λh) = -0.011836
x48 = 1.000000, 0.921818, 0.981807, 1.000880, 

1/(λ-λh) = -0.011836
x49 = 1.000000, 0.921825, 0.981796, 1.000896, 

1/(λ-λh) = -0.011836
x50 = 1.000000, 0.921832, 0.981787, 1.000910, 

1/(λ-λh) = -0.011836
x51 = 1.000000, 0.921839, 0.981779, 1.000923, 

1/(λ-λh) = -0.011836
x52 = 1.000000, 0.921846, 0.981773, 1.000935, 

1/(λ-λh) = -0.011836
x53 = 1.000000, 0.921852, 0.981767, 1.000945, 

1/(λ-λh) = -0.011836
x54 = 1.000000, 0.921858, 0.981763, 1.000955, 

1/(λ-λh) = -0.011836
x55 = 1.000000, 0.921864, 0.981759, 1.000963, 

1/(λ-λh) = -0.011836
x56 = 1.000000, 0.921869, 0.981756, 1.000971, 

1/(λ-λh) = -0.011836
x57 = 1.000000, 0.921874, 0.981754, 1.000977, 

1/(λ-λh) = -0.011836
x58 = 1.000000, 0.921879, 0.981752, 1.000983, 

1/(λ-λh) = -0.011836
x59 = 1.000000, 0.921883, 0.981750, 1.000989, 

1/(λ-λh) = -0.011836
x60 = 1.000000, 0.921887, 0.981749, 1.000993, 

1/(λ-λh) = -0.011836
x61 = 1.000000, 0.921891, 0.981748, 1.000997, 

1/(λ-λh) = -0.011836
x62 = 1.000000, 0.921894, 0.981747, 1.001001, 

1/(λ-λh) = -0.011836
x63 = 1.000000, 0.921897, 0.981746, 1.001004, 

1/(λ-λh) = -0.011836
x64 = 1.000000, 0.921900, 0.981746, 1.001007, 

1/(λ-λh) = -0.011836
x65 = 1.000000, 0.921903, 0.981746, 1.001010, 

1/(λ-λh) = -0.011836
x66 = 1.000000, 0.921905, 0.981746, 1.001012, 

1/(λ-λh) = -0.011836
x67 = 1.000000, 0.921907, 0.981746, 1.001014, 

1/(λ-λh) = -0.011836
x68 = 1.000000, 0.921909, 0.981746, 1.001016, 

1/(λ-λh) = -0.011836
x69 = 1.000000, 0.921911, 0.981746, 1.001017, 

1/(λ-λh) = -0.011836
x70 = 1.000000, 0.921913, 0.981746, 1.001019, 

1/(λ-λh) = -0.011836
x71 = 1.000000, 0.921914, 0.981746, 1.001020, 

1/(λ-λh) = -0.011836
x72 = 1.000000, 0.921915, 0.981746, 1.001021, 

1/(λ-λh) = -0.011836
x73 = 1.000000, 0.921917, 0.981746, 1.001022, 

1/(λ-λh) = -0.011836
x74 = 1.000000, 0.921918, 0.981746, 1.001022, 

1/(λ-λh) = -0.011836
x75 = 1.000000, 0.921919, 0.981747, 1.001023, 

λ = 10.833720

固有値・固有ベクトル確認
λ = 10.833720
x = 1.000000, 0.921926, 0.981748, 1.001028,
https://ideone.com/UMRONw
language:
C (gcc 8.3)
created:
22 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!