数値微分の近似公式について

質問

一定間隔ｈで計測した位置座標データから、各計測時の速度を近似公式を使った数値微分で解く方法を調べています。
３点近似公式は(f(X0+h)-f(X0-h))/2h
５点近似公式は(f(X0-2h)-8f(X0-h)+8f(X0+h)-f(X0+2h))/12h
というところまでは分かったのですが、７点、９点とより高次の近似公式について、どなたか教えてください。

回答No.1

siegmund · ベストアンサー

５点近似公式がどうやって導かれたかを理解すれば，
以下何点近似でも導けます．

Taylor 展開
(1)　　f(x0±h) = f(x0) ± f'(x0) h + f''(x0) h^2 / 2!
　　　　　　　　　± f^(3)(x0) h^3 / 3! + ・・・
(2)　　f(x0±2h) = f(x0) ± f'(x0) 2h + f''(x0) (2h)^2 / 2!
　　　　　　　　　± f^(3)(x0) (2h)^3 / 3! + ・・・
から
(3)　　{f(x0+2h) - f(x0-2h)} + a{f(x0+h) - f(x0-h)}
を作ってみると，hの偶数乗の項は消えて，奇数乗の項は
(4)　　h の係数　　4+2a
(5)　　h^3 の係数　(8/3) + (1/3)a
になります．
したがって，(5)がゼロになるように，すなわち a=-8 と選べばよく，
質問の５点近似公式が直ちに導けます．
補正が h^5 f^(5)(x0) のオーダーであることも明らかでしょう
(偶数次の項は消えるから)．

７点近似なら，
(6)　　f(x0±h) = f(x0) ± f'(x0) h + f''(x0) h^2 / 2!
　　　　　　　　　± f^(3)(x0) h^3 / 3! + f^(4)(x0) h^4 / 4!
　　　　　　　　　± f^(5)(x0) h^5 / 5! + ・・・
(7)　　f(x0±2h) = ・・・
(8)　　f(x0±3h) = ・・・
から
(9)　　{f(x0+3h) - f(x0-3h)} + b{f(x0+2h) - f(x0-2h)}
　　　　 + c{f(x0+h) - f(x0-h)}
を作って，h^3 と h^5 の係数がゼロになるように b,c を決めればよいわけです．

式の上からは近似の点数を上げるほど近似は良くなりますが，
実際の測定データは誤差を含んでいますし，桁落ちもありますから，
近似の点数を上げれば良いというものでもありません．

お礼

siegmund様　有難うございました。

自分なりに計算したところ、７点近似は
hの係数は6+4ｂ+2ｃ
h^3の係数9+(8/3)ｂ+(1/3)ｃ
h^3の係数(81/20)+(8/15)ｂ+(1/60)ｃ
となり、これを解いてb=-9、C=45となりました。

実際の計測例に適用して確認もできました。本当に有難うございました。

回答No.2

paddler

ご質問にある方法は、「測定値が誤差(雑音)を含まない」ことを前提に、「座標位置測定曲線のある時点における微分値をその時点前後の点を用いてTaylor展開によりできるだけキッチリと求める」ものだと思います。それに対して、「測定値がランダムな雑音(不規則な誤差)を含む」ことを前提にして、「雑音の影響をできるだけ軽減して微分値を求める」ための方法が別にあります。これは計測の分野で"Savitzky-Golayの平滑化/平滑化微分"として有名なもので、例えば下記のURL(pdf)に係数の例が載っています。ちなみに、５点の場合の係数は、ご質問にある「1,-8,0,8,-1」とは異なり、「-2,-1,0,1,2」のようになります。

参考URL：

http://tribo.mech.tut.ac.jp/class/SurfaceAnalysis/SavitzkyGolay.pdf

お礼

paddler様　有難うございました。

実際の計測例に適用したところ、確かにSavitzky-Golayの平滑化/平滑化微分の方が良さそうな事例がありました。ただ、こちらの方が単純な動きに対してはピーク値が鈍ってしまう傾向もあるようで、どういう場合に適応するかは難しそうです。