コンピュータ物理学 2 第 2 回 (2015.10.9) 第 1 回 10/ 2( 金 ) ガイダンス 第 2 回 10/ 9( 金 ) 数値表現と誤差 第 3 回 10/16( 金 ) 第 4 回 10/23( 金 ) 数値微分 積分 第 5 回 10/30( 木 ) 第 6 回 11/13( 金 ) 第 7 回 11/20( 金 ) 常微分方程式 第 8 回 11/27( 金 ) 第 9 回 12/ 4( 金 ) 第 10 回 12/11( 金 ) 偏微分方程式 第 11 回 12/18( 金 ) 第 12 回 12/25( 金 ) 第 13 回 1/ 8( 金 ) モンテカルロ法 第 14 回 1/ 22( 金 ) 量子力学 第 15 回 1/ 29( 金 ) 第 16 回 2/ 5( 金 )
達成したい流れ まずどういうプログラムを作りたいのかを想定するのが大事 第 1 回の続き ; Guss 分布データの計算 生成 ソースプログラムをコンパイルして > ++ guss. o guss.ee プログラムを実行すると >./guss.ee 平均値 標準偏差 を聞かれるので適当な値を入力すると > Men vlue? : 100 > Stndrd devition? : 30 o オプション : コンパイルして出来上がる実行ファイル名を guss.ee に指定 こう入力すると m=100, sigm=30 をもとに = 20 から =220 まで d=2.4 刻みでガウス関数値を算出し 出力ファイル (guss.dt) に書き出すようなプログラムを作る この値をもとにガウス関数を計算した値のデータファイルが生成される (=100 を中心に ±120 の範囲で 100 個のデータ生成する 刻み幅 =240/100=2.4) ファイル guss.dt の中身 20.000000 0.000004 17.600000 0.000006 15.200001 0.000008 12.800001 0.000011 210.399841 0.000015 212.799835 0.000011 215.199829 0.000008 217.599823 0.000006 f()
/* Guss funtion */ #inlude <stdio.h> #inlude <mth.h> #define n 100 /* point numer */ doule guss_f(doule, doule m, doule sigm){ doule g; g = 1.0/(sqrt(2.0*M_PI*sigm*sigm))*ep( ( m)*( m)/(2*sigm*sigm)); return g; int min() { /* lol vriles */ doule m, sigm; /* men, stndrd devition */ doule min, m; /* lultion region */ doule, d, guss; int i; FILE *outfile; ここでは分割数は define で定義 π はヘッダファイル mth.h の中で定義されている 日本語変換で π とかしないこと 後で見て分かりやすいようにコメント文を付けること プログラムを見やすくするために インデント ( 字下げ ) を活用する 使う変数の型宣言 プログラムでは日本語変換された文字は使わないこと σ ではなく 半角文字の sigm を使う 等 /* min progrm */ printf("men vlue? "); snf("%lf", &m); printf("stndrd devition? "); snf("%lf", &sigm); outfile = fopen("guss.dt", "w"); min = m sigm*4; m = m + sigm*4; = min; d = (m min)/n; for(i=1; i<=n; i++){ guss = guss_f(, m, sigm); fprintf(outfile, %f %11.8f n",, guss); += d; 1 回ループごとに に d を足す = +d と書いても良い プログラムを走らせたら Men vlue? と端末上に表示させて そこでキーボードから手入力した数字を変数 m に代入する (sigm も同様 ) 計算する範囲を最小値 最大値を指定して決定する =min から始まって d ごとに を増やして 100 回関数 guss を計算する snf で数値を読み取る場合は doule 型なら %lf とする flot 型なら %f で良い guss.dt という名前でファイルを作り 書き込めるモード w で開く 計算刻み幅 d は全データが 100 個になるように (m min) を 100 で割る guss 関数を計算する部分 ここでは min 関数とは別に用意した関数 guss にて計算をさせる 渡す変数は計算に必要な, m, sigm の 3 つ ファイルに出力するのは と計算値 guss 浮動小数点表示は %f だが 表示させる桁を指定したい場合は %11.8 等とする ( 小数点部 8 桁 全部で 11 桁で揃える場合 )
C++ 仕様の書き方 : #inlude <iostrem> #inlude <fstrem> #inlude <mth.h> using nmespe std; doule guss_f(doule, doule m, doule sigm){ doule g; g = 1.0/(sqrt(2.0*M_PI*sigm*sigm))*ep( ( m)*( m)/(2*sigm*sigm)); return g; int min() { doule m, sigm, min, m, guss,, d; int i, n=100; out << "Men vlue?"; in >> m; out << "Stndrd devition?"; in >> sigm; std::ofstrem ofs("guss.dt"); out << "m=" << m << endl; out << "sigm=" << sigm << endl; min = m sigm*4; m = m+sigm*4; = min; d = (m min)/n; for(i=1; i<=n; i++){ guss = guss_f(, m, sigm); ofs << << " " << guss << std::endl; +=d;
gnuplot で guss.dt をプロットする (1) スタートメニューから gnuplot を起動する (2) 現在の作業ディレクトリに移動する (3) guss.dt をプロットする gnuplot では uni コマンド d (hnge diretory) が使えるただし 移動先のディレクトリ名は で囲む この授業での環境なら d z: 等 plot guss.dt と入力するとデータの 1 列目を横軸 2 列目を縦軸にしてプロットしてくれる
(4) もうひとつ作った guss0.dt も並べてプロットする
第 2-3 回のテーマ 扱える数値の大きさ 小ささ 数値計算で生じる誤差 例題 レポート問題を通してこれらを定量的に評価できるようにする 級数和の計算 計算結果ファイルから数値を読み取っての比較
数値表現 it : メモリの最小記憶単位 (0, 1 の 2 進数 ) 1 B N it を使って表せる整数は 2 N 個 符号の正負区別に 1it 使うので 2 N 1 個 yte : 1 yte = 1 B = 8 its 1 kb = 2 10 B = 1024 ytes 1yte = 8it は 2 8 =256 1 yte で英数字 1 文字を記憶できる = 061 = 01100001 = 062 = 01100010 等 int 型整数 : 32 it 2 31 までの整数を表現できる 2 31 2,147,483, 648 許される限度を超えた大きな数計算機が扱える精度より小さな数 overflow underflow limits.h 内で INT_MAX 等で定義されている
浮動小数点は? = (-1) [ 符号部 ] [ 仮数部 ] 2 [ 指数部 ] [ バイアス ] 例 : (-1) 0 0.1234567 10 2 flot 型 1 つの実数変数を表すのに 32 it を使う 1it: 符号部 8it: 指数部 23it: 仮数部 [ バイアス ]=127 2 126 ~ 2 127 を実現 10 進数で 38 桁 m 1 2 3 23 1 2 m2 2 m3 2 m23 2 仮数部の有効桁数 = 2 進数で 23 桁 = 10 進数で log 2 23 7 桁 例 : 3.5 仮数 : 符号 = +, 絶対値 = 2.5 ( 仮数部の定義は [ 仮数部 1]) 基数 : 2 で固定 (IEEE 方式 ) 指数 : 0 = (-1) 0 3.5 2 0 = (-1) 0 1.75 2 1 = (-1) 0 (1+0.75) 2 128 127 符号部 = 0 仮数部 = 0.75 = 0.11000000000000000000000 指数部 = 128 = 2 7 = 10000000
doule 型 1 つの実数変数を表すのに 64 it を使う ( 倍精度 ) 符号部 1it 指数部 11 it [ バイアス ]=1023 2 1022 ~ 2 1023 を実現 10 進数で 308 桁 m 仮数部 52 it 1 2 3 52 1 2 m2 2 m3 2 m52 2 仮数部の有効桁数 = 2 進数で 52 桁 = 10 進数で log 2 52 16 桁 例題 1: 数値範囲の確認 doule 型の変数 に対して 1~180 までの階乗を計算し コンピュータ上での数値を確認せよ 1!, 2!, 3!, 4!, 5!,..., 180! また 各時点の と同時に y=1/ と z=*y も表示して 計算結果がどう表示されるか確認せよ 階乗は for ループの各回ごとに 1 づつ増える整数をかければ良い *=i 等 ( I は for 文のパラメータ ) 大きな数値や小さな数値を表示させるときは指数で出すのが良い ; printf( %e, ) 端末上にこう出力させたい i=1 ; =1.0, y=1.0, z=1.0 i=2 ; =2.0, y=0.5, z=1.0 i=3 ; =6.0, y=0.166667, z=1.0 : : : : : : : : 無限大 (inf) や非数値 (NN) がどういう状況で出てきてしまうかわかっておく ( 今後も計算結果に inf / NN が出て悩むことがあるかもしれないので )
端末上にこう出力させたい i=1 ; =1.0, y=1.0, z=1.0 i=2 ; =2.0, y=0.5, z=1.0 i=3 ; =6.0, y=0.166667, z=1.0 : : : : : : : : そのためには : For 文 for(i=1; i<=180; i++) を使って i=1 の時点では = 1 i=2 の時点では = 2 1 i=3 の時点では = 3 2 1 等々を各回ごとに出力する 1 回前のループでの の値に今回のループでの i の値をかけることで 階乗が得られる = * i ( に i をかけた値を新たに に代入する ) これと同じ演算を *= i と書ける 当然 i が増えていくと! は爆発的に増えていく そこで浮動小数点が扱える範囲を超えるのを確認する 同時に の逆数も各ループで見てみる y=1/ として こちらは i が増えていくと 極めて小さい数値になっていく *y は定義上 1 だが や y が正常値をとれなくなると 値が定義できなくなるので それも見る =1.0; 0! = 1 なので 初期値は=1. for(i=1; i<=n; i++){ *= i; y=1/; z=*y; printf("i=%d =%e y=%e z=%e n", i,, y, z); 浮動小数点を指数型 (10^) で表示させると大きな数値や小さな数値を精度良く見れる (%e)
数値計算における誤差 コンピュータ上で計算を行う際 2 種類の誤差を考慮する必要がある (1) 近似誤差 コンピュータで問題を解くために 数学の厳密性を簡単化したために生じる誤差 無限級数 有限の和 無限小の間隔 有限のステップ e n N! n f ( ) d lim f N n 0 n n 0! n 1 0 N n0 n N 1 N N n0 f n N 1 N (2) 丸め誤差 計算機で扱う数値は有限個のビットで表されるので ( 有限な桁数を扱う ) 精度に限界があり 誤差が生じる 計算内容によっては有限桁数しか扱えないために生じる誤差が蓄積して大きな誤差を生じることがある 引き算で生じる誤差 ( 桁落ち ) 掛け算で生じる誤差
演算で生じる誤差 計算機で行われる計算は解析的な解に対する近似でしかない 2 つの数値の和を計算する場合を考える : 計算機の中での数値は添え字 を付けて厳密値と異なる近似値だと考えられる : 計算機内の有限桁数により生じる誤差を ε, ε とすると 1 1 誤差 1% なら ε = 0.01 誤差 0.1% なら ε = 0.001 一方 引算では 相対誤差は 1 の平均的誤差は と の誤差を重み付きで平均したもの だとすると ( 一方の誤差が大きい場合 ) 1 大きい方の誤差で和の誤差が決まる 引算で生じる相対誤差は 1 のときは の相対誤差が大きくなる ほとんど同じ大きさの数の間で差をとる 大きい桁の部分がなくなる 残りの小さい桁の部分に占める誤差の割合が相対的に大きくなる 桁落ち
桁落ちの例 大きい値同士の差で小さい値が生じるとき 5001 70.71775 と 5000 70. 71068 の差を計算することを 考える ともに有効数字 7 桁の精度を保持している 5001 5000 70.71775 70.71068 0.00707 有効数字は 3 桁 小さい値同士の差で小さい値が生じるとき 51 7.141428と 50 7. 071068 の差を計算することを 考える ともに有効数字 7 桁の精度を保持している 51 50 7.141428 7.071068 0.070360 有効数字は 5 桁
例題 2: 桁落ちが見られる関数計算 球 Bessel 関数 球 Bessel 関数 j l () を l = 2,3,4,, 10 まで計算 ファイル出力してグラフで確認する 漸化式と既知の j 0 (), j 1 () から求める 2l 1 jl 1( ) jl ( ) jl 1( ) j 0 ( ) sin, j ( ) 1 sin os 2 キーボード入力で計算したい次数 l を決めて計算できるようにする = 0.1 (=min) から始めて 刻み幅 0.1 で = 30 (=m) までの範囲で計算する for 文で を増加するループを作り 各回において決めた l に対する j l () を漸化式から計算する 漸化式を計算するのも欲しい l に達するまで for 文でループさせる 漸化式計算は関数として別たてするとよい doule up(doule, int l){ min 関数内で for(=min; <=m; +=d){ u = up(, l); fprintf(outfile, "%f %f n",, u); ( 上昇漸化式なので関数名を up としてる ) 漸化式は前の 2 つの項から求めるので 2 つ前の項を one 1 つ前の項を two として three = (2.0*n+1.0)/ * two one として新たな項を求める (n を増加させていく ) l が大きくなってくると原点付近で計算に誤差が生じてくる j 0 () j 1 () これを計算する (j 2 ()) 同様に j 10 () まで