Transcript コンピュータ初級 (クラス2)

IT入門B2
ー 連立一次方程式（2） ー
授 業 内 容
• ガウスの消去法の問題点
– 浮動小数点数の特殊な値
– 数学関数
• ピボット選択つきガウスの消去法
• 演習
ガウスの消去法
前進消去と後退代入の組み合わせにより連立一次方程式の解
を求める手法
#include <stdio.h>
int main(void)
{
int i, j, k;
…
/* 前進消去 */
for (k=0; k<N-1; k++) {
for (i=k+1; i<N; i++) {
r = a[i][k]/a[k][k];
for (j=k+1; j<N; j++) {
a[i][j] = a[i][j] - a[k][j]*r;
}
b[i] = b[i] - b[k]*r;
}
}
/* 後退代入 */
for (i=N-1; i>=0; i--) {
r = 0.0;
for (j=i+1; j<=N-1; j++) {
r += a[i][j]*x[j];
}
x[i] = (b[i] - r)/a[i][i];
}
…
return 0;
}
ガウスの消去法
前回の授業で作成したガウスの消去法のプログラムを利用し
て，以下の連立一次方程式の解を計算してみよう：
Ax  b,
ただし
 2 4 1  3
0
 1  2 2 4 
10
, b   .
A
 4 2 3 5 
2
 5 4 3 1 
6


 
結果はどうなるか？
実行例
 2 4 1  3
0
 1  2 2 4 
10
, b   .
A
 4 2 3 5 
2
 5 4 3 1 
6


 
のとき， Ax  b を前回作成したガウスの消去法のプログラム
を利用して解を計算して見ると，
$ ./a.out
x[0] = nan
x[1] = nan
x[2] = nan
x[3] = nan
・"nan"とは何か？
・なぜ，正しい解が計算されないのか？
"nan"とは何か？
double型のような浮動小数点を扱う際，多くの計算機では
IEEE754と呼ばれる2進浮動小数点の規格に従って処理が行わ
れる．
<参考>IEEE754の規格書
http://www.validlab.com/754R/standards/754.html
IEEE754では，計算過程において発生する可能性がある数の中
で，通常の浮動小数点数では表現できない特別な値が定義さ
れている．
例 ・nan
・inf
"nan"とは何か？
非数 NaN（Not a Number）
0/0, 1 等の無効演算の結果として生成される特殊な数．
"nan"と表示される．
計算過程のある部分で一度NaNとして評価されると，それ
以降はすべてNaNとして評価される．
無限大
オーバーフローが起きた時に処理を続行可能とするため
に導入された特殊な数．
"inf"と表示される．
一度infとして評価されても，最終的な結果は浮動小数点
数になる場合がある．
"nan"とは何か？
プログラム例
#include <stdio.h>
int main(void)
{
double a, b;
a = 0.0/0.0;
printf("a = %f\n", a);
b = a + 1.0;
printf("b = %f\n", b);
a = 1.0/0.0;
printf("a = %f\n", a);
b = 1.0/a;
printf("b = %f\n", b);
return 0;
}
なぜ，正しい解が計算されないのか？
前進消去の第k列消去において akk  0 の場合，それ以降の計算
過程において無効演算が発生する：
akk xk + ak ,k +1xk +1 + L+ ak ,n1xn1  bk
ak +1,k xk + ak +1,k +1xk +1 + L+ ak +1,n1xn1  bk +1

an1,k xk + an1,k +1xk +1 + L+ an1,n1xn1  bn1
k +1  i  n 1




a
a
aik
ik
ik
 ai ,k +1  ak ,k +1

L
+
+



xk +1
ai ,n 1 ak ,n 1
xn 1 bi bk

akk 
akk 
akk





ai ,k +1
ai ,n 1
bi
なぜ，正しい解が計算されないのか？
そこで，aik k  i  n 1 の中で最大のものを探す，つまり
a pk  max aik
k i n1
となるpを見つけ，第k行と第p行を入れ替える：
akk xk + ak ,k +1xk +1 + L+ ak ,n1xn1  bk

a pk xk + a p,k +1xk +1 + L+ a p,n1xn1  bp

an1,k xk + an1,k +1xk +1 + L+ an1,n1xn1  bn1
a pk xk + a p,k +1xk +1 + L+ a p,n1xn1  bp

akk xk + ak ,k +1xk +1 + L+ ak ,n1xn1  bk

an1,k xk + an1,k +1xk +1 + L+ an1,n1xn1  bn1
この操作を部分ピボット選択という．
数学関数
指数関数や三角関数といった数学関数は，標準ライブラリと
して準備されている．
- 主な数学関数 sqrt() : double型数値の平方根を計算．
pow() : double型数値の任意の底と指数の累乗値を計算．
exp() : double型数値のeを底とする累乗値を計算．
log() : double型数値の自然対数を計算．
log10() : double型数値の常用対数を計算．
sin(), cos(), tan() : 三角関数値を計算．引数の単位はラジアン．
asin(), acos(), atan() : double型数値の逆三角関数値を計算．
abs() : int型整数値の絶対値を計算．
fabs() : double型数値の絶対値を計算．
数学関数
数学関数を利用するには，その関数が定義されたヘッダファ
イルをインクルードする必要がある：
abs() : <stdlib.h>
それ以外の関数 : <math.h>
#include <stdio.h>
#include <math.h>
…
また，ヘッダファイル<math.h>を利用する場合は，コンパイ
ル時に -lm オプションを付加する必要がある：
$ gcc
–lm
filename.c
数学関数
プログラム例
#include <stdio.h>
#include <math.h>
int main(void)
{
double a, b, c;
a = sqrt(5.0);
b = exp(10.0);
c = sin(M_PI/6.0);
printf("a = %.20f\n", a);
printf("b = %e\n", b);
printf("c = %f\n", c);
return 0;
}
← M_PI : πの値のマクロ
← 小数点以下を20桁表示
← 指数形式で表示
部分ピボット選択を行うプログラムの作成
(1) aik k  i  n 1 の中で最大のものを探す．
絶対値 a を計算する関数 ： fabs(a)
fabs(a[k+1][k]),fabs(a[k+2][k]), … ,fabs(a[n-1][k])
の中で最大の値fabs(a[p][k])を探し，pを記憶．
(2)第k行と第p行を入れ替える．
akk xk + ak ,k +1xk +1 + L+ ak ,n1xn1  bk

a pk xk + a p,k +1xk +1 + L+ a p,n1xn1  bp
a[k][j] と a[p][j] (j=k+1,k+2,…,n-1),
b[k] と b[p]
をそれぞれ入れ替える．
部分ピボット選択付きガウスの消去法
#include <stdio.h>
#include <math.h>
int main(void)
{
…
/* 前進消去 */
for (k=0; k<N-1; k++) {
/*
部分ピボット選択
/* 後退代入 */
for (i=N-1; i>=0; i--) {
r = 0.0;
for (j=i+1; j<=N-1; j++) {
r += a[i][j]*x[j];
}
x[i] = (b[i] - r)/a[i][i];
}
*/
for (i=k+1; i<N; i++) {
r = a[i][k]/a[k][k];
for (j=k+1; j<N; j++) {
a[i][j] = a[i][j] - a[k][j]*r;
}
b[i] = b[i] - b[k]*r;
}
}
…
return 0;
}
演 習
部分ピボット選択付きガウスの消去法を実現するプログラム
を完成させ，下記の連立一次方程式の解を計算してみよう：
Ax  b,
ただし
 2 4 1  3
0
 1  2 2 4 
10
, b   .
A
 4 2 3 5 
2
 5 4 3 1 
6


 