Microsoft PowerPoint コンピュータ物理2_第2回.pptx

コンピュータ物理学２

第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(金)      “

>  Mean value ? :  

100

> ./gauss.exe  

プログラムを実行すると、 平均値、標準偏差、を聞かれるので適当な値 を入力すると、

>  Standard deviation ? :  

30

この値をもとにガウス関数を計算した値のデータ ファイルが生成される (x=100を中心に ±120 の 範囲で100個のデータ生成する。刻み幅=240/100=2.4) ‐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 

・・・

ファイル

ga

uss.da

t

の中身

第1回の続き； Gauss 分布データの計算・生成

達成したい流れ

x

f(x)

こう入力すると m=100, sigma=30 をもとに x=‐20 から x=220 まで dx=2.4 刻みでガウス 関数値を算出し、出力ファイル (gauss.dat)  に書き出すようなプログラムを作る。

> c++ gauss.cc –o gauss.exe  

ソースプログラムをコンパイルして、 まずどういうプログラムを作りたいのかを想定するのが大事 ‐o オプション： コンパイルして出来上がる実行ファイル名 を gauss.exe に指定。

/*     Gauss function   */ #include <stdio.h> #include <math.h> #define n 100        /*  point number   */ double gauss_f(double x, double m, double sigma){ double g; g = 1.0/(sqrt(2.0*M_PI*sigma*sigma))*exp(‐(x‐m)*(x‐m)/(2*sigma*sigma)); return g; } int main() { /* local variables */ double  m, sigma;      /* mean, standard deviation  */ double  min, max;      /*  calculation region  */ double  x, dx, gauss;   int  i; FILE *outfile; /* main program */ printf("Mean value ? "); scanf("%lf", &m); printf("Standard deviation ? "); scanf("%lf", &sigma); outfile = fopen("gauss.dat", "w"); min = m ‐ sigma*4; max = m + sigma*4; x = min; dx = (max‐min)/n; for(i=1; i<=n; i++){ gauss = gauss_f(x, m, sigma); fprintf(outfile, “%f   %11.8f ¥n", x,  gauss); x += dx; } } プログラムを走らせたら“Mean value?” と端末上に 表示させて、そこでキーボードから手入力した数字 を変数 m に代入する。(sigmaも同様) scanf で数値を読み取る場合は double 型なら “%lf” とする。 float 型なら “%f” で良い。 π は ヘッダファイル math.h の中で 定義されている。日本語変換で “π” とかしないこと。 計算する範囲を最小値、最大値を指定して決定する x=min から始まって dx ごとに x を増やして 100 回関数 gauss を計算する gauss 関数を計算する部分。ここでは main 関数とは別に用意した 関数 gauss にて計算をさせる。渡す変数は計算に必要な x, m, sigma の３つ。 計算刻み幅 dx は全データが 100個 になるように (max – min) を 100 で割る ファイルに出力するのは x と 計算値 gauss。 浮動小数点表示は %f だが、表示させる 桁を指定したい場合は %11.8 等とする（小数点部 8桁、全部で 11桁で揃える場合） 使う変数の型宣言。プログラムでは日本語変換された文字は 使わないこと。”σ” ではなく、半角文字の ”sigma” を使う、等。 ・後で見て分かりやすいようにコメント文 を付けること。 ・プログラムを見やすくするために、インデント （字下げ）を活用する。 ここでは分割数はdefine で定義 １回ループごとに x に dx を足す。x = x+dx と 書いても良い。 gauss.dat という名前でファイルを作り、書き込めるモード “w” で開く。

#include <iostream> #include <fstream> #include <math.h> using namespace std; double gauss_f(double x, double m, double sigma){ double g; g = 1.0/(sqrt(2.0*M_PI*sigma*sigma))*exp(‐(x‐m)*(x‐m)/(2*sigma*sigma)); return g; } int main() { double m, sigma, min, max, gauss, x, dx; int i, n=100; cout << "Mean value ?"; cin >> m; cout << "Standard deviation ?"; cin >> sigma; std::ofstream ofs("gauss.dat"); cout << "m=" << m << endl; cout << "sigma=" << sigma << endl; min = m‐sigma*4; max = m+sigma*4; x = min; dx = (max‐min)/n; for(i=1; i<=n; i++){ gauss = gauss_f(x, m, sigma); ofs << x << " " << gauss << std::endl; x+=dx; } } C++ 仕様の書き方：

gnuplot で ‘gauss.dat’ をプロットする

(1) スタートメニューから gnuplot を起動する

(2) 現在の作業ディレクトリに移動する

・ gnuplot では unix コマンド “cd” (change directory) が使える ただし、移動先のディレクトリ名は “  “  で囲む。 ・ この授業での環境なら cd “z:” 等。 (3) gauss.dat をプロットする plot “gauss.dat” と入力すると データの 1列目を 横軸、 2列目を縦軸にしてプロッ トしてくれる

第2-3回のテーマ

・ 扱える数値の大きさ・小ささ

・ 数値計算で生じる誤差

例題・レポート問題を通してこれらを定量的に評価できる

ようにする。

・級数和の計算

・計算結果ファイルから数値を読み取っての比較

数値表現

・ bit   :  メモリの最小記憶単位

（0, 1 の2進数）

・ N bit を使って表せる整数は 2

N

_個

→ 符号の正負区別に 1bit 使うので 2

N‐1

_個

・ byte   :  

1 byte = 1 B = 8 bits

1 kB = 2

10

_{B = 1024 bytes}

1 byte で英数字１文字を記憶できる

1 B

int 型整数： 32 bit

2

31

_{までの整数を表現できる}

₂

31



_2,147,483,

₆₄₈

“a” = 0x61 = 01100001

“b” = 0x62 = 01100010  等

許される限度を超えた大きな数

_{→  overflow}

計算機が扱える精度より小さな数 →   underflow

1byte = 8bit は 2

8

₌₂₅₆

limits.h 内で INT_MAX 等で定義 されている

浮動小数点は？

(-1)

0

_{x 0.1234567 x 10}

2

float 型

1つの実数変数を表すのに 32 bit を使う

1bit: 符号部 8bit: 指数部 23bit:仮数部

仮数部の有効桁数 = 2進数で 23桁 = 10進数で log 223_{≈ 7 桁}

x

= (-1)

[符号部] 

_{× [仮数部] × 2}

[指数部]‐[バイアス] [バイアス]=127

2

‐126 

_{～ 2}

127

_を実現

→

10進数で 38桁

23 23 3 3 2 2 1 1

2

2

2

2

   

_

_

_

_

_

_

_



m

m

m

m



例：

例： 3.5 

仮数： 符号 = +,  絶対値 = 2.5 （仮数部の定義は [仮数部 ‐ 1]）

基数： 2 で固定 (IEEE方式)

指数： 0

x

= (-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}

double 型

1つの実数変数を表すのに 64 bit を使う（倍精度）

指数部 11 bit 

仮数部の有効桁数 = 2進数で 52桁 = 10進数で log 252 _{≈ 16 桁} 52 52 3 3 2 2 1 1

2

2

2

2

   

_

_

_

_

_

_

_



m

m

m

m



[バイアス]=1023

2

‐1022 

_{～ 2}

1023

_を実現

→

10進数で 308桁

符号部 1bit 

仮数部 52 bit

例題1：数値範囲の確認

double 型の変数 x に対して1～180 までの階乗を計算し、

コンピュータ上での数値を確認せよ。

!

180

,...,

!

5

,

!

4

,

!

3

,

!

2

,

!

1



x

また、各時点の x と同時に y=1/x と z=x*y も

表示して、計算結果がどう表示されるか確認せよ。

i=1 ;      x=1.0 ,  y=1.0,  z=1.0

i=2 ;      x=2.0 ,  y=0.5,  z=1.0

i=3 ;      x=6.0 ,  y=0.166667 , z=1.0

:       :       :       :

:       :       :       :

※ 階乗は for ループの各回ごとに 1づつ増える

整数をかければ良い。 x*=i 等( I はfor文のパラメータ)。 大きな数値や小さな数値を表示させるときは指数で 出すのが良い; printf(“%e”, x)。 無限大 (inf)や非数値 (NaN) がどういう 状況で出てきてしまうかわかっておく。 （今後も計算結果に inf / NaN が出て悩 むことがあるかもしれないので） 端末上にこう出力させたい

i=1 ;      x=1.0 ,  y=1.0,  z=1.0

i=2 ;      x=2.0 ,  y=0.5,  z=1.0

i=3 ;      x=6.0 ,  y=0.166667 , z=1.0

:       :       :       :

:       :       :       :

端末上にこう出力させたい

For 文 for(i=1; i<=180; i++) を使って

i=1 の時点では x = 1

i=2 の時点では x = 2×1

i=3 の時点では x = 3×2×1

等々

を各回ごとに出力する。

1回前のループでの x の値に今回のループでの

iの値をかけることで、階乗が得られる

x = x * i

（x に i をかけた値を新たに x に代入する）

これと同じ演算を x *= i と書ける。

そのためには：

x=1.0;

for(i=1; i<=n; i++){

x *= i;

y=1/x;

z=x*y;

printf("i=%d  x=%e y=%e z=%e¥n", i, x, y, z); 

}     

当然、i が増えていくと、x! は爆発的に増えていく。 そこで浮動小数点が扱える範囲を超えるのを確認 する。 同時に x の逆数も各ループで見てみる。 y=1/x として。 こちらは i が増えていくと、極めて小さい数値に なっていく。 x*y は定義上 1 だが、x や y が正常値をとれなく なると、値が定義できなくなるので、それも見る。 浮動小数点を指数型(ax10^b) で表示させると 大きな 数値や小さな数値を精度良く見れる。(%e)

0! = 1 なので、初期値はx=1.

数値計算における誤差

(1) 近似誤差

(2) 丸め誤差

コンピュータで問題を解くために、数学の厳密性を簡単化したために生じる誤差

・無限級数 → 有限の和

・無限小の間隔 → 有限のステップ





  





N n n n n x

n

x

n

x

e

0 0

!

!













































N n N n N

N

N

n

f

N

N

n

f

dx

x

f

0 0 1 0

1

1

lim

)

(

計算機で扱う数値は

有限個のビット

で表されるので（有限な桁数を扱う）

精度に限界があり、誤差が生じる。

計算内容によっては有限桁数しか扱えないために生じる誤差が蓄積

して大きな誤差を生じることがある。

・引き算で生じる誤差 （桁落ち）

・掛け算で生じる誤差

コンピュータ上で計算を行う際、2種類の誤差を考慮する必要がある。

演算で生じる誤差

c

b

a





2つの数値の和を計算する場合を考える：

計算機の中での数値は添え字c を付けて厳密値と 異なる近似値だと考えられる：

a

c



b

c



c

c



b

 

c

 

 

b c



c

b

c

b

c

b

c

a



1







1

















計算機内の有限桁数により生じる誤差を ε_b, ε_cとすると c b c

a

c

a

b

a

a

_{ 1}

_

_

_

_

c b

c

b







→ 相対誤差は だとすると（一方の誤差が大きい場合） b c

a

b

a

a

_{ 1}

_

_

誤差 1%   なら ε = 0.01 誤差 0.1%なら ε = 0.001 a の平均的誤差は b と c の誤差を重み付きで平均したもの 大きい方の誤差で和の誤差が決まる 計算機で行われる計算は解析的な解に対する近似でしかない。

c

b

a





a

_c



b

_c



c

_c 一方、引算では

a

c

a

b

a

a

c b c

 1









→ 引算で生じる相対誤差は b ≈ c   のときは a の相対誤差が大きくなる ほとんど同じ大きさの数の間で差をとる → 大きい桁の部分がなくなる → 残りの小さい桁の部分に占める誤差の割合が相対的に大きくなる

桁落ち

桁落ちの例

考える

の差を計算することを

と

₅₀₀₀

₇₀

_.

₇₁₀₆₈

71775

.

70

5001





ともに有効数字 7 桁の精度を保持している

00707

.

0

71068

.

70

71775

.

70

5000

5001









有効数字は3桁

考える

の差を計算することを

と

50

7

.

071068

141428

.

7

51





ともに有効数字 7 桁の精度を保持している

070360

.

0

071068

.

7

141428

.

7

50

51









有効数字は5桁

・大きい値同士の差で小さい値が生じるとき

・小さい値同士の差で小さい値が生じるとき

例題2：桁落ちが見られる関数計算 – 球 Bessel 関数 ‐

)

(

)

(

1

2

)

(

₁ 1

j

x

j

x

x

l

x

j

_l





_l



_l 球 Bessel 関数 j_l(x) を l = 2,3,4,・・・, 10 まで 計算・ファイル出力してグラフで確認する ・ 漸化式と既知の j₀(x), j₁(x) から求める ・ キーボード入力で計算したい次数 l を決めて計算できるようにする ・ x = 0.1 (=xmin) から始めて、刻み幅 0.1 で x = 30 (=xmax) までの範囲で計算する ・ ｆor文で x を増加するループを作り、各回において 決めた l に対する j_l(x) を漸化式から計算する。 漸化式を計算するのも欲しい l に達するまで for 文でループさせる。 ・ 漸化式計算は関数として別たてするとよい double up(double x, int l){

} main 関数内で for(x=min; x<=max; x+=dx){ u = up(x, l); fprintf(outfile, "%f  %f ¥n", x, u); } （上昇漸化式なので関数名を up としてる。） ・ 漸化式は前の２つの項から求めるので、２つ前の項を one １つ前の項を two として、 three = (2.0*n+1.0)/x * two ‐ one

として新たな項を求める(n を増加させていく)。 ・ l が大きくなってくると原点付近で計算に誤差が生じてくる。 2 1 0

cos

sin

)

(

,

sin

)

(

x

x

x

x

x

j

x

x

x

j







j

₀

(x)

j

₁

(x)

x

x

これを計算する(j

₂

(x))

同様に j

₁₀

(x) まで。