fork(1) download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define N 4
  5. #define MAX_ITER 150
  6. #define EPS 1e-7
  7.  
  8. void print_system(double A[N][N], double b[N]);
  9. void jacobi_solver(double A[N][N], double b[N]);
  10. void gauss_seidel_solver(double A[N][N], double b[N]);
  11. void verify_solution(double A[N][N], double b[N], double x[N]);
  12.  
  13. int main() {
  14.     double A[N][N] = {
  15.         { 5.0, -1.0, -2.0,  0.0},
  16.         { 0.0,  8.0,  1.0,  7.0},
  17.         {-2.0, -7.0, -4.0,  0.0},
  18.         { 4.0,  0.0, -3.0, -25.0}
  19.     };
  20.     double b[N] = {-1.32, 14.41, 0.70, -0.03};
  21.  
  22.     print_system(A, b);
  23.  
  24.     printf("\n======================================================\n");
  25.     printf("               1. ヤコビの反復法\n");
  26.     printf("======================================================\n");
  27.     jacobi_solver(A, b);
  28.  
  29.     printf("\n======================================================\n");
  30.     printf("               2. ガウス・ザイデル法\n");
  31.     printf("======================================================\n");
  32.     gauss_seidel_solver(A, b);
  33.  
  34.     return 0;
  35. }
  36.  
  37. void print_system(double A[N][N], double b[N]) {
  38.     printf("対角優位化した連立一次方程式:\n");
  39.     for (int i = 0; i < N; i++) {
  40.         printf("[ ");
  41.         for (int j = 0; j < N; j++) {
  42.             printf("%6.1f ", A[i][j]);
  43.         }
  44.         printf("] [x%d] = [ %6.2f ]\n", i + 1, b[i]);
  45.     }
  46. }
  47.  
  48. void jacobi_solver(double A[N][N], double b[N]) {
  49.     double x[N] = {0.0, 0.0, 0.0, 0.0};
  50.     double x_new[N];
  51.     int converged = 0;
  52.  
  53.     printf("回数\t残差2乗和\t\tx1\t\tx2\t\tx3\t\tx4\n");
  54.     printf("------------------------------------------------------------------------\n");
  55.  
  56.     for (int k = 1; k <= MAX_ITER; k++) {
  57.         for (int i = 0; i < N; i++) {
  58.             double sum = 0.0;
  59.             for (int j = 0; j < N; j++) {
  60.                 if (j != i) {
  61.                     sum += A[i][j] * x[j];
  62.                 }
  63.             }
  64.             x_new[i] = (b[i] - sum) / A[i][i];
  65.         }
  66.  
  67.         double rss = 0.0;
  68.         for (int i = 0; i < N; i++) {
  69.             rss += pow(x_new[i] - x[i], 2);
  70.         }
  71.  
  72.         printf("%d  \t%.2e \t%.5f \t%.5f \t%.5f \t%.5f \n", k, rss, x_new[0], x_new[1], x_new[2], x_new[3]);
  73.  
  74.         for (int i = 0; i < N; i++) {
  75.             x[i] = x_new[i];
  76.         }
  77.  
  78.         if (rss < EPS) {
  79.             printf("\n>> 収束判定条件を満足しました（計算回数: %d回）\n", k);
  80.             converged = 1;
  81.             break;
  82.         }
  83.     }
  84.  
  85.     if (!converged) {
  86.         printf("\n警告: 指定された計算回数以内に収束しませんでした。\n");
  87.     }
  88.  
  89.     printf("\n=== 最終解 (ヤコビの反復法) ===\n");
  90.     for (int i = 0; i < N; i++) {
  91.         printf("x[%d] = %10.6f\n", i + 1, x[i]);
  92.     }
  93.  
  94.     verify_solution(A, b, x);
  95. }
  96.  
  97. void gauss_seidel_solver(double A[N][N], double b[N]) {
  98.     double x[N] = {0.0, 0.0, 0.0, 0.0};
  99.     double x_old[N];
  100.     int converged = 0;
  101.  
  102.     printf("回数\t\t残差2乗和\t\tx1\t\tx2\t\tx3\t\tx4\n");
  103.     printf("------------------------------------------------------------------------\n");
  104.  
  105.     for (int k = 1; k <= MAX_ITER; k++) {
  106.         for (int i = 0; i < N; i++) {
  107.             x_old[i] = x[i];
  108.         }
  109.  
  110.         for (int i = 0; i < N; i++) {
  111.             double sum = 0.0;
  112.             for (int j = 0; j < N; j++) {
  113.                 if (j != i) {
  114.                     sum += A[i][j] * x[j]; 
  115.                 }
  116.             }
  117.             x[i] = (b[i] - sum) / A[i][i];
  118.         }
  119.  
  120.         double rss = 0.0;
  121.         for (int i = 0; i < N; i++) {
  122.             rss += pow(x[i] - x_old[i], 2);
  123.         }
  124.  
  125.         printf("%d \t%.2e \t%.5f \t%.5f \t%.5f \t%.5f \n", k, rss, x[0], x[1], x[2], x[3]);
  126.  
  127.         if (rss < EPS) {
  128.             printf("\n>> 収束判定条件を満足しました（計算回数: %d回）\n", k);
  129.             converged = 1;
  130.             break;
  131.         }
  132.     }
  133.  
  134.     if (!converged) {
  135.         printf("\n警告: 指定された計算回数以内に収束しませんでした。\n");
  136.     }
  137.  
  138.     printf("\n=== 最終解 (ガウス・ザイデル法) ===\n");
  139.     for (int i = 0; i < N; i++) {
  140.         printf("x[%d] = %10.6f\n", i + 1, x[i]);
  141.     }
  142.  
  143.     verify_solution(A, b, x);
  144. }
  145.  
  146. void verify_solution(double A[N][N], double b[N], double x[N]) {
  147.     printf("\n--- コード内での解の確認（検算: Ax == b） ---\n");
  148.     for (int i = 0; i < N; i++) {
  149.         double lhs_calc = 0.0;
  150.         for (int j = 0; j < N; j++) {
  151.             lhs_calc += A[i][j] * x[j];
  152.         }
  153.         double error = lhs_calc - b[i];
  154.         printf("%d行目: 計算値 Ax = %10.5f | 設定値 b = %10.5f | 誤差 = %e\n", 
  155.                i + 1, lhs_calc, b[i], error);
  156.     }
  157. }
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
対角優位化した連立一次方程式:
[    5.0   -1.0   -2.0    0.0 ] [x1] = [  -1.32 ]
[    0.0    8.0    1.0    7.0 ] [x2] = [  14.41 ]
[   -2.0   -7.0   -4.0    0.0 ] [x3] = [   0.70 ]
[    4.0    0.0   -3.0  -25.0 ] [x4] = [  -0.03 ]

======================================================
               1. ヤコビの反復法
======================================================
回数	残差2乗和		x1		x2		x3		x4
------------------------------------------------------------------------
1  	3.34e+00 	-0.26400 	1.80125 	-0.17500 	0.00120 
2  	9.21e+00 	0.02625 	1.82208 	-3.19519 	-0.02004 
3  	1.81e+00 	-1.17766 	2.21818 	-3.37676 	0.38882 
4  	1.50e-01 	-1.17107 	1.88312 	-3.46799 	0.21799 
5  	3.77e-01 	-1.27457 	2.04401 	-2.88494 	0.22999 
6  	1.38e-01 	-1.00917 	1.96063 	-3.11474 	0.14346 
7  	2.78e-02 	-1.11777 	2.06506 	-3.10151 	0.21350 
8  	2.15e-02 	-1.09159 	2.00213 	-3.22998 	0.19454 
9  	1.50e-02 	-1.15557 	2.03478 	-3.13292 	0.21414 
10  	4.03e-03 	-1.11021 	2.00549 	-3.15808 	0.19226 
11  	1.67e-03 	-1.12613 	2.02778 	-3.12950 	0.20253 
12  	1.41e-03 	-1.11024 	2.01522 	-3.16055 	0.19656 
13  	5.42e-04 	-1.12518 	2.02433 	-3.14651 	0.20283 
14  	1.96e-04 	-1.11774 	2.01709 	-3.15499 	0.19875 
15  	1.30e-04 	-1.12258 	2.02171 	-3.14604 	0.20096 
16  	6.52e-05 	-1.11807 	2.01866 	-3.15171 	0.19911 
17  	2.52e-05 	-1.12095 	2.02099 	-3.14863 	0.20051 
18  	1.31e-05 	-1.11925 	2.01938 	-3.15126 	0.19968 
19  	7.24e-06 	-1.12063 	2.02043 	-3.14929 	0.20027 
20  	3.13e-06 	-1.11963 	2.01967 	-3.15045 	0.19981 
21  	1.46e-06 	-1.12024 	2.02022 	-3.14962 	0.20011 
22  	7.84e-07 	-1.11980 	2.01985 	-3.15026 	0.19991 
23  	3.72e-07 	-1.12013 	2.02011 	-3.14984 	0.20006 
24  	1.70e-07 	-1.11992 	2.01993 	-3.15012 	0.19996 
25  	8.60e-08 	-1.12006 	2.02005 	-3.14991 	0.20003 

>> 収束判定条件を満足しました（計算回数: 25回）

=== 最終解 (ヤコビの反復法) ===
x[1] =  -1.120063
x[2] =   2.020050
x[3] =  -3.149911
x[4] =   0.200028

--- コード内での解の確認（検算: Ax == b） ---
1行目: 計算値 Ax =   -1.32054 | 設定値 b =   -1.32000 | 誤差 = -5.434164e-04
2行目: 計算値 Ax =   14.41069 | 設定値 b =   14.41000 | 誤差 = 6.884127e-04
3行目: 計算値 Ax =    0.69942 | 設定値 b =    0.70000 | 誤差 = -5.808020e-04
4行目: 計算値 Ax =   -0.03122 | 設定値 b =   -0.03000 | 誤差 = -1.219528e-03

======================================================
               2. ガウス・ザイデル法
======================================================
回数		残差2乗和		x1		x2		x3		x4
------------------------------------------------------------------------
1 	1.36e+01 	-0.26400 	1.80125 	-3.19519 	0.34238 
2 	9.66e-01 	-1.18182 	1.90106 	-2.91095 	0.16142 
3 	1.15e-01 	-1.04817 	2.02387 	-3.19270 	0.21662 
4 	1.29e-02 	-1.13630 	2.01080 	-3.12574 	0.19448 
5 	1.74e-03 	-1.11214 	2.02180 	-3.15708 	0.20211 
6 	2.22e-04 	-1.12247 	2.01904 	-3.14709 	0.19925 
7 	2.97e-05 	-1.11903 	2.02029 	-3.15099 	0.20027 
8 	3.88e-06 	-1.12034 	2.01988 	-3.14963 	0.19990 
9 	5.15e-07 	-1.11987 	2.02004 	-3.15013 	0.20004 
10 	6.77e-08 	-1.12005 	2.01999 	-3.14995 	0.19999 

>> 収束判定条件を満足しました（計算回数: 10回）

=== 最終解 (ガウス・ザイデル法) ===
x[1] =  -1.120045
x[2] =   2.019985
x[3] =  -3.149951
x[4] =   0.199987

--- コード内での解の確認（検算: Ax == b） ---
1行目: 計算値 Ax =   -1.32031 | 設定値 b =   -1.32000 | 誤差 = -3.084483e-04
2行目: 計算値 Ax =   14.40984 | 設定値 b =   14.41000 | 誤差 = -1.626182e-04
3行目: 計算値 Ax =    0.70000 | 設定値 b =    0.70000 | 誤差 = -6.661338e-16
4行目: 計算値 Ax =   -0.03000 | 設定値 b =   -0.03000 | 誤差 = -2.498002e-16
https://ideone.com/oYkmVi
language:
C (gcc 8.3)
created:
17 hours 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!