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\t%.2e\t%.4f\t%.4f\t\t%.4ft%.4f\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\t%.4f\t%.4f\t%.4f\t%.4f\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.5f\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. }
  158.  
Success #stdin #stdout 0s 5324KB
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.2640	1.8013		-0.1750t0.0012
2		9.21e+00	0.0262	1.8221		-3.1952t-0.0200
3		1.81e+00	-1.1777	2.2182		-3.3768t0.3888
4		1.50e-01	-1.1711	1.8831		-3.4680t0.2180
5		3.77e-01	-1.2746	2.0440		-2.8849t0.2300
6		1.38e-01	-1.0092	1.9606		-3.1147t0.1435
7		2.78e-02	-1.1178	2.0651		-3.1015t0.2135
8		2.15e-02	-1.0916	2.0021		-3.2300t0.1945
9		1.50e-02	-1.1556	2.0348		-3.1329t0.2141
10		4.03e-03	-1.1102	2.0055		-3.1581t0.1923
11		1.67e-03	-1.1261	2.0278		-3.1295t0.2025
12		1.41e-03	-1.1102	2.0152		-3.1606t0.1966
13		5.42e-04	-1.1252	2.0243		-3.1465t0.2028
14		1.96e-04	-1.1177	2.0171		-3.1550t0.1988
15		1.30e-04	-1.1226	2.0217		-3.1460t0.2010
16		6.52e-05	-1.1181	2.0187		-3.1517t0.1991
17		2.52e-05	-1.1210	2.0210		-3.1486t0.2005
18		1.31e-05	-1.1193	2.0194		-3.1513t0.1997
19		7.24e-06	-1.1206	2.0204		-3.1493t0.2003
20		3.13e-06	-1.1196	2.0197		-3.1504t0.1998
21		1.46e-06	-1.1202	2.0202		-3.1496t0.2001
22		7.84e-07	-1.1198	2.0199		-3.1503t0.1999
23		3.72e-07	-1.1201	2.0201		-3.1498t0.2001
24		1.70e-07	-1.1199	2.0199		-3.1501t0.2000
25		8.60e-08	-1.1201	2.0201		-3.1499t0.2000

>> 収束判定条件を満足しました（計算回数: 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.2640	1.8013	-3.1952	0.3424
2	9.66e-01		-1.1818	1.9011	-2.9109	0.1614
3	1.15e-01		-1.0482	2.0239	-3.1927	0.2166
4	1.29e-02		-1.1363	2.0108	-3.1257	0.1945
5	1.74e-03		-1.1121	2.0218	-3.1571	0.2021
6	2.22e-04		-1.1225	2.0190	-3.1471	0.1993
7	2.97e-05		-1.1190	2.0203	-3.1510	0.2003
8	3.88e-06		-1.1203	2.0199	-3.1496	0.1999
9	5.15e-07		-1.1199	2.0200	-3.1501	0.2000
10	6.77e-08		-1.1200	2.0200	-3.1500	0.2000

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

=== 最終解 (ガウス・ザイデル法) ===
x[1] =   -1.12005
x[2] =    2.01999
x[3] =   -3.14995
x[4] =    0.19999

--- コード内での解の確認（検算: 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/sTh9Dr
language:
C (gcc 8.3)
created:
18 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!