.
Ｃプログラミング (2014)
.
第5回
− 宿題１の解説と演習１ −
松田七美男
2014 年 10 月 23 日
. 宿題 1：概要
.
ガウスの消去法の確認
.
3 次の Hilbert 行列（係数が hij = 1/(i + j − 1)）を係数行列
とする 3 元連立一次方程式を Gauss の消去法を用いて解き
.なさい．
定数 bi は 1 ≤ bi ≤ 10 の整数となるように，疑似乱数発
生関数 rand() を用いて main 関数内で定義しなさい．
rand() を使うには，初期化関数 srand() を適当な整数
seed を引数に srand(seed) と呼び出す必要がありま
す．この seed には，各自の学籍番号を用いなさい．
例） 13YB034* ⇒ seed=34
Ax が b に等しいことを確かめる検算結果も表示しな
さい．
. 個人毎の方程式の定義
係数行列 A と定数ベクター b の要素を乱数を用いて定義
する．
#define N
3
#define SEED 300
...
srand(SEED);
for (i = 0; i < N; i++) {
b[i] = rand() % 10 + 1;
for (j = 0; i < N; j++) {
A[N*i + j] = 1.0/(i + j +1);
}
}
...
. ソース例
#include "lineq.h"
#define N
3
#define SEED 300
copy_mat(A, Ai, N, N);
print_mat("A =", A, N, N);
print_vec("b =", b, N);
mat_gauss(A, b, x, N);
print_vec("x =", x, N);
calc_ax(Ai, x, Ax, N);
print_vec("Ax = b?", Ax, N);
int main(int argc, char *argv[])
{
int i, j;
Mat A, Ai;
Vec b, x, Ax;
srand(SEED);
A = calloc( N*N, sizeof(double));
Ai = calloc( N*N, sizeof(double));
b = calloc( N, sizeof(double));
x = calloc( N, sizeof(double));
Ax = calloc( N, sizeof(double));
for (i = 0; i < N; i++) {
b[i] = rand() % 10 + 1;
for (j = 0; j < N; j++) {
A[N*i+j] = 1.0/(i+j+1);
}
}
.
free(A);
free(Ai);
free(b);
free(x);
free(Ax);
return 0;
}
. ソース例
#include "lineq.h"
#define N
3
#define SEED 300
copy_mat(A, Ai, N, N);
print_mat("A =", A, N, N);
print_vec("b =", b, N);
mat_gauss(A, b, x, N);
print_vec("x =", x, N);
calc_ax(Ai, x, Ax, N);
print_vec("Ax = b?", Ax, N);
int main(int argc, char *argv[])
{
int i, j;
Mat A, Ai;
Vec b, x, Ax;
srand(SEED); 乱数の初期化
A = calloc( N*N, sizeof(double));
Ai = calloc( N*N, sizeof(double));
b = calloc( N, sizeof(double));
x = calloc( N, sizeof(double));
Ax = calloc( N, sizeof(double));
for (i = 0; i < N; i++) {
b[i] = rand() % 10 + 1;
for (j = 0; j < N; j++) {
A[N*i+j] = 1.0/(i+j+1);
}
}
.
free(A);
free(Ai);
free(b);
free(x);
free(Ax);
return 0;
}
. ソース例
#include "lineq.h"
#define N
3
#define SEED 300
copy_mat(A, Ai, N, N);
print_mat("A =", A, N, N);
print_vec("b =", b, N);
mat_gauss(A, b, x, N);
print_vec("x =", x, N);
calc_ax(Ai, x, Ax, N);
print_vec("Ax = b?", Ax, N);
int main(int argc, char *argv[])
{
int i, j;
Mat A, Ai;
Vec b, x, Ax;
srand(SEED); 乱数の初期化
A = calloc( N*N, sizeof(double));
Ai = calloc( N*N, sizeof(double));
b = calloc( N, sizeof(double));
x = calloc( N, sizeof(double));
Ax = calloc( N, sizeof(double));
free(A);
free(Ai);
free(b);
free(x);
free(Ax);
return 0;
}
for (i = 0; i < N; i++) {
b[i] = rand() % 10 + 1;
for (j = 0; j < N; j++) {
A[N*i+j] = 1.0/(i+j+1);
aij
}
i,j が 0 から始まっている点に注意
}
.
. ソース例
#include "lineq.h"
#define N
3
#define SEED 300
copy_mat(A, Ai, N, N);
print_mat("A =", A, N, N);
print_vec("b =", b, N);
ガウスの消去法の呼出し
mat_gauss(A, b, x, N);
int main(int argc, char *argv[])
print_vec("x =", x, N);
{
calc_ax(Ai, x, Ax, N);
int i, j;
print_vec("Ax = b?", Ax, N);
Mat A, Ai;
Vec b, x, Ax;
free(A);
free(Ai);
srand(SEED); 乱数の初期化
free(b);
A = calloc( N*N, sizeof(double));
free(x);
Ai = calloc( N*N, sizeof(double));
free(Ax);
b = calloc( N, sizeof(double));
x = calloc( N, sizeof(double));
return 0;
Ax = calloc( N, sizeof(double));
}
for (i = 0; i < N; i++) {
b[i] = rand() % 10 + 1;
for (j = 0; j < N; j++) {
A[N*i+j] = 1.0/(i+j+1);
aij
}
i,j が 0 から始まっている点に注意
}
.
. 演習 1：課題
.
LU
. 分解法により逆行列を得る関数
LU 分解去を用いて，行列 A の逆行列 X = A−1 を
.求める関数 inv mat を作成しなさい．
宣言：inv mat(Mat A, int N)
テストプログラム lineq invpascal.c を作成して，5
次の Pascal 行列（Pij = Pi−1,j + Pi,j−1 , P1,j = Pi,1 = 1）
の逆行列を求めなさい．
. 逆行列
.
i 番目
.
とすると，
z}|{
bi = (0, · · · , 1 , · · · , 0)T
xi : Axi = bi の解
. 逆行列
.
i 番目
.
z}|{
bi = (0, · · · , 1 , · · · , 0)T
xi : Axi = bi の解
とすると，逆行列は以下のように計算される．
]
[
A−1 = A−1 I = A−1 b1 , b2 , · · · , bn
[
]
= A−1 b1 , A−1 b2 , · · · , A−1 bn
[
]
= x1 , x2 , · · · , xn
. 逆行列
.
i 番目
.
z}|{
bi = (0, · · · , 1 , · · · , 0)T
xi : Axi = bi の解
とすると，逆行列は以下のように計算される．
]
[
A−1 = A−1 I = A−1 b1 , b2 , · · · , bn
[
]
= A−1 b1 , A−1 b2 , · · · , A−1 bn
[
]
= x1 , x2 , · · · , xn
bi 毎に xi を計算し，それらを順番に並べればよい．
. 演習 1：実行画面
.
テストプログラム inv の実行画面
.
$ ./ex2014A-1
A =
1
1
1
2
1
3
1
4
1
5
1
3
6
10
15
1
4
10
20
35
1
5
15
35
70
X =
5
-10
10
-5
1
-10
30
-35
19
-4
10
-35
46
-27
6
-5
19
-27
17
-4
1
-4
6
-4
1
AX = I?
1
0
0
0
. 0
0
1
0
0
0
0
0
1
0
0
0
0
0
1
0
0
0
0
0
1
. 演習 1：実行画面
.
テストプログラム inv の実行画面
.
$ ./ex2014A-1
A =
1
1
1
2
1
3
1
4
1
5
1
3
6
10
15
1
4
10
20
35
1
5
15
35
70
X =
5
-10
10
-5
1
-10
30
-35
19
-4
10
-35
46
-27
6
-5
19
-27
17
-4
1
-4
6
-4
1
AX = I?
1
0
0
0
. 0
0
1
0
0
0
0
0
0
1
0
0
0
0
0
1
単位行列
0
0
1
0
0
. 演習課題提出に関する注意
Waseda-net Web メールサービスを用いて提
出しなさい．
宛先は [email protected]
件名には，学生番号（英数半角），名前，演習１をこの
順に記述しなさい．
例）1Y13B999，△○ ×□，演習 1
〆切は，10 月 23 日 10:30
送信控えの保存を忘れずに
ソース（inv mat.c と lineq invmat.c）と実行
結果の両方が必要．添付不可．
ソースには空白を挿入したり，字下げをほど
こすなど読みやすくなる工夫をしなさい．