コンピュータ初級 (クラス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,
ただし
4
1  3
2
0
 1  2 2

10
4
, b   .
A
2 3 5 
4
2
 5 4 3 1 
6


 
結果はどうなるか?
実行例
4
1  3
2
0
 1  2 2

10
4
, b   .
A
2 3 5 
4
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 +1 xk +1 + L + ak ,n 1 xn 1  bk
ak +1,k xk + ak +1,k +1 xk +1 + L + ak +1,n 1 xn 1  bk +1

an 1,k xk + an 1,k +1 xk +1 + L + an 1,n 1 xn 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 +1 xk +1 + L + ak ,n 1 xn 1  bk

a pk xk + a p ,k +1 xk +1 + L + a p ,n 1 xn 1  bp

an 1,k xk + an 1,k +1 xk +1 + L + an 1,n 1 xn 1  bn 1
a pk xk + a p ,k +1 xk +1 + L + a p ,n 1 xn 1  bp

akk xk + ak ,k +1 xk +1 + L + ak ,n 1 xn 1  bk

an 1,k xk + an 1,k +1 xk +1 + L + an 1,n 1 xn 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 +1 xk +1 + L + ak ,n 1 xn 1  bk

a pk xk + a p ,k +1 xk +1 + L + a p ,n 1 xn 1  b p
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,
ただし
4
1  3
2
0
 1  2 2

10
4
, b   .
A
2 3 5 
4
2
 5 4 3 1 
6


 