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 双曲線関数