プログラミング演習1
― 課題4 数値解析法 ―
1. Cによる数値解析のための準備
1.1 数学関数 Cには標準ライブラリの中に数学関数が用意されており、 ヘッダ <math.h> を インクルードすることで使用可能となる。
注: UNIX用のCコンパイラでは数学関数ライブラリのリンクが必要となる。 「undefined reference to `……'」などのリンクエラーが出たときはそれが原因である。ヘッダ <math.h> 内には使用できる関数が具体的に宣言されている。 Cで使用できる主な数学関数を以下に示す。
sin(x) 三角関数 sin(x) を計算する。 cos(x) 三角関数 cos(x) を計算する。 tan(x) 三角関数 tan(x) を計算する。 asin(x) 逆三角関数 sin-1(x) を計算する。(-π/2〜π/2) acos(x) 逆三角関数 cos-1(x) を計算する。(0〜π) atan(x) 逆三角関数 tan-1(x) を計算する。(-π/2〜π/2) atan2(y,x) 逆三角関数 tan-1(y/x) を計算する。(-π〜π) sinh(x) 双曲線関数 sinh(x) を計算する。 cosh(x) 双曲線関数 cosh(x) を計算する。 tanh(x) 双曲線関数 tanh(x) を計算する。 exp(x) 指数関数 ex を計算する。 log(x) 自然対数 ln(x) を計算する。(x > 0) log10(x) 10 を底とする対数 log10(x) を計算する。(x > 0) pow(x,y) べき乗 xy を計算する。 sqrt(x) 実数 x の平方根を求める。(x ≥ 0) ceil(x) 実数 x の小数部を切り上げた整数を返す。 floor(x) 実数 x の小数部を切り捨てた整数を返す。 fabs(x) 絶対値|x|を求める。 ldexp(x,n) x・2n を計算する。 frexp(x,&exp) 実数 x を仮数部と指数部(2のべき乗) exp に分解する。 modf(x,&ip) 実数 x を整数部 ip と小数部に分ける。 fmod(x,y) 実数の剰余(x % y に相当)を求める。 Microsoft Visual Studioのエディタ内でfabs等の関数名にカーソルを合わせ、 [F1]キーを押すと、そのキーワードのヘルプドキュメントが表示される。 おのおのの関数の詳しい使い方はここで調べることができる。 ページの下のほうの関連項目に「浮動小数点サポート」へのリンクがある。 ツールバーの[
練習1 pi.c (円周率を求める)ヘルプのお気に入り(P)]ボタン右側にある
(ヘルプのお気に入りに追加)ボタンで登録しておくとよい。
1: #include <stdio.h> 2: //#include <math.h> 3: 4: int main(void) 5: { 6: printf("PI = % .20g\n", 4.0*atan(1.0)); 7: return 0; 8: }2行目のコメントをはずさないと、どのようなエラーが出るか確認せよ。
ちなみに、通常は#define PI 3.14159265358979323846 または const double PI=3.1415926535897932384;をソースの #include の下に貼り付ければ、定数‘PI’に円周率を割り当てたことになる。
数学関数使用時は<math.h>をインクルード
1.2 数に関するユーティリティ ヘッダ<stdlib.h> には、いくつかの有用な数値変換関数や整数乱数発生関数が宣言されている。
atof(s) 文字列 s を double型浮動小数点数に変換する。 atoi(s) 文字列 s を int型整数に変換する。 atol(s) 文字列 s を long型整数に変換する。 rand() 0〜RAND_MAXの範囲の整数乱数を発生する。 srand(seed) rand() の種として seed をセットする。 abs(n) int型整数 n の絶対値 |n|を返す。 labs(n) long型整数 n の絶対値 |n|を返す。 練習2 rand10.c (乱数を発生する)
1: #include <stdio.h> 2: //#include <stdlib.h> 3: 4: int main(void) 5: { 6: int i; 7: 8: printf("RAND_MAX = % d\n", RAND_MAX); 9: for (i=0; i<10; i++) printf("% 2d: % d\n", i, rand()); 10: return 0; 11: }2行目のコメントをはずさないと、どのようなエラーが出るか確認せよ。
実行のみを繰り返したとき、結果が変わるかどうか確認せよ。
数値変換、乱数関数使用時は<stdlib.h>をインクルード
1.3 整数に関する制限 ヘッダ <limits.h> には整数型のサイズを表す定数が定義されている。
CHAR_BIT char型のビット数 CHAR_MAX char型の最大値 CHAR_MIN char型の最小値 INT_MAX int型の最大値 INT_MIN int型の最小値 LONG_MAX long型の最大値 LONG_MIN long型の最小値 SHRT_MAX short型の最大値 SHRT_MIN short型の最小値
UCHAR_MAX unsigned char型の最大値 UINT_MAX unsigned int型の最大値 ULONG_MAX unsigned long型の最大値 USHRT_MAX unsigned short型の最大値 Misrosoft Visual Studioのエディタ内でINT_MAX等の定数マクロ名にカーソルを合わせ、 [F1]キーを押すと「データ型定数」のヘルプドキュメントが現れる。 おのおのの定数の詳しい意味と具体的な数値はここで調べることができる。
練習3 chkint.c (標準的な整数のビットサイズと制限の確認)(ヘルプのお気に入りに追加) ボタンで登録しておくとよい。
1: #include <stdio.h> 2: //#include <limits.h> 3: 4: int main(void) 5: { 6: printf("int: %u bit, unsigned: %u bit\n", sizeof(int)*8, sizeof(unsigned)*8); 7: printf("INT_MAX = % d\n", INT_MAX); 8: printf("INT_MIN = % d\n", INT_MIN); 9: printf("UINT_MAX = %u\n", UINT_MAX); 10: return 0; 11: }2行目のコメントをはずさないと、どのようなエラーが出るか確認せよ。
1.4 浮動小数点数に関する定数 ヘッダ <float.h> には浮動小数点数に関する様々な定数が定義されている。
単精度 _FLT_RADIX 指数表現における基数b _FLT_ROUNDS 丸めのモード FLT_DIG 精度の10進けた数 FLT_EPSILON 計算機イプシロン FLT_MANT_DIG 浮動小数点数の有効数字部(_FLT_RADIXを基数とする)のけた数p FLT_MAX 表現可能な最大の有限浮動小数点数 FLT_MAX_EXP 表現可能な数の指数の最大値 FLT_MIN 最小の正規化された正の浮動小数点数 FLT_MIN_EXP 最小正規化数の指数の値
倍精度 _DBL_RADIX 指数表現における基数b _DBL_ROUNDS 丸めのモード DBL_DIG 精度の10進けた数 DBL_EPSILON 計算機イプシロン DBL_MANT_DIG 浮動小数点数の有効数字部(_DBL_RADIXを基数とする)のけた数p DBL_MAX 表現可能な最大の有限浮動小数点数 DBL_MAX_EXP 表現可能な数の指数部の最大値 DBL_MIN 最小の正規化された正の浮動小数点数 DBL_MIN_EXP 最小正規化数の指数部の値 Microsoft Visual Studioのエディタ内でDBL_MAX等の定数マクロ名にカーソルを合わせ、 [F1]キーを押すと「データ型定数」のヘルプドキュメントが現れる。 おのおのの定数の詳しい意味と具体的な数値はここで調べることができる。
練習4 chkflt.c (浮動小数点数のビットサイズと制限の確認)(ヘルプのお気に入りに追加) ボタンで登録しておくとよい。
1: #include <stdio.h> 2: //#include <float.h> 3: 4: int main(void) 5: { 6: printf("float: %u bit, double: %u bit\n", sizeof(float)*8, sizeof(double)*8); 7: printf("FLT_MAX = % e\n", FLT_MAX); 8: printf("FLT_MIN = % e\n", FLT_MIN); 9: printf("FLT_EPSILON = % e\n", FLT_EPSILON); 10: printf("DBL_MAX = % e\n", DBL_MAX); 11: printf("DBL_MIN = % e\n", DBL_MIN); 12: printf("DBL_EPSILON = % e\n", DBL_EPSILON); 13: return 0; 14: }2行目のコメントをはずさないと、どのようなエラーが出るか確認せよ。
1.5 関数引数 非線形方程式や数値積分、微分方程式のサブルーチンを作成する場合、 対象となる関数や方程式が替わるたびにサブルーチン内の関数名を書き換える必要がある。
しかし、サブルーチンの引数の一つに関数名を指定できるようにすれば、 サブルーチンを呼び出すメイン側に関数名の決定権を譲ることになり、 非常に汎用性の高いサブルーチンとなる。
それを実現するために使われるのが関数引数である。 (熊谷・玉城・白川「例題で学ぶC言語」の第11章 11.5 関数引数 pp.106-108 を参照のこと)
例えば後述する「二分法」のサブルーチン bisection で方程式 f (x)=0 の左辺の関数 f (x) を任意の名前で呼び出せるようにするには、 仮の関数名 f の前に*(アスタリスク)を付けてカッコで閉じた 関数名 (*f) (これを関数ポインタという)を用い、 これをサブルーチンの引数リスト内に関数の引数(型のみでよい)をつけて double (*f)(double) のように定義しておく。
double bisection(double (*f)(double), double x1, double x2) { 二分法の本体 }さらに、このサブルーチン内部で関数を呼び出すときはif (f(x1) * f(x3) < 0) x2 = x3;というように引数リスト内で定義した仮の関数名(この場合は f( ) )で呼び出すようにしておく。方程式の左辺の関数が実際には func1(x) という名前で定義されているなら、 main関数で二分法のサブルーチンを呼び出すときに
x = bisection(func1, x1, x2);と実際の関数名を引数なしで指定して呼び出せばよい。 練習5 funcp.c (関数引数の使用)
1: #include <stdio.h> 2: 3: double func1(double x) // f(x) = x を返す関数 4: { 5: printf("# func1(% g) is called", x); 6: return x; 7: } 8: 9: double func2(double x) // f(x) = x^2 を返す関数 10: { 11: printf("# func2(% g) is called", x); 12: return x*x; 13: } 14: 15: double func_p(double (*f)(double), double x) // 関数引数を用いたサブルーチン 16: { 17: double y; 18: 19: y = f(x); // ここでは仮の関数名 f() を用いる 20: printf(" in func_p().\n"); 21: return y; 22: } 23: 24: int main(void) 25: { 26: printf("func1(1) = % g\n", func_p(func1,1)); // func1 を指定して func_p を呼び出す 27: printf("func2(2) = % g\n", func_p(func2,2)); // func2 を指定して func_p を呼び出す 28: return 0; 29: }![]()
上記のプログラムを実行するとどうなるかを事前に予想せよ。 実際に実行してみて関数引数の使い方を確認せよ。
1.6 for文による等間隔実数列の生成 数値解析では一定間隔に刻んだ実数列を用いて計算することが多い。 以下は、柴田望洋「新版 明解C言語 入門編」の第7章にある、 0 から1 まで 0.01 間隔の実数列を生成するプログラム List 7-9 を少し手直ししたもの(forflt.c)である。
1: /* 2: 0 から1 まで 0.01 間隔の実数列を生成する。 3: */ 4: #include <stdio.h> 5: 6: int main(void) 7: { 8: float x; // double にするとどうなるか? 9: 10: for (x=0; x<=1; x+=0.01) printf(" x = % g\n", x); 11: 12: return 0; 13: }実行例x = 0 x = 0.01 x = 0.02 ... x = 0.979999 x = 0.989999 x = 0.999999このように、実数(単精度浮動小数点数)を直接、for文で増減させると、 刻み幅(上記の場合は0.01)を2進数にしたときの丸め誤差(表現誤差) が累積されて次第に値を正確に刻んでくれなくなる場合が多い。 このようなことを避けるためにも、for文は必ず整数で制御するべきである。 以下のプログラム(forflt2.c)は、 この考え方を適用した「明解C言語」の第7章 List 7-10と同じものである。1: /* 2: 0 から1 まで 0.01 間隔の実数列を生成する。(整数で制御) 3: */ 4: #include <stdio.h> 5: 6: int main(void) 7: { 8: int i; 9: float x; 10: 11: for (i=0; i<=100; i++) { 12: x = i / 100.0; // 100 にするとどうなるか? 100. では? 13: printf(" x = % g\n", x); 14: } 15: return 0; 16: }また、実数(浮動小数点数)宣言を double にすることで表現誤差が非常に小さくなる。 float を用いていて結果がおかしいときには double にするだけで直ることがある。練習6 (等間隔実数列の生成)
- 上記の2つのプログラム(forflt.cと forflt2.c)の結果を比較せよ。
- 2つめのプログラム(forflt2.c)で12行めの‘100.0’を‘100’にしてみよ。 また、‘100.’にしてみよ。
- 1つめのプログラム(forflt.c)で8行目の‘float’を‘double’にしてみよ。
for文は必ず整数で制御 実数定数は小数点をつけておく 実数は float よりも double で宣言
1.7 Excelによるグラフ描画 出力データをCSV形式で保存すれば、 Excelでただちに読み出すことができる。 CSV形式は一組のデータを単にコンマで区切り、組ごとに改行すればよい。 ファイル保存で名前をつける際、拡張子を“.csv” とする。 先頭行にデータの名前を入れておくと後で便利である。
以下の例は後述の2分法の出力ファイルで、 名前を“sample.csv”とする。
# n, x, err 1, 1.500000000000000e+000, 5.000000000000000e-001 2, 1.250000000000000e+000, 2.500000000000000e-001 3, 1.375000000000000e+000, 1.250000000000000e-001 4, 1.437500000000000e+000, 6.250000000000000e-002 5, 1.406250000000000e+000, 3.125000000000000e-002 6, 1.421875000000000e+000, 1.562500000000000e-002 7, 1.414062500000000e+000, 7.812500000000000e-003 8, 1.417968750000000e+000, 3.906250000000000e-003 9, 1.416015625000000e+000, 1.953125000000000e-003 10, 1.415039062500000e+000, 9.765625000000000e-004 11, 1.414550781250000e+000, 4.882812500000000e-004 12, 1.414306640625000e+000, 2.441406250000000e-004 13, 1.414184570312500e+000, 1.220703125000000e-004 14, 1.414245605468750e+000, 6.103515625000000e-005 15, 1.414215087890625e+000, 3.051757812500000e-005 16, 1.414199829101563e+000, 1.525878906250000e-005 17, 1.414207458496094e+000, 7.629394531250000e-006 18, 1.414211273193359e+000, 3.814697265625000e-006 19, 1.414213180541992e+000, 1.907348632812500e-006 20, 1.414214134216309e+000, 9.536743164062500e-007 # sqrt(2) = 1.414213657379150e+000
#(シャープ)記号はGNUPLOTで描画するときに コメント行となるようにするためのもので、 Excelでしか描画しないのであればなくてもよい。
- “sample.csv”のアイコンをダブルクリックする。
- ファイル名 *.csv には Excel 2003 が関連付けられているので、 Excel 2003 が起動し、sample.csv が読み込まれる。
- グラフを作成するために、データの範囲を選択する。 以下の操作により、データ範囲をすばやく選択できる。
- セル A1 が選択されていることを確認。
- [Ctrl] キーと [Shift] キーを同時に押しながら [→] キーを押す。
- [Ctrl] キーと [Shift] キーを同時に押しながら [↓] キーを押す。
- 数値の範囲が反転し、選択されていることを確認する。
- グラフウィザードボタン
を押す。
- グラフウィザードダイアログが表示されるので、以下の操作を行う。
- [グラフの種類(C)] で「散布図」を選択。
- [形式(T)] で左下の「データポイントを折れ線でつないだ散布図です。」を選択。
- [次へ > ] ボタンを2回クリックし、[タイトルとラベル] タブの 「X/数値軸(A)」 に“n” を、 「Y/数値軸(V)」 に“x, err” を記入する。
- [目盛線] タブを選択し、「Y/数値軸の目盛線(O)」 のチェックをはずす。
- [完了(F)] ボタンをクリック。
- 好みに応じて軸の目盛の範囲や凡例の位置、背景などを変更する。
- グラフをコピーして Word 2003 など他のアプリケーションに貼り付けることも可能である。
練習7 (Excelによるグラフの描画)
- まず“sample.csv”を右クリックし、 「対象をファイルに保存(A)...」で演習1課題4用のフォルダ(例えば“Z:\proen1\kadai4”)に保存する。
- 実際に上記の操作を行い、グラフを描いてみよ。
1.8 GNUPLOTによるグラフ描画 GNUPLOTは、
ボタンの「すべてのプログラム(P)
」の [
gnuplot] から起動できる。 [
gnuplot] を右クリックして「送る(N)」→「デスクトップ(ショートカットを作成)」 を選択しておけば、デスクトップのショートカットアイコンから起動できるようになる。
- 起動した後、[ファイル]→[ディレクトリの移動]で カレントディレクトリを、例えば‘Z:\proen1\kadai4’など、 CSVファイルがあるディレクトリ(フォルダに相当)に変える。
- ショートカットアイコンを右クリックしてプロパティ(R)を開き、 「作業フォルダ(S)」にCSVファイルがあるフォルダ名を登録しておけば、 ディレクトリ移動は不要になる。
- gnuplot> plot 'sample.csv' title 'x' with lines, 'sample.csv' using 1:3 title 'err' with lines と入力する。
- 意味は、‘sample.csv’(の1, 2列め)を使って‘x’という凡例titleで lineにてplotし、さらに‘sample.csv’の1列めと3列めを使って‘err’という 凡例titleでlineにてplotせよ、ということ。
- gnuplot> p 'sample.csv' t 'x' w l, 'sample.csv' u 1:3 t 'err' w l と省略できる。
- これにより、別ウィンドウにグラフが描画される。
- gnuplot> set xlabel 'n' と入力する。
- 意味はx軸のlabelを‘n’にset。
- gnuplot> se xl 'n' でもよい。
- gnuplot> set ylabel 'x, err' と入力する。
- 意味はy軸のlabelを‘x, err’にset。
- gnuplot> se yl 'x, err' でもよい。
- 変更分をグラフに反映させるため
gnuplot> replot と入力する。
- 意味は再plot。
- gnuplot> rep のみでもよい。
- グラフウィンドウ上部のタイトルバーを右クリックし、 [Options]から[Choose Font...], [Line Styles...]を選択し、 好みに応じてフォント(Font)や線種(Line Style)を変更する。
- グラフウィンドウ上部のタイトルバーを右クリックし、 [Options]から[Copy to Clipboard] を選択すれば、 他のアプリケーションに貼り付けることができる。
- 左下に現れるマウスポインタの座標表示を消すには
gnuplot> unset mouse と入力すればよい ( [M] キーを押してもよい)。- マウスの機能を復活させるには
gnuplot> set mouse と入力すればよい (あるいは再度 [M] キーを押す)。- 終了は
gnuplot> quit あるいは単に
gnuplot> q と入力する。詳しくは GNUPLOTのオンラインヘルプ を参照のこと。練習8 (GNUPLOTによるグラフの描画)
- デスクトップ上に gnuplot のショートカットアイコンを作成せよ。
- 演習でよく用いるフォルダを gnuplot ショートカットの[プロパティ(R)]の[作業フォルダ(S):]に登録せよ。
- 実際に上記の操作を行ってグラフを描いてみよ。
1.9 PADについて PADとは Problem Analysis Diagram の略である。 構造化フローチャートの一種で、処理の方向は「上から下へ」流れる。 具体的な処理内容が書かれた「処理箱」と呼ばれる 横長の四角い箱の左端を「縦線で結合(連接)」することによりプログラムのフローをあらわす。 Cにおいては、処理箱は文に、連結されたものはブロックに相当する。
![]()
「判断(選択)」、「反復」などの処理の中の文・ブロックはそれらの箱の右側に置き、横線で結ぶ。 これらの処理ではフローをいったん右側の文・ブロックに移す。 ブロックの場合は上からはじまり、一番下で左側に戻る。
「判断(if文)」は、箱の右辺が中折れしており、 上下の頂点から横に1本ずつ2本の分岐線を伸ばすことができ、 条件によっていずれか1つの分岐のみが実行される。
![]()
「前判定反復(while文)」は、箱の左側を二重線とする。
![]()
「後判定反復(repeat〜until文)」、「問題向き反復(for文)」は箱の右側を二重線とする。 Cでの後判定反復の実装にはrepeat〜until文の代わりにdo〜while文を用いる。 このとき、条件文が逆(否定)になることに注意する。
![]()
「定義」、「引用(参照)」は箱の内部に手続き名などを書き、両側あるいは境界を二重線とする。 定義は、箱の右側から二重線を横に出して処理ブロックに結ぶ。
![]()
引用(参照)は箱の両側または境界が二重線であること以外は処理箱と同様の扱いとなる。
![]()
PADのテンプレート
Wordなどで挿入できるファイル形式(jpeg, gif, tiffなど)にしておけば、 レポート作成に利用できる。
- ここ(pad.ppt)を右クリックし、 「対象をファイルに保存」で演習1_課題4用のフォルダに保存する。
- pad.pptを開き、ファイル内のPADをPowerPointの新規スライドに コピー・貼り付けし、グループ解除し、加工する。
- 加工し終わったらペイントツールを立ち上げ、新規ページを小さめにしておく。
- 完成したPADをPowerPoint上でドラッグ選択し、コピーする。
- ペイントツール(Adobe Photoshopなど)上で新規ページに貼り付け、 画像ファイルとして名前をつけて保存する。
また、Word上で書いたPADのサンプルも用意した。 次のように、判断(選択)の箱をフローチャートの「記憶データ」で代用している。 (上から処理箱を重ねればもう少し見栄えが良くなるが…。) レポートにはこちらの方法で描いてもよいことにする。
![]()
【目次】 | 【2.】| 【3.】| 【付録1】 | 【付録2】 | 【付録3】