2005/11/08 数値解析学と演習 福井 先生 数値微分と数値積分 <=> 数式処理による微分と積分 微分 -> 簡単 エントロピーが増大 積分 -> 難しい エントロピーが現象 数式微積分ができれば、数値的な手法は不要なようだが.. 数式処理では解けない場合が沢山ある 解けた場合でも効率が良い場合がある [数式処理では解けない 例] 関数が、測定値として与えられている コピー機のエッジチェック 数値微分している エッジ内部は、平坦化する処理が入っている 写真モードの場合は、データに忠実 自動車のギアボックスの雑音 雑音を平滑化するために、変化のおおきい所をチェック => (論理的にやってもよかったが) 数値的にやった 多変数関数 理論的にもできないわけではないが.. 極値が欲しい場合は、数値的にやった方が簡単 == 数値微分 微分の定義 \frac{df(x)}{dx} = \lim_{\delta x \rightarrow \infty} \frac{f(x+\delta x) - f(x)}{\delta x} # 計算機では、「有限で諦める」必要があるので近似 '=, \frac{f(x+\delta x) - f(x)}{\delta x} # 前進差分 \delta x が十分に小さい値 # 数学で「十分に」とあったら、ようするに「解らない」ってことだ.. # 一般には、「十分に」は「解らない」だが、関数の性質が解っていれば、「解る」 ## 数値計算の一番良い方法は、答がわかっている場合に作れる ## => 矛盾している.. # 連続と、離散は異る ( 風呂の水と洗面器で汲むか、ざるで汲むかの違いがある...) # 連続の場合は、完備性がある # [例] 音声のサンプリング # y = \sin{x} なのに、サンプリングを π周期で取ると結果的に y = 1 になる # => サンプリングが悪い # y が \sin{x} だとわかっていれば、π周期にしなかった # サンプリングを巧くするには、ある程度、関数の性質を知る必要がある 上記の近似計算は、一次 (直線) の近似 関数が連続であることがわかっていれば、 \frac{df(x)}{dx} '=, \frac{f(x) - f(x-\delta x)}{\delta x} # 後退差分 も使える。これを利用すれば、二次の近似式も作れる。 [関数の微分の一次近似の誤差評価] # テーラ展開を利用して.. f(x+Δx) = f(x) + \frac{Δx}{1!)f'(x) + \frac{Δx^2}{2!)f''(x) .. これから f'(x) = \frac{f(x+Δx) - f'(x)}{Δx} + \frac{Δx}{2!)f''(x) .. 上記の近似は、後の項を捨てている ( 打ち切り誤差 ) つまり、 \frac{Δx}{2!)f''(x) .. が誤差項 # f''(x) が零 ( つまり変曲点 ) では、f'''(x) の項目が誤差になる つまり、誤差は、f''(x) が理解れば、解るわけだが... # f'(x) が知りたいに、f''(x) が解るわけがない !! ## 数学的なアプローチは、いつもこんなかんじ ## 一般には、f''(x) の上界などを調べる 一般に、Δx を小くすれば、打ち切り誤差は 0 に近付く => 無限桁の計算であれば、これで Okey [数値微分に於ける、数値的な振舞い] f(x+Δx) - f(x) の計算を考える この二つの数は値が近い 差を取ると、有効桁数が減ってしまう ( 桁落問題 ) Δx が大きい f(x+Δx)= 1.345678 <= 8 桁の有効数字 -) f(x) = 1.234567 ------------------------- 0.111111 <= 7 桁の有効数字 Δx が小さい f(x+Δx)= 1.234568 <= 8 桁の有効数字 -) f(x) = 1.234567 ------------------------- 0.000001 <= 1 桁の有効数字 # 世に言う「丸め誤差」 ## これは、「情報落」であって、「丸め誤差」とはいいたくない 数値微分 Δxを小くする作業 結果的に f(x+Δx) と f(x) の値が近くなる 誤差の振舞い 打ち切り誤差 右上り 丸め誤差 右下がり 数値積分の誤差 V 字 # Δx を大きくしても小くしても、誤差が大きくなる # 二つのグラフ ( X 型 ) の上の部分 ## 誤差がキャンセルして、V の先が X の交点の下に行くことがある Δx をどのようにすれば、良いか ? # 永坂先生と福井先生の結果 d(E_t + E_R)/dx = 0 として考える Δx_{opt} は、計算機の精度の半分の精度 ただし、f''(x) がほぼ 1 の時 # E_t : 打ち切り誤差 と E_R : 丸め誤差の値が同じ、程度の精度になる。 単精度の場合 24bit の精度があるので、12 bit 位の精度をもつΔx がよい。 例 x = 4 => Δx = 4/2^12 = 4/4096 = 0.001 x = 4000 => Δx = 1 # Δx を無闇に小くしても精度があがるわけではない。 後退差分の誤差評価 Δx を - Δx にするだけ 誤差項の符号が変る [二次の近似 ( 二次関数による近似 )] # やはり、テーラ展開を用いて.. f(x+Δx) = f(x) + \frac{Δx}{1!)f'(x) + \frac{Δx^2}{2!)f''(x) .. f(x-Δx) = f(x) - \frac{Δx}{1!)f'(x) + \frac{Δx^2}{2!)f''(x) .. 上から、下を引いてあげると、 f(x+Δx)-f(x-Δx) = 2(Δx f'(x) + \frac{Δx^3}{3!)f'''(x) ..) # 以下、奇数の次数の項目だけがでる これより、 f'(x) = \frac{f(x+Δx)-f(x-Δx)}{2Δx) - \frac{Δx^2}{3!)f'''(x) ..) # 一般に、このようにテーラ展開を繰り返すと、何次元の公式でも作れる [二次の近似の計算量] 一次の近似では、関数の値を二度利用している 二次の近似でも、関数の値を二度しか計算していない 対称性があるので 数値計算では対称性があると、(同じ計算量に対して)精度がよくなる # 「悪さ」をする場合もある # 実際の計算では.. # Δx を変えながら、連続的に計算すると.. # 一次の場合 => 始点は同じ => 終点が違うだけ # 二次の場合 => 前後が計算しなおし => 両方計算するので計算量が増える 公式の次数が増えると、E_t は、エッジがきつくなる ( Δx が大きい方にずれる ) Δx は大きくする必要がある 二次の場合は、Δx_{opt} は精度の 1/3 程度大きさ 計算機の精度が増えると E_R のエッジがきつくなる ( Δx が小さい方にずれる ) Δx は小くする必要がある 倍精度の場合は、Δx_{opt} は精度が大きくなる == [数値微分の演習問題] ID : 11/08 締切 : 11/21 課題: 前進差分で、関数の微分を近似する f'(x) = \frac{f(x+Δx)-f(x)}{Δx} 関数: y = e^x # 微分しても結果が同じなので、答が楽 値: x = 1 # exp の場合は、どこでもオーバフローしない限り同じ Δx = 1, 0.5, 0.05, ..., 2^{-24} 結果: 誤差の結果 ( 誤差のグラフが出るとよい ) # まるめ誤差は、コーディングによって異るので、グラフも人それぞれになる == E_t : 打ち切り誤差 ( トランケーション誤差 ) 〜 Δx^k ( k は近似の次数 ) 式の近似によって生れる E_R : 丸め誤差 ( ラウンド誤差 ) 〜 \frac{Δx} 有効桁数によって生れる == 数値積分の応用例 有限要素法 多重積分 測量 積分 不定積分 定積分 <= 数値積分ができるのはこちらのみ \int_{a}^{b} f(x) dx = F(b) - F(a) 数値積分 => 発想は、「面積を求める」ということ ( その意味では簡単 ) あとは、近似をどのように精度を上て行うか ? 数値積分の方法 ニュートン = コースの積分公式 等間隔 # 積分区間を等間隔に分割してその面積の総和を計算する 分割した領域の面積は、台形として面積を求める => 台形公式 # 意外と、これで、答えが出てしまう !! ## 個々の台形に入る誤差は、関数の値の誤差 ## 連続した台形の前後では、関数の値の誤差が逆符号になる ## => 誤差同士がキャンセルしてしまう ## => 結局は、全体の区間の両端での誤差だけが入る ガウス (= ルジャンドル)の積分公式 不等間隔 # 区間が無限大でも扱える [台形公式(一次近似)] S = \int_{a}^{b} f(x) dx の時、区間 [a,b] をΔx の幅で、分割すると.. S '=, Δx(f(a) + f(a+Δx))/2 + Δx(f(a+Δx) + f(a+2Δx))/2 + .. = Δx(f(a)/2 + f(a+Δx) + f(a+2Δx) + .. + f(b-Δx) + f(b)/2) となる。 # [余談] # f(a) と f(b) が相手がいない # 相手ができる 「f(a)=f(b)[周期関数]」ならどうだろうか ? # 端点の誤差もキャンセルする可能性がある # 周期関数をその周期の区間で、台形公式を使うと凄く精度がでる ## f(a)=f(b) となる条件はむずかしいが.. ## f(a)=f(b) = 0 ( 無限区間で、両側が 0 になる場合.. ) も、一種の周期関数 ## 更に、関数の形を変数変換して、周期関数に変形して関数する ## 森、高橋の理論 ( 二重指数関数法 ) ## 複素関数論でのローラン展開すると良くわかる。 # 公式の詳細は text (p. 405) 参照 区間の幅(Δx)を半分にした値を利用(補外)すると、高次元の公式が得られる => ロンバーグ積分法 # 精度が、「数値的」にでる 誤差 2次 4 次 8 次 Δx 1 F_0^0 1/2 F_1^0 F_1^1 1/4 F_2^0 F_2^1 F_2^2 2^{-n} F_n 一般に、台形公式の誤差は、区間 Δx の二乗の大きさに比例する。 # 本来の次数は 1 だが、積分するので二乗になった 台形公式の誤差の主要項目は、 Δx^2 f''(x) の形 Δx^2 f''(x) - 4 ((Δx/2)^) f''(x) の形で、f''(x) の誤差項目を消すことができる。 => 補外 「補外」を行うことにより、(台形方法を使うだけで..)高次の公式と同じ精度が出る [シンプソン公式(二次近似)] 台形法では、2 点を利用した ( 直線で近似 ) /| / | | | 1 1 1 1 1 1 1 1 .. 1 1 1 1 1 2 2 .. 2 1 シンプソン公式では、3 点を利用 ( 放物線で近似 ) /| / | /| | / | | | | | 1 4 1 1 4 1 1 4 1 1 4 1 .. 1 2 2 .. 2 1 # 積分方法では、二つの台形を集めて次数をあげる。感覚としては、一つの台形を分割する方が自然だが ## 高次元の公式では、台形公式とちがって、(単に)周期だからといって精度は上らない ### これは、関数の値だけでなく、微係数も異るから.. [ニュートン = コースの積分公式の一般形(text p.457)] m C_0 C_1 C_2 C_3 C_4 台形 1 1/2 1/2 シンプソン 2 1/3 4/3 1/3 3 3/8 9/8 9/8 3/8 4 14/45 64/45 24/45 64/45 14/45 # 法則性がある # 係数の和が、m (次数)の値に一致する # 係数の和が区間のサイズになる # => 次数をふやすと区間が増える # テキストでは、M の大きいものがかいてあるが.. ## 数値計算的には、「M = 7 までしか役にたたない」というのが常識 ## 数値計算では、負の係数が出ることを嫌がる ## 負の計算があると誤差が累積するので... 計算量的には、同じ次数にすると、ロンバーグ積分の方が有利 # 収束判定が必要で、そこに少し難しい問題がある。 ニュートン = コース系では、関数の値への重み付けによって、公式が定まる 次数が変われば係数が変るだけ == 直交多項式 P_n(x) : 次数 n の多項式とする 区間 [a,b] の積分で、 \int_{a}^{b} W(x)P_n(x)P_j dx = 0 ( j = 0, .., n-1 ) # W(x) は重み関数 を満す、P_n(x) を「直交多項式系列」と呼ぶ 直交系は、無駄な情報を持たない # cf. 固有ベクトル , 直交座標系, 内積が 0 なベクトル ルジャンドル多項式 ( Text 参照 ) 直交多項式系の典型的な例 P_0(x) = 1 P_1(x) = x P_2(x) = (3x^2-1)/2 P_3(x) = (5x^3-3)/2 これらの関数の 0 点での関数の値を利用して積分を行う => ガウス=ルジャンドル法 係数に平方根などがでて複雑になるが.. 関数の値の量が少くて、精度の良い積分が得られる # いまは、係数を手で入力することがない # ライブラリィでは、係数をプログラムに埋め込みしている # ## 実際には、係数を計算するプログラムがあり、それで生成した係数を、計算式に埋め込んでいる。 ## 「何故、直交多項式を使うと精度が上るか」は時間があったときにでも.. == ニュートン = コース系 区間を等間隔に分けて、多項式で近似 ガウス = ルジャンドル系 直交多項式の零点での関数を計算 一つ一つの関数の計算は大変だが、計算回数は減る # 関数の計算が大変な時に利用する。 # 元の関数を、「都合のよい関数系を作って、それで表現する」ことにより、必要な計算を行う。 ## 係数は、人間の手を介さないようにしないと間違える !! == [数値積分の演習問題] ID : 11/15 締切 : 11/21 課題: 台形法で、積分 \int_{a}^{b}f(x)dx = \frac{f(a)}{2}+\sum_{i=1}^{n-1}f(a+\frac{i(b-a)}{n}) + \frac{f(b)}{2} 積分: \int_{0}^{1}\frac{4}{1+x^2}dx 値: π 結果: 区間の分割数をふやすと、誤差がどうなるかを調べる [おまけ] モンテカルロ法 x, y の値 ( [0,1] 区間 ) の乱数を N 回数発生し、 x^2 + y^2 < 1 となる個数 n を数える n/N が、π/4 に近付く => 計算量に比して精度が殆どでない 多重積分などのように、台形積分が大変な場合に利用する == まとめ 積分公式は、ニュートン = コース系の 台形公式 と シンプソン がメイン。 これの次元を上げるならば、 ロンバーグ積分 その他に、 ガウス系 もある。 数値積分は、有限要素法のベース関数を作るのに利用する位