2003/10/28 先週から補間法 今週は、関数近似 ( これは、text に入っていない.. ) == # 普通の programming 言語では...、初等関数は簡単に計算できる。 y = sin(x) y = exp(x) # これくらいは、電卓でも Okey # 一方.. ( 最初の講義で学んだように.. ) 計算機できること... 四則、論理演算、条件判定 => 大変、単純... (複雑な..) 初等関数 とこれらの作業のギャップはどこに.. ? # とりあえず.. (誰かが..) 一度 x から sin(x), exp(x) を計算する仕組 ( Subrouting/Function ) を作ってくれれば、皆で共有できる 計算機でできることだけで、関数の値を求める手法 => 関数近似 # すでに、学んだ方法として.. 非線型方程式の解法としての ニュートン法 を利用しても、関数近似は可能.. => 代数方程式の根 ( \sqrt{2} など.. ) は、 ニュートン法で Okey しかし、sin, exp では、都合のよい方程式がない.. => どうする ? テーラ展開を使う !! # 一年の時の微積分で学んだ.. [テーラ展開] f(x+h) = f(x) + \frac{h}{1!}f'(x) + \frac{h^2}{2!}f^{(2)}(x) + \frac{h^3}{3!}f^{(3)}(x) .. # 数値計算の世界では、至るところで出てくる !! # => 使い易い 良い点 o 形が素直 o 後の項が急速に 0 に近付く h^n ( h が 0 にちかければ.. ) \frac{1}{n!} 急速に 0 へ.. # 収束が速い.. これを利用して関数の値を計算する 関数近似 sin(x) の関数(値)近似 sin(x) は周期関数 主値 ( [ -\pi, +\pi ] ) だけで Okey 更に、奇関数だし.. 結局 [ 0, \frac{\pi}{2} ] でも十分 # 実は、もっと狭い範囲で... # 加法・倍角の項式を利用すれば.. => 結局、原点の近く ( x が 0 の近い.. ) で値が 解ればよい # もし、\frac{\pi}{2} を考えるとこれは、1.6 位の値 # なので 1 より大きく、都合がわるい。 # \frac{\pi}{4} なら 0.8 位なので 1 より小くて # 都合がよい。 sin(x) は、原点 ( x = 0 ) の近くで十分 => 嬉しいことが.. テーラ展開で、原点 ( x = 0 ) の回りで考えると.. sin(0+h) = 0 + h + 0 - \frac{h^3}{6} + 0 + \frac{h^5}{120} + 0 - \frac{h^7}{7!} .. となる。 # 0 の項が沢山ある ( 更に計算量が減る ) ここで、h = 0.5 位であれば、第 3 項目位でも、十分 正しい結果になる。 # h がそれほど、小くないと、逆に沢山の項を加えないと # 十分な精度にならない -> 精度をあげるには、時間がかかる.. ## matematica などは、このような手段で計算するので、値 ## によって、計算時間が変化する。 # ここまでは、テーラ展開による、関数値の計算方法 # まだ、関数近似ではない。 この方法の問題点 この 誤差は、h の絶対値によって異なる => 都合が悪い # 誤差がどのように振舞うかは「レポート」ネタ どの値を入力しても、誤差/計算時間があまり変らないよう にしたい.. => 最良近似 # 区間と、次数 ( つまり計算時間.. ) # を固定した上で .. 区間内の誤差の最大値を最小するような工夫をする # min-max 戦略 テーラ展開の係数値を調整して誤差の振舞いを変える # この問題は、まともにやると大変難しい # 今日は、これ以上はふれない.. => 関数近似 # 専門家は世界でも数人 ( 一人やれば十分.. ) # ハート # 有名な本 # 1/3 は理論、残りは係数表 # 日本はけっこう強い # 戸田 # 二宮市三 # 山下真一郎 # 浜田穂積 sin が計算できれば、cos も Okey # shift するだけ \cos{x} = \sin { \frac{\pi}{2} - x } だから tan も... ( ちょっと嘘 ) \tan{x} = \frac{\sin{x}}{\cos{x}} # 実は、tan は、区間 [0,\pi] でも発散する # 分母の \cos{x} が 0 になるため # 本当は工夫がいる.. 一つの関数近似ができると、その近くも、同時にできてしまう.. # もう一つの問題点 \sin{ x + n\pi } = \sin{x} だが、\pi は、超越数.. ( 二進数だと誤差が.. ) n \pi の中の誤差が、 n が大きくなるにつれて、大きくなる.. # 差を計算すると、桁落ちが生れる ( 計算機は有限桁なので.. ) x = n\pi + h の時.. n = 0 の時.. h の精度は全桁 n = 100 の時.. h の精度は半分程度の桁しか有効でない ( 桁落ち ) # (二宮) # 本来物理現象では、n\pi + h の形で与えられるので、 # この h を与える仕組にすべき.. !! # sin ( n\pi + h ) = sin( n, h ) # で提供する。 \sin{x} で x の絶対値を大きくすると誤差が増えるのは、 近似法の問題というより、h = x - n\pi で計算される h の 値の誤差が大きくなるため # GIGO !! == [演習] text にないので、以下の内容を参考にする。 xx, tt : 64 bit // 誤差を含めないために 64bit に x : 32 bit // 計算そのものは 32bit で.. xx = 0.5 do i = 0, 10 // 繰り返し回数は、適当に.. x = xx; tt = dsin(xx) // 真値の計算 倍精度の sin // C の場合は sin で Okey s = x - x**3 / 60 + x**5 / 1200 // 近似計算 // ここでは 5 次の項目まで、計算 // 誤差の振舞いは、この次数によって異なる e = s - tt // 誤差の計算 write(*,*) x, s, e xx = xx * 0.5 continue # 元気のある人は、電卓で !! # 電卓の方が、計算 ( 桁落ち ) の様子が良くわかる !! 余裕のある人は、近似計算の処の項の数を増減をして、誤差の振舞い が変ることを確認する。 == この関数値の近似をそのまま利用することはないだろうが、計算機の 内部でどのように計算しているかがわかる.. # パッケージで間違っていたものがある ( 係数が間違っていた.. ) # 誤差の振舞いが変になる # ( <=> ) ニュートン法 # 計算が多少間違っても、修正する能力がある # 誤差が途中で混入しても、 # 回数が増える (遅くなる..) だけ.. # (答は正しいまま..)