プログラミング論 I
2008年06月26日,2008年07月03日
連立１次方程式を解く
消去法：ガウスの消去法
反復法：ヤコビ法
http://www.ns.kogakuin.ac.jp/~ct13140/Prog.2008/
F-1
概要
• 連立１次方程式を解く
– 消去法：ガウスの消去法
– 反復法：ヤコビ法
やや難しいです．
F-2
C言語
F-3
2重のfor文
for(i=0; i<2; i++){
for(j=0; j<3; j++){
printf("Hello\n");
}
}
for(i=0; i<2; i++){
j=0;
printf("Hello\n");
j=1;
printf("Hello\n");
j=2;
printf("Hello\n");
}
i=0;
j=0;
printf("Hello\n");
j=1;
printf("Hello\n");
j=2;
printf("Hello\n");
i=1;
j=0;
printf("Hello\n");
j=1;
printf("Hello\n");
j=2;
printf("Hello\n");F-4
2重のfor文
for(i=0; i<2; i++){
for(j=0; j<3; j++){
printf("Hello\n");
}
}
i=0;
for(j=0; j<3; j++){
printf("Hello\n");
}
i=1;
for(j=0; j<3; j++){
printf("Hello\n");
}
i=0;
j=0;
printf("Hello\n");
j=1;
printf("Hello\n");
j=2;
printf("Hello\n");
i=1;
j=0;
printf("Hello\n");
j=1;
printf("Hello\n");
j=2;
printf("Hello\n");F-5
2重のfor文
for(i=0; i<10; i++){
printf("#");
}
↓
##########
F-6
2重のfor文
for(j=0; j<3; j++){
for(i=0; i<10; i++){
printf("#");
}
printf("\n");
}
↓
##########
##########
##########
F-7
C言語：配列の宣言と使用
• 配列
int型が3個の配列．
data[0]～data[2]
が使用可能．
int data[3];
int i;
"0から始めて，3未満"
data[0] = 3;
でfor文を使うと
data[1] = 4;
やりやすい．
data[2] = 5;
for(i=0; i<3; i++){
printf("data[%d]", i);
printf("=%d\n", data[i]);
}
プログラム
data[0]=3
data[1]=4
data[2]=5
実行結果
F-8
C言語：2次元配列
int data[3][2];
int型が2×3=6個の配列．
int i, j;
data[0][0]～data[2][1]
data[0][0] = 3; data[0][1] = 4;
が使用可能．
data[1][0] = 5; data[1][1] = 6;
「長さ2個の配列」が3個．
data[2][0] = 7; data[2][1] = 8;
(長さ3の配列が2個でない)
for(i=0; i<3; i++){
for(j=0; j<2; j++){
data[0][0]=3
printf("data[%d][%d]", i, j);
data[0][1]=4
printf("=%d\n", data[i][j]);
}
data[1][0]=5
}
data[1][1]=6
プログラム
data[2][0]=7
data[2][1]=8
実行結果
F-9
配列のメモリ格納方法
int x[4][3];
は，
「長さ3の配列」が4個 であり，
「長さ4の配列」が3個 ではない？
計算機のメモリ上には，以下の順で格納されている．
ｘ[0][0] ｘ[0][1] ｘ[0][2] ｘ[1][0] ｘ[1][1] ｘ[1][2] ｘ[2][0] ｘ[2][1]
これは重要ではない.無理して理解する必要はない.
F-10
練習 0
int x[100]のx[0]～x[99]には値が格
納されている.
この100個をint y[100]のy[0]～
y[99]にコピーするプログラムは?
ヒント:
y[0] = x[0];
y[1] = x[1];
:
をfor文を使って作り上げれば良い.
F-11
解答 0
int i;
for(i=0; i<100; i++){
y[i] = x[i];
}
F-12
17
18
19
20
21
22
23
24
25
/* sum.c 行列の加算 */
26
#define COL 3
27
#define ROW 2
28
29
#include <stdio.h>
30
31
void print_matrix(double mat[COL][ROW]){
32
int r, c;
33
for(r=0; r<ROW; r++){
34
for(c=0; c<COL; c++){
35
printf("%6.2lf ", mat[r][c]);
36
}
37
printf("¥n");
38
}
39
}
40
41
42
行列の加算
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void main(){
double a[COL][ROW], b[COL][ROW];
double s[COL][ROW];
int r, c;
a[0][0]=0;
a[1][0]=3;
b[0][0]=6;
b[1][0]=9;
a[0][1]= 1;
a[1][1]= 4;
b[0][1]= 7;
b[1][1]=10;
a[0][2]=2;
a[1][2]=5;
b[0][2]=8;
b[1][2]=11;
printf("matrix a=\n");
print_matrix(a);
printf("matrix b=\n");
print_matrix(b);
/* begin : sum, s=a+b */
for(r=0; r<ROW; r++){
for(c=0; c<COL; c++){
s[r][c] = a[r][c] + b[r][c];
}
}
/* end : sum, s=a+b */
printf("matrix s=\n");
print_matrix(s);
}
F-13
行列の乗算
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
24
30
31
32
33
#define L 2
34
#define M 3
36
#define N 4
37
38
#include <stdio.h>
39
40
void print_matrix
41
(double *mat, int row, int col){
42
int r, c;
43
for(r=0; r<row; r++){
44
for(c=0; c<col; c++){
45
printf("%6.2lf ", *(mat+col*r+c));
46
}
48
printf("¥n");
49
}
50
}
51
/* pro.c 行列の乗算
(LxM 行列)×(MxN 行列) */
void main(){
double a[L][M], b[M][N], p[L][N];
int l, m, n;
a[0][0]= 0; a[0][1]= 1; a[0][2]= 2;
:省略
printf("matrix a=\n");
print_matrix(a, L, M);
printf("matrix b=\n");
print_matrix(b, M, N);
for(l=0; l<L; l++){
for(n=0; n<N; n++){
/* begin : calc p[l][n] */
double sum = 0.0;
for(m=0; m<M; m++){
sum += a[l][m]*b[m][n];
}
p[l][n] = sum;
/* end : calc p[l][n] */
}
}
printf("matrix s=\n");
print_matrix(p, L, N);
}
F-14
連立方程式の解法
F-15
連立一次方程式 (消去法)
x - 2y + 3z = 6 ‥ ‥(1a)
3x - 4y + 5z = 10 ‥ ‥(2a)
-2x + 5y - 7z = -13 ‥ ‥(3a)
式(2a)，式(3a)から x を消去
(2a)の両辺から，3×(1a)を引き x を消去
(3a)の両辺に，2×(1a)を足し x を消去
x - 2y + 3z = 6 ‥ ‥(1a)
2y - 4z = -8 ‥ ‥(2b) = (2a) - 3×(1a)
y - z = -1 ‥ ‥(3b) = (3a) + 2×(1a) F-16
連立一次方程式 (消去法)
式(3b)から y を消去
(3b)の両辺から，0.5×(2b)を引き x を消去
x - 2y + 3z = 6 ‥ ‥(1a)
2y - 4z = -8 ‥ ‥(2b)
z = 3 ‥ ‥(3c) = (3b) - 0.5×(2b)
z が求まった．
三角形になっている．
F-17
連立一次方程式 (消去法)
x - 2y + 3z = 6 ‥ ‥(1a)
2y - 4z = -8 ‥ ‥(2b)
z = 3 ‥ ‥(3c)
式(2b)から z を消去
(2b)の両辺に，4×(3c)を足し z を消去
x - 2y + 3z = 6 ‥ ‥(1a)
2y
= 4 ‥ ‥(2d)
y が求まった．
z = 3 ‥ ‥(3c) x - 2y + 3z = 6 ‥ ‥(1a)
(2d)の両辺を2で割る
y
= 2 ‥ ‥(2e)
z = 3 ‥ ‥(3c)
F-18
連立一次方程式 (消去法)
式(1a)から y,z を消去
(1a)の両辺に，2×(2e)を足し y を消去
(1a)の両辺から，3×(3c)を引いて z を消去
x
y
= 1 ‥ ‥(1f)
= 2 ‥ ‥(2e)
z = 3 ‥ ‥(3c)
解けた
F-19
連立方程式の記述方法/とらえ方
数学的な記述方法
x - 2y + 3z = 6 ‥ ‥(1a)
3x - 4y + 5z = 10 ‥ ‥(2a)
-2x + 5y - 7z = -13 ‥ ‥(3a)
係数だけを書いた
6
1 -2 3
3 -4 5
10
- 2 5 - 7 - 13
F-20
連立方程式の記述方法/とらえ方
行列として考える/記述することも可能
1 - 2 3 x 6
3 - 4 5 y 10
- 2 5 - 7 z - 13
数学的な記述方法
x - 2y + 3z = 6
3x - 4y + 5z = 10
-2x + 5y - 7z = -13
1 -2 3
x
6
A 3 - 4 5 ， x y ， b 10 として
- 2 5 - 7
z
- 13
連立方程式は，Ax=b と考える/記述することも可能．
Ax=b の両辺に，Aの逆行列A-1を左からかけると，
A-1Ax=A-1b → x=A-1b となり，x が求まる．
つまり，「逆行列を求める」と「連立方程式を解く」は非常に似ている．
F-21
消去法
• 未知数を1個ずつ消去していき，連立方程
式を解く手法
– ガウスの消去法，ガウス・ジョルダン法，LU分
解法，コレスキー法など
– 基本は，中学数学で習った消去法と同じ
– 基本的に(解くことのできる方程式であれば)必
ず解ける．
F-22
ガウス(Gauss)の消去法
• 前進消去フェーズ
– 1回の試行で，未知数を1個消していく．
– 左下が0(ゼロ)の三角形を作る．
• 後退代入フェーズ
– 1回の試行で，未知数の解を1個求める．
F-23
ガウスの消去法 (前進消去)
n個の未知数 x0, x1, ‥‥,xn-1のn元1次連立方程式
a0(0,0) x0 + a0(0,1) x1 + a0(0, n) -1xn -1 b0(0) 式 (0.0)
a1(,00) x0 + a1(,01) x1 + a1(,0n)-1xn -1 b1(0) 式 (0.1)
n元1次
連立方程式
an(0-)1,0 x0 + an(0-)1,1x1 + an(0-)1, n -1xn -1 bn(0-)1 式 (0.n - 1)
式(0.0)を用いて，式(0.1)～式(0.n-1)から，x0を消去する．
( 0)
具体的には，式(0.j)から，「 a (j01) a00
×式(0.0)」を引く．
ただし，j=1～n-1
( x)
( y)
注意： a0,0 と a0,0 は別の値．右上の添え字が別なら別の値．
右上の添え字は，変数の世代を表している．
F-24
ガウスの消去法 (前進消去)
結果，方程式は以下のように変形される．
a0( 0,0) x0 + a0( 0,1) x1 + a0( 0, 2) x2 + a0( 0,n)-1 xn -1 b0( 0) 式 (0.0)
a1(,11) x1 + a1(,12) x2 + a1(,1n)-1 xn -1 b1(1) 式 (1.1)
a2(1,1) x1 + a2(1, )2 x2 + a2(1,n) -1 xn -1 b2(1) 式 (1.2)
n-1元1次
連立方程式
an(1-)1,1 x1 + an(1-)1, 2 x2 + an(1-)1,n -1 xn -1 bn(1-)1 式 (1.n - 1)
( 0)
ただし， m j
a (j0,0)
a0(0,0)
，
a (j1,)k a (j0, k) - m(j0) a0(0, k)
b(j1) b(j0) - m(j0)b0(0)
j =1～n-1，k =1～n-1
( j が縦座標，k が横座標)
F-25
ガウスの消去法 (前進消去)
a0( 0,0) x0 + a0( 0,1) x1 + a0( 0, 2) x2 + a0( 0,n)-1 xn -1 b0( 0) 式 (0.0)
a1(,11) x1 + a1(,12) x2 + a1(,1n)-1 xn -1 b1(1) 式 (1.1)
a2(1,1) x1 + a2(1, )2 x2 + a2(1,n) -1 xn -1 b2(1) 式 (1.2)
n-1元1次
連立方程式
an(1-)1,1 x1 + an(1-)1, 2 x2 + an(1-)1,n -1 xn -1 bn(1-)1 式 (1.n - 1)
a0( 0,k)
n元方程式が，
n-1元方程式に
なった．
0
0
a (j0,k)
(1)
a
j ,k
0
：
0
式(1.1)～式(1.n-1)から，x0が消去された．
F-26
ガウスの消去法 (前進消去)
同様にして，
a0( 0,0) x0 + a0( 0,1) x1 + a0( 0, 2) x2 + a0( 0,n)-1 xn -1 b0( 0) 式 (0.0)
a1(,11) x1 + a1(,12) x2 + a1(,1n)-1 xn -1 b1(1) 式 (1.1)
a2( 2, 2) x2 + a2( 2,n)-1 xn -1 b2( 2) 式 (2.2)
an( 2-)1, 2 x2 + an( 2-)1,n -1 xn -1 bn( 2-1) 式 (2.n - 1)
(1)
ただし， m j
a (j1,1)
a1(,11)
，
n-2元1次
連立方程式
a (j2, k) a (j1,)k - m(j1) a1(,1k)
b(j2) b(j1) - m(j1)b1(1)
j =2～n-1，k =2～n-1
( j が縦座標，k が横座標)
F-27
ガウスの消去法 (前進消去)
式(2.2)～式(2.n-1)から，x1が消去された．n-2元方程式になった．
a0( 0, k)
a0( 0, k)
a1(,1k)
0
0
0 0
0
(1)
a j ,k
0 0
0
a (j2,k)
： ：
：
0 0
0
以下同様に，「xjを消去」を繰り返す．消去は，j =0～n-2 のn-1回行う．
0
0
0
：
0
a0( 0, k)
a1(,1k)
a 2( 2, k)
0
0 0
( 3)
a
j ,k
： ：
0 0
0
0
0
：
0
0
0 0
： ： 0
0 0 0 0
F-28
ガウスの消去法 (後退代入)
a0(0,0) x0
+ a0(0,1) x1 + a0(0,2) x2
a1(,11) x1 + a1(,12) x2
+
+
a0(0, n) -1xn -1
b0(0)
式(0.0)
+
+
a1(,1n) -1xn -1
b1(1)
式(1.1)
an( n--22, n) - 2 xn - 2
式(n-1,n-1)より，xn -1
bn( n--11)
an( n--1,1n) -1
+ an( n--22, n) -1xn -1 bn( n--22)
an( n--1,1n) -1xn -1 bn( n--11)
式(n - 2, n - 2)
式(n - 1, n - 1)
となり xn-1 が求まる．
xn-1を式(n-2,n-2)に代入して，xn - 2
bn( n--22) - an( n--22, n) -1xn -1
an( n--22, n) - 2
xn-2，xn-1を式(n-3,n-3)に代入して， xn-3を求める．
以下同様に，xj+1～xn-1を式( j, j)に代入し，xjを求める．
F-29
ガウスの消去法のまとめ
連立n元1次方程式の解き方
• 前進消去
L=0, 1,‥‥, n-2の順に，以下を行う．
j =L+1, L+2,‥‥, n-1の順に以下を行う．
m(jL)
a (jL, L)
aL( L, L)
k =L+1, L+2,‥‥, n-1順に以下を行う．
a (jL, k+1) a (jL, k) - m (jL ) a L( L, k)
b(jL +1) b(jL) - m(jL)bL( L)
F-30
ガウスの消去法のまとめ
連立n元1次方程式の解き方
• 後退代入
j=n-1, n-2,‥‥, 0の順に，以下を行う．
1 ( j ) n -1 ( j )
xj
b j - a j , k xk
a (j ,j )j
k j +1
F-31
ガウスの消去法 ポビット
( L)
• 第 L 段の消去で除算に使う aL, L をポビットと呼ぶ．
– ポビットが 0 である場合，計算を継続できない．
– ポビットの絶対値が非常に小さい場合，桁落ちな
どにより誤差が大きい可能性がある．
通常は，方程式を入れ替えて，
( L)
ai, L (ただし，i=L～n-1)の中で絶対値が最大のものを，
ポビットとして使用する．
F-32
前進消去
6
1 -2 3
3 -4 5
10
- 2 5 - 7 - 13
6
1 - 2 3
0 2 - 4 - 8
0 1 - 1 - 1
6
1 - 2 3
0 2 - 4 - 8
0 0
1
3
F-33
練習 1
• ガウスの消去法前進消去で，
上三角型に直せ
2 - 1 3 - 2
- 1 1 - 2 3
1 5
1
6
F-34
解答 1
2 - 1 3 - 2
- 1 1 - 2 3
1 5
1
6
2 - 1 3
1
1
0 2 - 2
0 0
5
- 2
2
7
2 - 1 3
1
1
0
2
2
0 11 - 1
2
2
- 2
2
7
F-35
反復法
• 計算を繰り返し，徐々に解に近い値を求め
ていく手法．
– ヤコビ法，ガウス・ザイデル法，SOR法など
– 必ずしも収束しない．
– 元数の大きい連立方程式の解を求めるのに
用いることが多い．
F-36
ヤコビ法 Jacobi法
n個の未知数 x0, x1, ‥‥,xn-1のn元1次連立方程式
a0,0 x0
a1,0 x0
a2,0 x0
an -1,0 x0
+
+
+
a0,1x1
a1,1x1
a2,1x1
+
+
+
a0,2 x2
a1,2 x2
a2,2 x2
+ an -1,1x1 + an -1,2 x2
+ + a0, n -1xn -1 b0
+ + a1, n -1xn -1 b1
+ + a2, n -1xn -1
b2
+ + an -1, n -1xn -1 bn -1
式(0)
式(1)
式(2)
式(n - 1)
式(0)をx0について解くと
x0
1
a0,0
b0 - a0,1x1 - a0,2 x2 - - a0,n-1xn-1
となる．
ただし x1, x2,…xn-1が未知なので結局，x0は求まらない？
F-37
ヤコビ法
同様にして
1
b0 - a0,1x1 - a0,2 x2 - - a0,n-1xn-1
x0
a0,0
1
x1
b1 - a1,0 x0 - a1,2 x2 - - a1, n -1xn -1
a1,1
xn -1
:
1
an -1, n -1
bn-1 - an-1,0 x0 - an-1,1x1 - - an -1,n-2 xn-2
となり xk を x0, x1,…,xk-1, xk+1,…,xn-1より 求めることができる．
ヤコビ法
x0, x1, …,xn-1に具体的な初期値を与え，
「x0, x1, … ,xk-1 ,xk+1, … ,xn-1からxkを求める」
ことを繰り返し，解に近づいていく．
F-38
ヤコビ法
n -1
1
( L)
b x (jL +1)
a
x
i, k k
j
a j, j
k 0, k j
第L世代の解より，
第L+1世代の解を求める
F-39
ヤコビ法の例
10x - y + z = 13
-2x +20y - z = -24
2x - 2 y +5z = 14
解は，
x=1, y=-1, z=2
3
2.5
2
解の候補
1.5
1
x
y
z
0.5
0
0
2
4
6
8
10
12
14
16
18
20
-0.5
-1
-1.5
世代(繰り返し回数)
F-40
ヤコビ法の例
3
2.5
2
解の候補
1.5
1
x
y
z
0.5
0
0
2
4
6
8
10
12
14
16
18
20
-0.5
-1
-1.5
世代(繰り返し回数)
10
2.5
1
0.1
0
5
10
15
20
0.01
1.5
x誤差
y誤差
z誤差
1
解の候補の誤差
解の候補の誤差
2
0.001
0.0001
x誤差
y誤差
z誤差
1E-05
1E-06
1E-07
1E-08
0.5
1E-09
1E-10
0
0
5
10
世代(繰り返し回数)
15
20
1E-11
世代(繰り返し回数)
F-41
© Copyright 2024 ExpyDoc
