連続系アルゴリズム演習第3回 DE変換を用いた数値積分法

数値積分 積分を数値的に行う いくつかの点から元の関数を推察 サンプル点の値に対して、線形な変換 ニュートン・コーツ型 ガウス型 等間隔サンプル点 台形公式やシンプソン公式等 ガウス型 直行多項式の和で近似 サンプル点の値に対して、線形な変換

元の関数を等間隔点で数値積分する うまくいかない場合がある (当然だけれども)精度が悪い！ 格子点数 値 10 2.76098424523760455074 100 3.02067393649180804260 1000 3.10333743010756357705 10000 3.12949474161108387449 100000 3.13776694058186977898 1000000 3.14038285636361580444 10000000 3.14121008209783481036 中心部分(良い性質を持ってる部分)でも 差分化を細かくしているのが勿体ない 誤差は大体点の数に対して 1/sqrt(n)くらいで減る

変数変換 φ(t)をうまく選ぶことで、数値積分をやりやすくする！ 今回扱うDE変換はその一手法 無限区間積分に変換する

DE変換 元の関数 この部分がうまくいかなかった 変換後の関数 無限区間の積分がいるけど こっちはうまくいきそう！

DE変換の良さ 原点から遠ざかると、二重指数的に減少 程よく精度が出る 実質有限区間の積分でOK 差分の細かさによって、指数的に誤差減少 台形公式を使う 無限区間積分では、格子幅hに対してO(exp(1/h))の誤差 最適公式

双曲線関数 参考:

双曲線関数 指数的に 増加 指数的に ±1に漸近 指数的に 増加

有限区間から無限区間へのDE変換 有限区間[-1,1]の積分を、無限区間積分にする 変数変換をx=φ(t)として、 [-1,1]に収まらない場合は、収まるように事前に変数変換をす る 変数変換をx=φ(t)として、 tが-∞の時には、xは-1に近づき、tが+∞の時にはxは+1 に近づく

実際にやってみる 関数 を積分する 変数変換として、 変換後の関数は 指数で増大 exp(-exp(t))で1に近づく 関数 を積分する 変数変換として、 変換後の関数は 指数で増大 exp(-exp(t))で1に近づく exp(exp(t))^2で増大

そのまま計算する この例ではInfが出る(こともある) この例は既知の事実としてf(1)=Inf ある程度のtに対して 端点に特異点を持つ場合、特殊な処理が必要となる 特異点とは、Taylor展開不可能な点 この場合だと無限大に発散する点

解決案その1 tanh(π/2 sinh(t)) t=3あたりで丸め誤差の問題 1 1.0 - Δ t 1 t 1 t=3あたりでの丸め誤差の問題はない underflowはある(大体t=6あたり) t

解決案その1 問題はtanhの引数が1に近すぎるとき 関数もそれ用に変形して解いてみる tanhではなくて(1-tanh)という関数を求める この関数を積分すると、ましになる

解決案その2 もっと式変形をがんばる うまく項と項を削除する 具体的には、 　　を変形して解いていく

式変形続き

解決案1の結果 Δx = 0.1で、(-4.0, 4.0)を計算してみる 結果： 3.141592653589793116 ほとんど誤差がない 4*atan(1.0)=3.141592653589793116 場合によって丸め誤差程度は出る可能性がある これは避けられない

解決案2の結果 0.1刻みで[-4, 4]の範囲を計算 2.000000000000000000000000

DE変換の二つのポイント 刻み幅 打ち切り誤差 刻み幅をhとすると、誤差はO(e(-1/h))となる 二重指数的に減少する曲線を数値積分する 区間を固定して、点の数を増やすと、点の数に対して誤差が指数的 に減少する 打ち切り誤差 二重指数的に減少する曲線を数値積分する ものすごい速度で減るので、大体の大きさは見積もれる

刻み幅 台形公式の誤差 Δx→0で真の解に近づくとする 無限区間(-∞,∞)での積分誤差はO(exp(1/Δx)) 誤差 ブレがεよりも十分小さくなったら 誤差が十分小さいとみなす。 例えば、精度として1E-12程度の 誤差を許容するのであれば、 Δx=hの時の解と、Δx=h/2の時の 解の差が、1E-6よりも十分小さい ならば良い h 解の差が1E-6より十分小さい→Δx=hの時の誤差が1E-6より十分小さい →Δx=h/2の時の解が1E-12より小さいとみなせる

区間を固定した場合の誤差の推移 ちなみに、課題でも このグラフを提出すること 誤差 このあたりは丸め誤差レベル なので、これ以上は無理 点の数

打ち切り誤差 打ち切り誤差 無限区間積分を有限区間で切ることの問題 二重指数的に減少することを利用して見積もる 境界 拡大 Δに対して ほぼ0 Δ 境界での値Δがεより小さいならば、 誤差はhεより小さいと考えられる h

積分区間の決定 広くもなく狭くもない範囲で 大雑把に区間を広げていって、打ち切り誤差がなさそう な区間で止める 狭い場合 広い場合 打ち切り誤差が出る 台形公式の精度も出ない 広い場合 Infやら0やらNaNやらを扱う可能性 誤差減少の傾きが緩くなるかも 大雑把に区間を広げていって、打ち切り誤差がなさそう な区間で止める その区間で点を増やしていけばいい

ということで課題 以下の問題から2題以上選択する 問題1 問題2 必ず刻み幅と誤差の関係を示したグラフも出すこと このスライド中でも行った積分 二つの方法を両方実装すること 問題2 積分 問題1とまったく同じ方法で出来る 但し、式変形は自分でやること そんなに難しくはないはず 真の解はわからないので、ある程度細かく差分化したものを 真の解として、その差を使うことでグラフを作る

課題続き 問題3 問題4 積分 区間[0,∞)におけるDE変換は φ'(t)は自分で計算すること 区間[0,∞)において、元の関数が指数的に減少する場合は変 換として

課題続き 課題5 締切は、11月19日(月)とします 積分 区間(-∞,∞)におけるDE変換は φ'(t)は自分で(ry 本来の出題は11/5なので、そこから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 双曲線関数