fork 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.         {10.0,  2.0, -3.0,  1.0},
  16.         { 1.0,  9.0,  2.0, -1.0},
  17.         { 2.0, -1.0,  8.0,  3.0},
  18.         { 3.0,  2.0,  1.0,  7.0}
  19.     };
  20.     double b[N] = {13.0, -5.0, 41.0, 35.0};
  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残差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 0.01s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
対角優位化後の方程式:
[   10.0    2.0   -3.0    1.0 ] [x1] = [  13.00 ]
[    1.0    9.0    2.0   -1.0 ] [x2] = [  -5.00 ]
[    2.0   -1.0    8.0    3.0 ] [x3] = [  41.00 ]
[    3.0    2.0    1.0    7.0 ] [x4] = [  35.00 ]

======================================================
 1. ヤコビの反復法
======================================================
回数	残差2乗和		x1		x2		x3		x4
------------------------------------------------------------------------
1  	5.33e+01 	1.30000 	-0.55556 	5.12500 	5.00000 
2  	8.28e+00 	2.44861 	-1.28333 	2.85556 	3.86944 
3  	2.45e-01 	2.02639 	-1.03225 	2.90139 	3.90933 
4  	2.88e-02 	1.98593 	-0.99109 	3.02337 	4.01199 
5  	1.13e-03 	2.00403 	-1.00230 	3.00013 	4.00014 
6  	1.97e-05 	2.00049 	-1.00046 	2.99865 	3.99891 
7  	4.77e-06 	1.99980 	-0.99988 	3.00023 	4.00012 
8  	1.28e-07 	2.00003 	-1.00002 	3.00002 	4.00002 
9  	3.23e-09 	2.00001 	-1.00001 	2.99998 	3.99999 

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

=== 連立方程式の最終解 (ヤコビの反復法) ===
x[1] =   2.000008
x[2] =  -1.000007
x[3] =   2.999983
x[4] =   3.999987

--- 求まったものが解である確認（コード内での検算: Ax == b） ---
1行目: 計算値 Ax =   13.00011 | 設定値 b =   13.00000 | 誤差 = 1.057448e-04
2行目: 計算値 Ax =   -5.00007 | 設定値 b =   -5.00000 | 誤差 = -7.238616e-05
3行目: 計算値 Ax =   40.99985 | 設定値 b =   41.00000 | 誤差 = -1.516174e-04
4行目: 計算値 Ax =   34.99991 | 設定値 b =   35.00000 | 誤差 = -9.467652e-05

======================================================
 2. ガウス・ザイデル法
======================================================
回数	残差2乗和		x1		x2		x3		x4
------------------------------------------------------------------------
1  	4.01e+01 	1.30000 	-0.70000 	4.71250 	3.96964 
2  	5.37e+00 	2.45679 	-1.43468 	2.84285 	3.95088 
3  	4.09e-01 	2.04470 	-0.97550 	3.01031 	3.97237 
4  	3.57e-03 	2.00096 	-1.00547 	3.00944 	3.99980 
5  	1.32e-04 	2.00394 	-1.00256 	2.99877 	3.99922 
6  	2.39e-05 	2.00022 	-0.99984 	3.00026 	3.99982 
7  	1.62e-07 	2.00006 	-1.00008 	3.00004 	3.99999 
8  	8.29e-09 	2.00003 	-1.00001 	2.99999 	3.99999 

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

=== 連立方程式の最終解 (ガウス・ザイデル法) ===
x[1] =   2.000030
x[2] =  -1.000013
x[3] =   2.999994
x[4] =   3.999992

--- 求まったものが解である確認（コード内での検算: Ax == b） ---
1行目: 計算値 Ax =   13.00028 | 設定値 b =   13.00000 | 誤差 = 2.813554e-04
2行目: 計算値 Ax =   -5.00009 | 設定値 b =   -5.00000 | 誤差 = -9.294962e-05
3行目: 計算値 Ax =   41.00000 | 設定値 b =   41.00000 | 誤差 = 1.570910e-06
4行目: 計算値 Ax =   35.00000 | 設定値 b =   35.00000 | 誤差 = 0.000000e+00
https://ideone.com/jGlsZX
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!