回答編集履歴
1
あ
answer
CHANGED
|
@@ -1,54 +1,151 @@
|
|
|
1
|
-
|
|
1
|
+
ガウスの消去法
|
|
2
2
|
|
|
3
|
-
行列 A の上三角行列
|
|
4
|
-
|
|
3
|
+
アルゴリズムについては [こちらのpdf](http://www.cs.tsukuba.ac.jp/~takahito/algebra1/system.pdf) がわかりやすいので参考にしてください。
|
|
4
|
+
ピボット選択は零割り防止の他、丸め誤差が大きくなる計算を防ぐ目的があります。 ピボット選択については、[こちら](http://nkl.cc.u-tokyo.ac.jp/12s/3D/pivoting.pdf) を参考にしてください。
|
|
5
5
|
|
|
6
|
-
|
|
6
|
+
```c
|
|
7
7
|
|
|
8
|
-
|
|
8
|
+
#include <math.h>
|
|
9
|
-
|
|
10
|
-
```c
|
|
11
9
|
#include <stdio.h>
|
|
12
10
|
|
|
11
|
+
#define n 4
|
|
12
|
+
|
|
13
|
+
void print_matrix(double A[][n])
|
|
14
|
+
{
|
|
15
|
+
for (int i = 0; i < n; ++i) {
|
|
16
|
+
for (int j = 0; j < n; ++j)
|
|
17
|
+
printf("%f ", A[i][j]);
|
|
18
|
+
printf("\n");
|
|
19
|
+
}
|
|
20
|
+
}
|
|
21
|
+
|
|
13
22
|
int main()
|
|
14
23
|
{
|
|
24
|
+
// clang-format off
|
|
15
|
-
|
|
25
|
+
double A[4][4] = {
|
|
16
26
|
{ 2, -1, 6, 3},
|
|
17
27
|
{ 2, -2, 4, 2},
|
|
18
28
|
{-1, 3, -4, 4},
|
|
19
29
|
{ 3, -2, 12, 12}};
|
|
30
|
+
// clang-format on
|
|
31
|
+
printf("original\n");
|
|
20
|
-
|
|
32
|
+
print_matrix(A);
|
|
21
33
|
|
|
22
|
-
// 上三角行列
|
|
23
|
-
printf("upper triangular matrix\n");
|
|
24
34
|
for (int i = 0; i < n; ++i) {
|
|
35
|
+
printf(
|
|
36
|
+
"step %d"
|
|
37
|
+
"\n----------------------------\n",
|
|
38
|
+
i);
|
|
39
|
+
|
|
40
|
+
// 1. <ピボット選択> k行i列の中で最大の要素を調べる。
|
|
41
|
+
// ----------------------------------
|
|
42
|
+
double max_val = fabs(A[i][i]);
|
|
43
|
+
int max_row = i;
|
|
25
|
-
for (int
|
|
44
|
+
for (int row = i + 1; row < n; ++row) {
|
|
26
|
-
|
|
45
|
+
if (fabs(A[row][i]) > max_val) {
|
|
46
|
+
max_val = fabs(A[row][i]);
|
|
27
|
-
|
|
47
|
+
max_row = row;
|
|
28
|
-
|
|
48
|
+
}
|
|
29
49
|
}
|
|
30
|
-
}
|
|
31
50
|
|
|
32
|
-
|
|
51
|
+
// 2. <ピボット選択> 現在の行と値が最も大きい行を入れ替える。
|
|
33
|
-
|
|
52
|
+
// ----------------------------------
|
|
53
|
+
if (i != max_row) {
|
|
34
|
-
|
|
54
|
+
for (int col = i; col < n; ++col) {
|
|
35
|
-
|
|
55
|
+
double tmp = A[max_row][col];
|
|
56
|
+
A[max_row][col] = A[i][col];
|
|
57
|
+
A[i][col] = tmp;
|
|
58
|
+
}
|
|
36
|
-
printf("%d ",
|
|
59
|
+
printf("swap rows: %d <--> %d\n", i, max_row);
|
|
37
|
-
if (j == i - 1)
|
|
38
|
-
|
|
60
|
+
print_matrix(A);
|
|
39
61
|
}
|
|
62
|
+
|
|
63
|
+
// 3. <前進消去>
|
|
64
|
+
// ----------------------------------
|
|
65
|
+
for (int row = i + 1; row < n; ++row) {
|
|
66
|
+
double tmp = A[row][i] / A[i][i];
|
|
67
|
+
printf("---> A[%d][%d] / A[%d][%d] = %f\n", row, i, i, i, tmp);
|
|
68
|
+
|
|
69
|
+
/* 第 i 行を -A[row][i] / A[i][i] 倍して、第 j 行に加える */
|
|
70
|
+
A[row][i] = 0; // A[row][i] -c * A[i][i] = 0 は確定しているので、
|
|
71
|
+
// 計算しなくてもよい。
|
|
72
|
+
for (int col = i + 1; col < n; ++col)
|
|
73
|
+
A[row][col] += -tmp * A[i][col];
|
|
74
|
+
printf("---> eliminate: row %d + %f * row %d\n", row, tmp, i);
|
|
75
|
+
print_matrix(A);
|
|
76
|
+
}
|
|
40
77
|
}
|
|
78
|
+
|
|
79
|
+
// 結果
|
|
80
|
+
printf("result\n");
|
|
41
|
-
|
|
81
|
+
print_matrix(A);
|
|
42
82
|
}
|
|
43
83
|
```
|
|
44
84
|
|
|
45
85
|
```
|
|
86
|
+
original
|
|
46
|
-
|
|
87
|
+
2.000000 -1.000000 6.000000 3.000000
|
|
88
|
+
2.000000 -2.000000 4.000000 2.000000
|
|
89
|
+
-1.000000 3.000000 -4.000000 4.000000
|
|
90
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
47
|
-
|
|
91
|
+
step 0
|
|
48
|
-
|
|
92
|
+
----------------------------
|
|
49
|
-
4
|
|
50
|
-
|
|
93
|
+
swap rows: 0 <--> 3
|
|
94
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
51
|
-
2
|
|
95
|
+
2.000000 -2.000000 4.000000 2.000000
|
|
96
|
+
-1.000000 3.000000 -4.000000 4.000000
|
|
52
|
-
-1 3
|
|
97
|
+
2.000000 -1.000000 6.000000 3.000000
|
|
98
|
+
---> A[1][0] / A[0][0] = 0.666667
|
|
99
|
+
---> eliminate: row 1 + 0.666667 * row 0
|
|
100
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
101
|
+
0.000000 -0.666667 -4.000000 -6.000000
|
|
102
|
+
-1.000000 3.000000 -4.000000 4.000000
|
|
103
|
+
2.000000 -1.000000 6.000000 3.000000
|
|
104
|
+
---> A[2][0] / A[0][0] = -0.333333
|
|
105
|
+
---> eliminate: row 2 + -0.333333 * row 0
|
|
106
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
107
|
+
0.000000 -0.666667 -4.000000 -6.000000
|
|
108
|
+
0.000000 2.333333 0.000000 8.000000
|
|
109
|
+
2.000000 -1.000000 6.000000 3.000000
|
|
110
|
+
---> A[3][0] / A[0][0] = 0.666667
|
|
111
|
+
---> eliminate: row 3 + 0.666667 * row 0
|
|
112
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
113
|
+
0.000000 -0.666667 -4.000000 -6.000000
|
|
114
|
+
0.000000 2.333333 0.000000 8.000000
|
|
115
|
+
0.000000 0.333333 -2.000000 -5.000000
|
|
53
|
-
|
|
116
|
+
step 1
|
|
117
|
+
----------------------------
|
|
118
|
+
swap rows: 1 <--> 2
|
|
119
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
120
|
+
0.000000 2.333333 0.000000 8.000000
|
|
121
|
+
0.000000 -0.666667 -4.000000 -6.000000
|
|
122
|
+
0.000000 0.333333 -2.000000 -5.000000
|
|
123
|
+
---> A[2][1] / A[1][1] = -0.285714
|
|
124
|
+
---> eliminate: row 2 + -0.285714 * row 1
|
|
125
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
126
|
+
0.000000 2.333333 0.000000 8.000000
|
|
127
|
+
0.000000 0.000000 -4.000000 -3.714286
|
|
128
|
+
0.000000 0.333333 -2.000000 -5.000000
|
|
129
|
+
---> A[3][1] / A[1][1] = 0.142857
|
|
130
|
+
---> eliminate: row 3 + 0.142857 * row 1
|
|
131
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
132
|
+
0.000000 2.333333 0.000000 8.000000
|
|
133
|
+
0.000000 0.000000 -4.000000 -3.714286
|
|
134
|
+
0.000000 0.000000 -2.000000 -6.142857
|
|
135
|
+
step 2
|
|
136
|
+
----------------------------
|
|
137
|
+
---> A[3][2] / A[2][2] = 0.500000
|
|
138
|
+
---> eliminate: row 3 + 0.500000 * row 2
|
|
139
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
140
|
+
0.000000 2.333333 0.000000 8.000000
|
|
141
|
+
0.000000 0.000000 -4.000000 -3.714286
|
|
142
|
+
0.000000 0.000000 0.000000 -4.285714
|
|
143
|
+
```
|
|
144
|
+
|
|
145
|
+
```
|
|
146
|
+
result
|
|
147
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
148
|
+
0.000000 2.333333 0.000000 8.000000
|
|
149
|
+
0.000000 0.000000 -4.000000 -3.714286
|
|
150
|
+
0.000000 0.000000 0.000000 -4.285714
|
|
54
151
|
```
|