Transcript 授業スライド
連続系アルゴリズム演習第3回
DE変換を用いた数値積分法
数値積分
積分を数値的に行う
いくつかの点から元の関数を推察
ニュートン・コーツ型
ガウス型
等間隔サンプル点
台形公式やシンプソン公式等
直行多項式の和で近似
サンプル点の値に対して、線形な変換
元の関数を等間隔点で数値積分する
うまくいかない場合がある
(当然だけれども)精度が悪い!
1
dx
1
1 x
2
格子点数
値
10
2.76098424523760455074
100
3.02067393649180804260
1000
3.10333743010756357705
10000
3.12949474161108387449
100000
3.13776694058186977898
1000000
3.14038285636361580444
10000000
3.14121008209783481036
中心部分(良い性質を持ってる部分)でも
差分化を細かくしているのが勿体ない
誤差は大体点の数に対して
1/sqrt(n)くらいで減る
変数変換
b
a
1 ( b )
f ( x)dx 1
(a)
f ( (t )) ' (t )dt
φ(t)をうまく選ぶことで、数値積分をやりやすくする!
今回扱うDE変換はその一手法
無限区間積分に変換する
DE変換
元の関数
この部分がうまくいかなかった
変換後の関数
無限区間の積分がいるけど
こっちはうまくいきそう!
DE変換の良さ
原点から遠ざかると、二重指数的に減少
実質有限区間の積分でOK
程よく精度が出る
差分の細かさによって、指数的に誤差減少
台形公式を使う
無限区間積分では、格子幅hに対してO(exp(1/h))の誤差
最適公式
双曲線関数
e x e x
x3
sinh x
x
2
3!
e x ex
x2
cosh x
1
2
2!
sinh x
tanhx
cosh x
e e
4
sinh2 x cosh2 x 1
sinh2 x cosh2 x
2x
2 x
x5
5!
x4
4!
(sinh x)' cosh x
(cosh x)' sinh x
sinh2 x cosh2 x
1
(tanhx)'
cosh2 x
cosh2 x
e ix e ix
x3 x5
x
参考: sin x
2
3! 5!
e ix e ix
x2 x4
cos x
1
2
2! 4!
sin x
tan x
cos x
双曲線関数
指数的に
増加
指数的に
±1に漸近
指数的に
増加
有限区間から無限区間へのDE変換
有限区間[-1,1]の積分を、無限区間積分にする
[-1,1]に収まらない場合は、収まるように事前に変数変換をす
る
変数変換をx=φ(t)として、
(t ) tanh( sinh t )
2
tが-∞の時には、xは-1に近づき、tが+∞の時にはxは+1
に近づく
実際にやってみる
関数
1
1
1
1 x
を積分する
dx
2
変数変換として、
(t ) t anh( sinh t )
2
' (t )
2
cosh t
cosh2 ( sinh t )
2
変換後の関数は
f ( x)
1
1 x2
cosh(t )
指数で増大
2
f (tanh(sinh( t )))
dt
2
2
cosh ( sinh(t ))
2
exp(-exp(t))で1に近づく
exp(exp(t))^2で増大
そのまま計算する
この例ではInfが出る(こともある)
この例は既知の事実としてf(1)=Inf
f (tanh(sinh ( t )))
2
ある程度のtに対して
1 tanh(sinh(
t )) 1 M
2
具体的には、t=3.1あたり
端点に特異点を持つ場合、特殊な処理が必要となる
特異点とは、Taylor展開不可能な点
この場合だと無限大に発散する点
解決案その1
tanh(π/2 sinh(t))
t=3あたりで丸め誤差の問題
1.0 - Δ
1
0
1 tanh( sinh(t )) 1
2
t
e2
e
2
sinh( t )
sinh( t )
e
e
sinh( t )
2
sinh( t )
2
2e
e
2
1
0
sinh( t )
2
sinh( t )
e
sinh( t )
2
t=3あたりでの丸め誤差の問題はない
underflowはある(大体t=6あたり)
t
解決案その1
問題はtanhの引数が1に近すぎるとき
tanhではなくて(1-tanh)という関数を求める
関数もそれ用に変形して解いてみる
sinh( t )
2
2e
, y 1 x, y 1 tanh( sinh(t ))
sinh( t )
sinh( t )
2
1 x2
e2
e 2
1
1
g ( y)
1 (1 y ) 2
2 y y2
f ( x)
1
0
f ( y)
2
cosh(t )
cosh2 ( sinh(t ))
2
dt g ( y )
0
2
cosh(t )
dt
cosh2 ( sinh(t ))
2
この関数を積分すると、ましになる
解決案その2
もっと式変形をがんばる
うまく項と項を削除する
具体的には、
t t
2 sinh( t ) 2 sinh( t )
(e e )
(
e
e
)
2
4
f (tanh( sinh t ))
f
sinh( t )
sinh( t )
sinh( t )
sinh( t )
2
2
cosh ( sinh t )
(e 2
e 2
) (e 2
e 2
)2 / 4
2
cosh t
を変形して解いていく
式変形続き
t t
2 sinh( t ) 2 sinh( t )
(e e )
(e
e
)
4
f
sinh( t )
sinh( t )
sinh( t )
sinh( t )
2
(e 2
(e 2
2
2
e
)
e
)
/4
1
(e t e t )
sinh( t )
sinh( t )
sinh( t )
sinh( t )
2
e 2
)2
(e 2
e 2
) 2 (e
1
sinh( t )
sinh( t )
(e 2
e 2
)2
e2
sinh( t )
2
e
2
sinh( t )
2
(e t e t )
(e 2
e t e t
e2
sinh( t )
e
sinh( t )
2
sinh( t )
e
sinh( t )
2
2
)
解決案1の結果
Δx = 0.1で、(-4.0, 4.0)を計算してみる
結果: 3.141592653589793116
ほとんど誤差がない
4*atan(1.0)=3.141592653589793116
場合によって丸め誤差程度は出る可能性がある
これは避けられない
解決案2の結果
2
et e t
e
2
sinh( t )
e
sinh( t )
2
dt
2
et e t
e
2
sinh( t )
0.1刻みで[-4, 4]の範囲を計算
2.000000000000000000000000
e
sinh( t )
2
dt
DE変換の二つのポイント
刻み幅
刻み幅をhとすると、誤差はO(e(-1/h))となる
区間を固定して、点の数を増やすと、点の数に対して誤差が指数的
に減少する
打ち切り誤差
二重指数的に減少する曲線を数値積分する
ものすごい速度で減るので、大体の大きさは見積もれる
刻み幅
台形公式の誤差
Δx→0で真の解に近づくとする
無限区間(-∞,∞)での積分誤差はO(exp(1/Δx))
誤差
ブレがεよりも十分小さくなったら
誤差が十分小さいとみなす。
h
例えば、精度として1E-12程度の
誤差を許容するのであれば、
Δx=hの時の解と、Δx=h/2の時の
解の差が、1E-6よりも十分小さい
ならば良い
解の差が1E-6より十分小さい→Δx=hの時の誤差が1E-6より十分小さい
→Δx=h/2の時の解が1E-12より小さいとみなせる
区間を固定した場合の誤差の推移
ちなみに、課題でも
このグラフを提出すること
誤差
このあたりは丸め誤差レベル
なので、これ以上は無理
点の数
打ち切り誤差
打ち切り誤差
無限区間積分を有限区間で切ることの問題
二重指数的に減少することを利用して見積もる
境界
拡大
Δ
境界での値Δがεより小さいならば、
誤差はhεより小さいと考えられる
h
Δに対して
ほぼ0
積分区間の決定
広くもなく狭くもない範囲で
狭い場合
広い場合
打ち切り誤差が出る
台形公式の精度も出ない
Infやら0やらNaNやらを扱う可能性
誤差減少の傾きが緩くなるかも
大雑把に区間を広げていって、打ち切り誤差がなさそう
な区間で止める
その区間で点を増やしていけばいい
ということで課題
以下の問題から2題以上選択する
問題1
必ず刻み幅と誤差の関係を示したグラフも出すこと
1
このスライド中でも行った積分 1
二つの方法を両方実装すること
dx
1 x2
問題2
1
1
dx
積分 (2 x)(1 x) (1 x) 1.94
問題1とまったく同じ方法で出来る
1/ 4
3/ 4
但し、式変形は自分でやること
そんなに難しくはないはず
真の解はわからないので、ある程度細かく差分化したものを
真の解として、その差を使うことでグラフを作る
課題続き
問題3
dx
積分 (1 x) 1
区間[0,∞)におけるDE変換は
φ'(t)は自分で計算すること
0
2
(t ) exp(
2
sinh( t ))
問題4
積分 x exp( x)dx 1
区間[0,∞)において、元の関数が指数的に減少する場合は変
換として
(t ) exp(t exp(t ))
φ'(t)は自分で計算すること
0
課題続き
課題5
dx
積分 1 x
区間(-∞,∞)におけるDE変換は
φ'(t)は自分で(ry
2
(t ) sinh( sinh( t ))
締切は、11月19日(月)とします
本来の出題は11/5なので、そこから2週間
2
References
[1]数値解析 第2版, 森正武, 共立出版株式会社, 2002年
[2]数値計算法の数理, 杉原正顯 室田一雄, 岩波書店, 1994
年
[3]FORTRAN 77 数値計算プログラミング 増補版, 森正武, 岩
波書店, 1987年
[4]Quadrature formulas obtained by variable transformation
and the DE-rule, Masatake MORI, Journal of Computational
and Applied Mathematics 12&13(1985) 119-130
[5]The double-exponential transformation in numerical
analysis, Masatake Mori Masaaki Sugihara, Journal of
Computational and Applied Mathematics 127(2001) 287-296
[6]Wikipedia 双曲線関数