IT 入門 B2 ー 連立一次方程式（ 2 ） ー

ガウスの消去法の問題点 – 浮動小数点数の特殊な値 – 数学関数 ピボット選択つきガウスの消去法 演習 授 業 内 容 授 業 内 容

ガウスの消去法 前進消去と後退代入の組み合わせにより連立一次方程式の解 を求める手法 #include 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; }

ガウスの消去法 前回の授業で作成したガウスの消去法のプログラムを利用し て，以下の連立一次方程式の解を計算してみよう： ただし 結果はどうなるか？

実行例 のとき， を前回作成したガウスの消去法のプログラムを利 用して解を計算して見ると， ・ "nan" とは何か？ ・なぜ，正しい解が計算されないのか？ $./a.out x[0] = nan x[1] = nan x[2] = nan x[3] = nan

"nan" とは何か？ double 型のような浮動小数点を扱う際，多くの計算機では IEEE754 と呼ばれる 2 進浮動小数点の規格に従って処理が行わ れる． IEEE754 の規格書 IEEE754 では，計算過程において発生する可能性がある数の 中で，通常の浮動小数点数では表現できない特別な値が定義 されている． 例 ・ nan ・ inf

"nan" とは何か？ 非数 NaN （ Not a Number ） 0/0, 等の無効演算の結果として生成される特殊な数． "nan" と表示される． 計算過程のある部分で一度 NaN として評価されると，それ 以降はすべて NaN として評価される． 無限大 オーバーフローが起きた時に処理を続行可能とするため に導入された特殊な数． "inf" と表示される． 一度 inf として評価されても，最終的な結果は浮動小数点 数になる場合がある．

"nan" とは何か？ プログラム例 #include 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; }

 i kk ik kin ni kk ik nknik ki kk ik kkki b a a bbx a a a aax a a a aa                        1 1, 1,1,1 1, 1,1, なぜ，正しい解が計算されないのか？ 前進消去の第 k 列消去において の場合，それ以降の計算過 程において無効演算が発生する：

なぜ，正しい解が計算されないのか？ そこで， の中で最大のものを探す，つまり となる p を見つけ，第 k 行と第 p 行を入れ替える： この操作を部分ピボット選択という．

数学関数 指数関数や三角関数といった数学関数は，標準ライブラリと して準備されている． - 主な数学関数 - sqrt() : double 型数値の平方根を計算． pow() : double 型数値の任意の底と指数の累乗値を計算． exp() : double 型数値の e を底とする累乗値を計算． log() : double 型数値の自然対数を計算． log10() : double 型数値の常用対数を計算． sin(), cos(), tan() : 三角関数値を計算．引数の単位はラジアン． asin(), acos(), atan() : double 型数値の逆三角関数値を計算． abs() : int 型整数値の絶対値を計算． fabs() : double 型数値の絶対値を計算．

数学関数 数学関数を利用するには，その関数が定義されたヘッダファ イルをインクルードする必要がある： abs() : それ以外の関数 : また，ヘッダファイル を利用する場合は，コンパイ ル時に -lm オプションを付加する必要がある： #include … $ gcc –lm filename.c

数学関数 #include int main(void) { double a, b, c; a = sqrt(5.0); b = exp(10.0); c = sin(M_PI/6.0); ← M_PI : π の値のマクロ printf("a = %.20f\n", a); ← 小数点以下を 20 桁表示 printf("b = %e\n", b); ← 指数形式で表示 printf("c = %f\n", c); return 0; } プログラム例

部分ピボット選択を行うプログラムの作成 (1) の中で最大のものを探す． 絶対値 を計算する関数 ： 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 行を入れ替える． a[k][j] と a[p][j] (j=k+1,k+2,…,n-1), b[k] と b[p] をそれぞれ入れ替える．

#include int main(void) { … /* 前進消去 */ 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; } 部分ピボット選択付きガウスの消去法

部分ピボット選択付きガウスの消去法を実現するプログラム を完成させ，下記の連立一次方程式の解を計算してみよう： ただし 演 習 演 習