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\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("回数  残差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.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.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.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/fJsOMH
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!