2004/11/16 解析学とコンピュータ 補間法の応用 前回は、数値微分だった。今回は、数値積分。 == 数学的には、微分するより、積分は大変 数学に於ては、一般に、エントロピーを増大させるのは、簡単だが、 エントロピーを減少させるのは大変 例 : 数式の展開 ( エントロピー増大 ) より、因数分解( エントロピー減少 )の方が大変 数値計算では、微分より積分の方が簡単。 例 : ランプ関数 ( Step Function ) を微分しようとする発散。積分は容易 数学の積分 不定積分を扱う 計算機の積分 (数値積分) 定積分を扱う ( 不定積分は苦手 ) 関数 f(x) が与えらたときに、f(x) を a から b まで、定積分するには ? x 軸 ( y = 0 ) と、 x = a, y = b, そして f(x) で囲まれた面積を計算する 面積をどう計算するか ? 縦に細かくきって短冊を作り、その短冊 (羊羹) の和を計算する 短冊の上の部分 ( f(x) により決る ) をどのように近似するか ? 直線 ( 一次式 ) で近似する : 台形法 放物線 ( 二次式 ) で近似する : シンプソン法 理論的には.. 短冊の幅を小くする ( 打ち切り誤差を下げる ) 短冊の上辺の近似の次数を挙げる # 微分の場合は、h を小くすると、計算誤差 ( 丸め ? ) が生じた # 積分の場合は、それほど多くの誤差が生じない ( Why ? : 加算だから.. ) # ただ、もちろん、計算量は増える。 近似関数の次数 d は、h に対する誤差の傾きになる 実数が大きければ、同じ h でも、精度が良くなる Text (p.112) ニュートン=コーツの積分方式 n = 1 の時に台形法 n = 2 の時にシンプソン法 n = 3 ... # これまでの方法 (台形法、シンプソン法) を、一つの形にまとめたのがニュートン=コーツの積分方式 text の表の意味 計算を行う時の各関数値の係数 ( 係数の和は常に n ) 例えば、ラグランジェの補間を使えば解る。 # 係数が n になる理由 # 定値関数の積分を考えればよい # 関数が定値なので、その時に、f(a) (b-a) にするには、少くても、係数の和が n ( 次元 ) でなければならない。 # # この性質 ( 係数の和が n ) は公式のチェックにも利用できる == 数値積分の誤差評価 cf. text p.112 結果の一番大きな誤差で評価する。 誤差は、h^3/12 になる 一般的には、関数が判らないと、より細かい誤差は判らない。 短冊の加算をすると誤差がキャンセルされることもある。 特に周期関数で、台形方式だと、誤差がキャンセルされて誤差が出ない 急速に零に収束する関数は、周期関数とみなせる ( 両側で零 ) 誤差無しで計算できる 元の関数を、変数変換して、上の形に変換 二重指数関数変換 ( 森毅先生 ) 複素関数論を利用して、零点のローラント変換を利用して証明できる # 電気屋は、複素関数論が必須 # on/off <-> 実数/虚数 # 交流 <-> e^{x+y i} がニュートン=コーツの積分方式は、n = 1, 2, .. と沢山ある 確かに、n が大きいと、誤差カーブが急になり、理論上の誤差は小くなる しかし、数値計算故は、 n = 1, 2 で十分 == ロンバーグ積分 (p.114) はじめは、台形法で区間全体を h ( = b-a )として値を求める ( これは誤差が多い ) 誤差は、M h^3/12 h を更に半分にすると 誤差は、M (h/2)^3/12 x 2 = 1/4 M h^3/12 つまり、誤差の関係が解るので、最初の値の誤差を後の誤差の値で、キャンセルさせる 同様に h/4 だとすると更に、誤差に関する情報がわかるので、 0 1 2 0 R_00 \ 1 R_10 - R_11 \ \ 2 R_20 - R_21 - R_22 \ \ \ 3 R_30 - R_31 - R_32- R_33 台形方法だけで誤差をキャンセルすると、次数が一つ高い結果が得られる その方法なんどか繰り返すと、より次数の高い計算結果が得られる # 段数は、数段で、十分に精度が出る # この手法は他にも利用されている ( 展開の加速 ) : 加速法 数値積分への応用の場合は、途中の点の計算も必要なので、計算量が増える (2^n) 数値微分への応用の場合は、途中の点の計算がそれほど不要 (2n) # 積分と微分の違い == ここまでは、積分区間が有限の場合の例だった。一般には、無限区間の積分もありうる これは、これまでの方法では積分できない 一般の場合は解けない ( なぜなら収束しないから.. ) 特別な場合 ( 収束性が明らかな例 ) はできる \int_0^{\infity}e^{-x}g(x) dx ラゲール=ガウスの公式 \int_{-\infity}^{\infity}e^{-x^2}g(x) dx エルミート=ガウスの公式 # 今日は、時間がないので、詳しい話はなし.. == 数値積分は、それ単独で、積分値そのものを求めることはあまりない ( 既に、ライブラリになっているので ) 有限要素法の要素値を数値積分法で求める形で利用することが多い == [演習] 11/16 台形法で次の関数を数値積分 \int_0^1\frac{4}{1+x^2}dx # 答は π (円周率) になるはず 計算する時には、f の値を加えてから h と 1/2 を掛けるようにする # この計算は電卓でも可能 区間を 4,8,16 分割をして、誤差が少くなることも確認する