回答編集履歴
1
あ
test
CHANGED
|
@@ -1,32 +1,52 @@
|
|
|
1
|
-
上三角行列、下三角行列の定義に基づき、for 文の条件式を設定すればよいです。
|
|
2
|
-
|
|
3
|
-
|
|
4
|
-
|
|
5
|
-
|
|
1
|
+
ガウスの消去法
|
|
6
|
-
|
|
2
|
+
|
|
3
|
+
|
|
4
|
+
|
|
7
|
-
|
|
5
|
+
アルゴリズムについては [こちらのpdf](http://www.cs.tsukuba.ac.jp/~takahito/algebra1/system.pdf) がわかりやすいので参考にしてください。
|
|
8
|
-
|
|
9
|
-
|
|
10
|
-
|
|
11
|
-
|
|
6
|
+
|
|
12
|
-
|
|
13
|
-
|
|
14
|
-
|
|
15
|
-
|
|
7
|
+
ピボット選択は零割り防止の他、丸め誤差が大きくなる計算を防ぐ目的があります。 ピボット選択については、[こちら](http://nkl.cc.u-tokyo.ac.jp/12s/3D/pivoting.pdf) を参考にしてください。
|
|
16
8
|
|
|
17
9
|
|
|
18
10
|
|
|
19
11
|
```c
|
|
20
12
|
|
|
13
|
+
|
|
14
|
+
|
|
15
|
+
#include <math.h>
|
|
16
|
+
|
|
21
17
|
#include <stdio.h>
|
|
22
18
|
|
|
23
19
|
|
|
24
20
|
|
|
21
|
+
#define n 4
|
|
22
|
+
|
|
23
|
+
|
|
24
|
+
|
|
25
|
+
void print_matrix(double A[][n])
|
|
26
|
+
|
|
27
|
+
{
|
|
28
|
+
|
|
29
|
+
for (int i = 0; i < n; ++i) {
|
|
30
|
+
|
|
31
|
+
for (int j = 0; j < n; ++j)
|
|
32
|
+
|
|
33
|
+
printf("%f ", A[i][j]);
|
|
34
|
+
|
|
35
|
+
printf("\n");
|
|
36
|
+
|
|
37
|
+
}
|
|
38
|
+
|
|
39
|
+
}
|
|
40
|
+
|
|
41
|
+
|
|
42
|
+
|
|
25
43
|
int main()
|
|
26
44
|
|
|
27
45
|
{
|
|
28
46
|
|
|
47
|
+
// clang-format off
|
|
48
|
+
|
|
29
|
-
|
|
49
|
+
double A[4][4] = {
|
|
30
50
|
|
|
31
51
|
{ 2, -1, 6, 3},
|
|
32
52
|
|
|
@@ -36,49 +56,109 @@
|
|
|
36
56
|
|
|
37
57
|
{ 3, -2, 12, 12}};
|
|
38
58
|
|
|
39
|
-
int n = 4;
|
|
40
|
-
|
|
41
|
-
|
|
42
|
-
|
|
43
|
-
//
|
|
59
|
+
// clang-format on
|
|
44
|
-
|
|
60
|
+
|
|
45
|
-
printf("
|
|
61
|
+
printf("original\n");
|
|
62
|
+
|
|
63
|
+
print_matrix(A);
|
|
64
|
+
|
|
65
|
+
|
|
46
66
|
|
|
47
67
|
for (int i = 0; i < n; ++i) {
|
|
48
68
|
|
|
69
|
+
printf(
|
|
70
|
+
|
|
71
|
+
"step %d"
|
|
72
|
+
|
|
73
|
+
"\n----------------------------\n",
|
|
74
|
+
|
|
75
|
+
i);
|
|
76
|
+
|
|
77
|
+
|
|
78
|
+
|
|
79
|
+
// 1. <ピボット選択> k行i列の中で最大の要素を調べる。
|
|
80
|
+
|
|
81
|
+
// ----------------------------------
|
|
82
|
+
|
|
83
|
+
double max_val = fabs(A[i][i]);
|
|
84
|
+
|
|
85
|
+
int max_row = i;
|
|
86
|
+
|
|
49
|
-
for (int
|
|
87
|
+
for (int row = i + 1; row < n; ++row) {
|
|
50
|
-
|
|
88
|
+
|
|
51
|
-
|
|
89
|
+
if (fabs(A[row][i]) > max_val) {
|
|
90
|
+
|
|
52
|
-
|
|
91
|
+
max_val = fabs(A[row][i]);
|
|
92
|
+
|
|
53
|
-
|
|
93
|
+
max_row = row;
|
|
54
|
-
|
|
94
|
+
|
|
55
|
-
|
|
95
|
+
}
|
|
56
96
|
|
|
57
97
|
}
|
|
58
98
|
|
|
99
|
+
|
|
100
|
+
|
|
101
|
+
// 2. <ピボット選択> 現在の行と値が最も大きい行を入れ替える。
|
|
102
|
+
|
|
103
|
+
// ----------------------------------
|
|
104
|
+
|
|
105
|
+
if (i != max_row) {
|
|
106
|
+
|
|
107
|
+
for (int col = i; col < n; ++col) {
|
|
108
|
+
|
|
109
|
+
double tmp = A[max_row][col];
|
|
110
|
+
|
|
111
|
+
A[max_row][col] = A[i][col];
|
|
112
|
+
|
|
113
|
+
A[i][col] = tmp;
|
|
114
|
+
|
|
115
|
+
}
|
|
116
|
+
|
|
117
|
+
printf("swap rows: %d <--> %d\n", i, max_row);
|
|
118
|
+
|
|
119
|
+
print_matrix(A);
|
|
120
|
+
|
|
121
|
+
}
|
|
122
|
+
|
|
123
|
+
|
|
124
|
+
|
|
125
|
+
// 3. <前進消去>
|
|
126
|
+
|
|
127
|
+
// ----------------------------------
|
|
128
|
+
|
|
129
|
+
for (int row = i + 1; row < n; ++row) {
|
|
130
|
+
|
|
131
|
+
double tmp = A[row][i] / A[i][i];
|
|
132
|
+
|
|
133
|
+
printf("---> A[%d][%d] / A[%d][%d] = %f\n", row, i, i, i, tmp);
|
|
134
|
+
|
|
135
|
+
|
|
136
|
+
|
|
137
|
+
/* 第 i 行を -A[row][i] / A[i][i] 倍して、第 j 行に加える */
|
|
138
|
+
|
|
139
|
+
A[row][i] = 0; // A[row][i] -c * A[i][i] = 0 は確定しているので、
|
|
140
|
+
|
|
141
|
+
// 計算しなくてもよい。
|
|
142
|
+
|
|
143
|
+
for (int col = i + 1; col < n; ++col)
|
|
144
|
+
|
|
145
|
+
A[row][col] += -tmp * A[i][col];
|
|
146
|
+
|
|
147
|
+
printf("---> eliminate: row %d + %f * row %d\n", row, tmp, i);
|
|
148
|
+
|
|
149
|
+
print_matrix(A);
|
|
150
|
+
|
|
151
|
+
}
|
|
152
|
+
|
|
59
153
|
}
|
|
60
154
|
|
|
61
155
|
|
|
62
156
|
|
|
63
|
-
//
|
|
157
|
+
// 結果
|
|
64
|
-
|
|
65
|
-
|
|
158
|
+
|
|
66
|
-
|
|
67
|
-
for (int i = 0; i < n; ++i) {
|
|
68
|
-
|
|
69
|
-
for (int j = 0; j < i; ++j) {
|
|
70
|
-
|
|
71
|
-
printf("%d ", a[i][j]);
|
|
72
|
-
|
|
73
|
-
if (j == i - 1)
|
|
74
|
-
|
|
75
|
-
|
|
159
|
+
printf("result\n");
|
|
76
|
-
|
|
77
|
-
|
|
160
|
+
|
|
78
|
-
|
|
79
|
-
}
|
|
80
|
-
|
|
81
|
-
r
|
|
161
|
+
print_matrix(A);
|
|
82
162
|
|
|
83
163
|
}
|
|
84
164
|
|
|
@@ -88,20 +168,134 @@
|
|
|
88
168
|
|
|
89
169
|
```
|
|
90
170
|
|
|
171
|
+
original
|
|
172
|
+
|
|
91
|
-
|
|
173
|
+
2.000000 -1.000000 6.000000 3.000000
|
|
174
|
+
|
|
92
|
-
|
|
175
|
+
2.000000 -2.000000 4.000000 2.000000
|
|
176
|
+
|
|
177
|
+
-1.000000 3.000000 -4.000000 4.000000
|
|
178
|
+
|
|
179
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
180
|
+
|
|
93
|
-
|
|
181
|
+
step 0
|
|
94
|
-
|
|
182
|
+
|
|
95
|
-
|
|
183
|
+
----------------------------
|
|
96
|
-
|
|
97
|
-
|
|
184
|
+
|
|
98
|
-
|
|
99
|
-
|
|
185
|
+
swap rows: 0 <--> 3
|
|
186
|
+
|
|
100
|
-
|
|
187
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
188
|
+
|
|
101
|
-
2
|
|
189
|
+
2.000000 -2.000000 4.000000 2.000000
|
|
190
|
+
|
|
102
|
-
|
|
191
|
+
-1.000000 3.000000 -4.000000 4.000000
|
|
192
|
+
|
|
103
|
-
-1 3
|
|
193
|
+
2.000000 -1.000000 6.000000 3.000000
|
|
194
|
+
|
|
104
|
-
|
|
195
|
+
---> A[1][0] / A[0][0] = 0.666667
|
|
196
|
+
|
|
197
|
+
---> eliminate: row 1 + 0.666667 * row 0
|
|
198
|
+
|
|
199
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
200
|
+
|
|
201
|
+
0.000000 -0.666667 -4.000000 -6.000000
|
|
202
|
+
|
|
203
|
+
-1.000000 3.000000 -4.000000 4.000000
|
|
204
|
+
|
|
205
|
+
2.000000 -1.000000 6.000000 3.000000
|
|
206
|
+
|
|
207
|
+
---> A[2][0] / A[0][0] = -0.333333
|
|
208
|
+
|
|
209
|
+
---> eliminate: row 2 + -0.333333 * row 0
|
|
210
|
+
|
|
211
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
212
|
+
|
|
213
|
+
0.000000 -0.666667 -4.000000 -6.000000
|
|
214
|
+
|
|
215
|
+
0.000000 2.333333 0.000000 8.000000
|
|
216
|
+
|
|
217
|
+
2.000000 -1.000000 6.000000 3.000000
|
|
218
|
+
|
|
219
|
+
---> A[3][0] / A[0][0] = 0.666667
|
|
220
|
+
|
|
221
|
+
---> eliminate: row 3 + 0.666667 * row 0
|
|
222
|
+
|
|
223
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
224
|
+
|
|
225
|
+
0.000000 -0.666667 -4.000000 -6.000000
|
|
226
|
+
|
|
227
|
+
0.000000 2.333333 0.000000 8.000000
|
|
228
|
+
|
|
229
|
+
0.000000 0.333333 -2.000000 -5.000000
|
|
230
|
+
|
|
105
|
-
|
|
231
|
+
step 1
|
|
232
|
+
|
|
106
|
-
|
|
233
|
+
----------------------------
|
|
234
|
+
|
|
235
|
+
swap rows: 1 <--> 2
|
|
236
|
+
|
|
237
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
238
|
+
|
|
239
|
+
0.000000 2.333333 0.000000 8.000000
|
|
240
|
+
|
|
241
|
+
0.000000 -0.666667 -4.000000 -6.000000
|
|
242
|
+
|
|
243
|
+
0.000000 0.333333 -2.000000 -5.000000
|
|
244
|
+
|
|
245
|
+
---> A[2][1] / A[1][1] = -0.285714
|
|
246
|
+
|
|
247
|
+
---> eliminate: row 2 + -0.285714 * row 1
|
|
248
|
+
|
|
249
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
250
|
+
|
|
251
|
+
0.000000 2.333333 0.000000 8.000000
|
|
252
|
+
|
|
253
|
+
0.000000 0.000000 -4.000000 -3.714286
|
|
254
|
+
|
|
255
|
+
0.000000 0.333333 -2.000000 -5.000000
|
|
256
|
+
|
|
257
|
+
---> A[3][1] / A[1][1] = 0.142857
|
|
258
|
+
|
|
259
|
+
---> eliminate: row 3 + 0.142857 * row 1
|
|
260
|
+
|
|
261
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
262
|
+
|
|
263
|
+
0.000000 2.333333 0.000000 8.000000
|
|
264
|
+
|
|
265
|
+
0.000000 0.000000 -4.000000 -3.714286
|
|
266
|
+
|
|
267
|
+
0.000000 0.000000 -2.000000 -6.142857
|
|
268
|
+
|
|
269
|
+
step 2
|
|
270
|
+
|
|
271
|
+
----------------------------
|
|
272
|
+
|
|
273
|
+
---> A[3][2] / A[2][2] = 0.500000
|
|
274
|
+
|
|
275
|
+
---> eliminate: row 3 + 0.500000 * row 2
|
|
276
|
+
|
|
277
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
278
|
+
|
|
279
|
+
0.000000 2.333333 0.000000 8.000000
|
|
280
|
+
|
|
281
|
+
0.000000 0.000000 -4.000000 -3.714286
|
|
282
|
+
|
|
283
|
+
0.000000 0.000000 0.000000 -4.285714
|
|
284
|
+
|
|
107
|
-
```
|
|
285
|
+
```
|
|
286
|
+
|
|
287
|
+
|
|
288
|
+
|
|
289
|
+
```
|
|
290
|
+
|
|
291
|
+
result
|
|
292
|
+
|
|
293
|
+
3.000000 -2.000000 12.000000 12.000000
|
|
294
|
+
|
|
295
|
+
0.000000 2.333333 0.000000 8.000000
|
|
296
|
+
|
|
297
|
+
0.000000 0.000000 -4.000000 -3.714286
|
|
298
|
+
|
|
299
|
+
0.000000 0.000000 0.000000 -4.285714
|
|
300
|
+
|
|
301
|
+
```
|