Document

コンピュータ物理学２
第5回
(2015.10.30)
第1回
第2回
第3回
第4回
第5回
第6回
第7回
第8回
10/ 2(金) ガイダンス
10/ 9(金) 数値表現と誤差
10/16(金) “
10/23(金) 数値微分・積分
10/30(木) “
11/13(金) “
11/20(金) 常微分方程式
11/27(金)
“
第9回
第10回
第11回
第12回
第13回
第14回
第15回
第16回
12/ 4(金) “
12/11(金) 偏微分方程式
12/18(金) “
12/25(金) “
1/ 8(金) モンテカルロ法
1/ 22(金) 量子力学
1/ 29(金) “
2/ 5(金) “
解答例
課題 No.1 ：特殊関数の数値計算と誤差
(1) 例題２の結果をグラフで示してまとめよ。
次数 l = 2, 3, 4 について３つの計算結果
を右図に示す。（正しく計算できている）
次数 l = 9, 10 についの計算結果。
j2(x)
j3(x)
j4(x)
j9(x)
j10(x)
l=15 について、最初の 6行の数値を見てみると：
0.100000
0.200000
0.300000
0.400000
0.500000
0.600000
3.239599e+14
-4.051041e+09
1.519395e+06
-3.377090e+05
-6.500883e+03
4.624043e+02
Bessel関数は次数が上がると、より 0 に近づく
値をとるので、x=0~1 程度で 0.1 より大きな値
をとってるのは正しく計算できてないことをあら
わしている
j(x)~0 付近を拡大して見てみると、
j9(x)
j10(x)
j15(x)
double 型と float 型での計算値の違い
double型での計算
float型での計算
x
0.100000
0.200000
0.300000
0.400000
0.500000
0.600000
0.700000
0.800000
0.900000
1.000000
1.100000
1.200000
1.300000
1.400000
1.500000
double 型
での結果
5.909052e-09
3.067328e-08
2.329596e-07
9.790376e-07
2.977466e-06
7.377563e-06
1.586612e-05
3.075513e-05
5.505920e-05
9.256116e-05
1.478647e-04
2.264326e-04
3.346105e-04
4.796341e-04
6.696206e-04
float 型
での結果
-1.389741e-03
-7.208532e-04
7.136416e-04
-3.794678e-04
-1.573712e-04
3.550287e-05
1.258652e-05
5.608677e-05
5.353427e-05
9.220998e-05
1.512924e-04
2.286922e-04
3.328141e-04
4.778758e-04
6.696834e-04
x=0.1 では 6桁も結果が食い違う
x=1.4 では 3.7% の食い違いまで収まっている
（それでも無視はできないが）
桁落ちについて考えてみる
x=0.1 における値に注目する
有効桁数が 8桁ある場合
j0 ( x  0.1)  9.9833416  10 1
j1 ( x  0.1)  3.3300012  10  2
3
 j1 ( x  0.1)  j0 ( x  0.1)
0 .1
 9.9900036  10 1  9.9833416  10 1
j2 ( x  0.1) 
有効桁数が 4桁しかない場合
j0 ( x  0.1)  9.983  10 1
j1 ( x  0.1)  3.330  10  2
3
 j1 ( x  0.1)  j0 ( x  0.1)
0.1
 9.990  10 1  9.983  10 1
j2 ( x  0.1) 
 0.006662  10 1
 0.007  10 1
 6.662  10  4
 7  10  4
5
 j2 ( x  0.1)  j1 ( x  0.1)
0.1
 3.331  10  2  3.3300012  10  2
j3 ( x  0.1) 
5
 j2 ( x  0.1)  j1 ( x  0.1)
0.1
 3.50  10  2  3.330  10  2
j3 ( x  0.1) 
 0.0009988  10  2
 0.17  10  2
 9.9988  10  6
 1.7  10  3
(2) 漸化式を逆向きに使って計算して、(1) での誤差の改善を確認せよ；
特に l = 9, 10 等大きい場合について
l=10 について見てみる
後進漸化式
前進漸化式
l=20 について見てみるとより違いは顕著
後進漸化式
前進漸化式
前進漸化式
前進漸化式での x~0 で大きな誤差
後進漸化式では x~0 でも十分 0 に近い値
後進漸化式
0.100000
0.200000
0.300000
0.400000
0.500000
0.600000
0.700000
0.800000
0.900000
1.000000
3.426744e+02
-1.371593e-01
3.908256e-04
-3.662896e-04
-2.153573e-05
3.815484e-06
-7.701172e-07
-3.546635e-07
2.769321e-08
4.001015e-08
7.271511e-21
7.441173e-18
4.286293e-16
7.599905e-15
7.064124e-14
4.363467e-13
2.032691e-12
7.701463e-12
2.491675e-11
7.116553e-11
2.000000
2.100000
2.200000
2.300000
6.834199e-08
1.101869e-07
1.738001e-07
2.684388e-07
6.825301e-08
1.101832e-07
1.738026e-07
2.684239e-07
23桁
10‐5
(2) の program 例
/* Bessel function */
/* up and down recursion relation */
#include <stdio.h>
#include <math.h>
int main(){
int l;
double x, dx, u, d, min, max;
FILE *outfile;
outfile = fopen("bessel2.dat", "w");
printf("Bessel order = ");
scanf("%d", &l);
#define start 50
double up(double x, int l){
int n;
double one, two, three;
min = 0.1;
max = 30.0;
dx = 0.1; one = sin(x)/x;
two = (sin(x)‐x*cos(x))/(x*x);
for(n=2; n<=l; n++){
three = (2.0*(n‐1.0)+1.0)/x * two ‐ one; /* 前進漸化式 */
for(x=min; x<=max; x+=dx){
u = up(x, l); /* 前進漸化式で計算 */
d = down(x, l); /* 後進漸化式で計算 */
fprintf(outfile, "%f %e %e¥n", x, u, d);
}
fclose(outfile);
one = two;
two = three;
}
return(three); /* 最終的に得られた j の値を返す */
}
double down(double x, int n){
int l;
double scale, j[start+2];
j[start+1] = 0.0;
j[start] = 1.0e‐3;
for(l=start; l>0; l‐‐){ /* j[start] から j[0] まで計算する */
j[l‐1] = ((2.0*l+1.0)/x)*j[l] ‐ j[l+1]; /* 後進漸化式 */
}
scale = ((sin(x))/x)/j[0]; /* 規格化のための補正係数 */
return (j[n]*scale); /* 補正した j[n] を返す */
}
}
数値微分
微分の定義：
df ( x)
f ( x  h)  f ( x )
 lim
h 0
dx
h
幅 h を小さくとると厳密な解に近づくが、数値的にこの通り微分
をすると、桁落ち誤差が発生し、正しくない答が得られる
手法１：前進差分
手法２：中心差分
後に微分方程式の数値計算を行う際に
この差分の計算を頻繁に利用する。
手法１：前進差分
最も直接的なやり方。
テーラー級数に展開する：
h2
h 3 ( 3)
f ( x  h)  f ( x)  hf ' ( x) 
f ' ' ( x) 
f ( x)  
2
6
f ( x  h)  f ( x ) h
h 2 ( 3)
f ' ( x) 
 f ' ' ( x) 
f ( x)  
h
2
6
数値計算での微分
f c ' ( x) 
f c ' ( x)
と置く。
h
f ( x  h)  f ( x )
 f ' ( x)  f ' ' ( x)  
h
2
前進差分
プログラムで：
diff = (f(x+h)‐f(x))/h;
近似誤差は ~ h の程度
・ x ～ x+h の区間における関数をその区間の両端の
２点を用いて直線で表す
・ h を小さくすると近似誤差を減らせるが、小さくし過
ぎると、桁落ちにより精度が落ちる。
x
x
h
2
xh
手法２：中心差分
前進差分：１ステップ h だけ前進する代わりに
中心差分： h/2 だけ前進し h/2 だけ後退する
f c ' ( x) 
f ( x  h / 2)  f ( x  h / 2)
 Dc f ( x, h)
h
f(x±h/2) をテーラー展開して、代入する：
h
h
(h / 2) 2
(h / 2)3 (3)
f ( x  )  f ( x)  f ' ( x) 
f ' ' ( x) 
f ( x)  
2
2
2
6
h
h
(h / 2) 2
(h / 2)3 (3)
f ( x  )  f ( x)  f ' ( x) 
f ' ' ( x) 
f ( x)  
2
2
2
6
h 2 ( 3)
Dc f ( x, h)  f ' ( x) 
f ( x)  
2
3!2
プログラムで：
近似誤差は
~ h2 の程度
中心差分
diff = (f(x+h*0.5)‐f(x‐h*0.5))/h;
前進差分より誤差が小さい
→ 前進差分ほど h を小さくしなくてもよい
（桁落ちを防げる）
x
x
h
2
xh
例題
関数
f ( x)  cos x
の微分を数値計算してプロットする
・刻みの数 n を適当にとって(n=100,1000, 等)、 x = [xa, xb] = [0, 2π] の各点
での微分値を求め、プロットするプログラムをつくる。微分計算で使う刻み幅
h は h = (xb‐xa)/n .
・前進差分、中心差分の両方を計算し、一行ごとに両結果を出力させる
x
f’(x): 前進差分
f’(x): 中心差分
0
‐0.00314158 ‐2.0793e‐15
0.00628319 ‐0.00942462 ‐0.00628313
0.0125664 ‐0.0157073 ‐0.012566
0.0188496 ‐0.0219893 ‐0.0188484
for (i=0; i<n; i++){
x = xa + i*h;
diff1 = ****
diff2 = *****
fprintf(***);
}
・・・
・・・
・・・
→ f’(x) = -sin x になっていることを確かめる
→ 前進差分と中心差分の結果がそれぞれどれくらいの誤差になっているかを
plot せよ。 ( 横軸を x, 縦軸を | 計算値 – (‐sin(x)) | でplot )
微分の計算値、その誤差各2種類を１つのファイルに書き出せばよい
fprintf(outfile, "%g %e %e %e %e ¥n", x, diff1, diff2, error1, error2);

Download Report