二階微分の数値計算法について

質問

関数(sinとかcosとか)それ自身は、Excelの関数に、必要に応じ四則演算や合成を用いて定義されているものとします。
この関数を、一階微分しようとした場合、たとえば

dt=0.00001
fplus=f(t+dt)
ft=f(t)

dfdt=(fplus-ft)/dt

のように計算すれば、そこそこの精度で計算結果がえられ、爆発のような不審な挙動は
なさそうな感じです。

一方、二階微分は数式で書くと

f''(t)=(f'(t+dt)-f'(t))/(dt)

で、

f'(t)=(f(t+dt)-f(t))/dt
f'(t+dt)=(f((t+dt)+dt)-f(t+dt))/dt

なので、結局

f''(t)=(f(t+2*dt)-f(t))/(dt*dt)
で計算できそうな気がするのですが、こうすると、
(dt*dt)の影響が大きくなりすぎて、計算結果が爆発的に大きくなってしまいます。
恐らく無限小の次数のようなものの影響なのではないかと思うのですが…。

対処方法を教えてください。よろしくお願いします。

■関連質疑
http://okwave.jp/qa/q1903228.html

回答No.1

中村 拓男（@tknakamuri） · ベストアンサー

>f'(t)=(f(t+dt)-f(t))/dt
>f'(t+dt)=(f((t+dt)+dt)-f(t+dt))/dt

>なので、結局

>f''(t)=(f(t+2*dt)-f(t))/(dt*dt)

はおかしいです。

f''(t)= (f'(t+dt) - f'(t)) / (dt) = (f(t+2*dt) + f(t) - 2*f(t+dt))/(dt*dt)

が正しい。

お礼

回答ありがとうございました。

おっしゃる通りです。大変恥ずかしいミスをしていました。
このまま、このスレを閉じてしまうと、単なるスレ汚しになってしまうので、
後で「補足」のところで、サンプルプログラム（VBA）と、簡単な考察を書きます。
のちに二階微分を計算しなきゃいけなくなった人のため…。

補足

中心差分、前進差分、後退差分を含むＥｘｃｅｌＶＢＡのプログラムを作成しました。遅くなりもうしわけありません。

【Ｅｘｃｅｌのワークシート上の設定】
ExcelのCellには、
B4セルにtの値が、
Ｃ４には関数が、例えば=sin(B4)のように入力されているものとします。
D４セルには、dtの値が、例えば0.00001のように入力されているものとする。

【表示結果】
下記マクロを実行すると、

Ｅ4セルに中心差分による一階微分の値が表示されます。
F4セルに中心差分による二階微分の値が表示されます。

Ｇ4セルに前進差分による一階微分の値が表示されます。
Ｈ4セルに前進差分による一階微分の値が表示されます。

Ｉ4セルに前進差分による一階微分の値が表示されます。
J4セルに前進差分による一階微分の値が表示されます。

【ソースコード】

Sub 二階微分()

t = Cells(4, 2)
dt = Cells(4, 4)

ft = Cells(4, 3)

'============一階微分============
Cells(4, 2) = t + dt
Lft = Cells(4, 3)

Cells(4, 2) = t - dt
Bft = Cells(4, 3)

dCft = (Lft - Bft) / (2 * dt)
dLft = (Lft - ft) / dt
dBft = (Lft - ft) / dt

Cells(i + 10, 4) = dft
Cells(i + 10, 7) = dgt

'============二階微分============
Cells(4, 2) = t + 2 * dt
LLft = Cells(4, 3)

Cells(4, 2) = t - 2 * dt
BBft = Cells(4, 3)

ddCft = (Lft - 2 * ft + Bft) / (dt ^ 2)
ddLft = (LLft - 2 * Lft + ft) / (dt ^ 2)
ddBft = (ft - 2 * Bft + BBft) / (dt ^ 2)

'============表示============
Cells(4, 5) = dCft
Cells(4, 6) = ddCft

Cells(4, 7) = dLft
Cells(4, 8) = ddLft

Cells(4, 9) = dBft
Cells(4, 10) = ddBft

End Sub

回答No.2

Tacosan

きちんと紙の上で計算すれば #1 の通り.

実際には中心差分
f''(t) = [f(t+dt) - 2f(t) + f(t-dt)]/dt^2
で計算する方がちょっとだけうれしい.

お礼

回答ありがとうございました。
大変恥ずかしいミスをしていました。後に、「補足入力」にて、
ExcelVBAのサンプルプログラム（中心差分版）を作成し、
後の人のためになるかならないかということで、公開します。

補足

下に、中心差分、後退差分も含んだソースコードを書きました。
遅くなりもうしわけありません。

後の人のために、参考にしたものを、以下に補足的に記載します。

■参考にしたWeb上の記事
http://ja.wikipedia.org/wiki/%E5%B7%AE%E5%88%86%E6%B3%95
http://d.hatena.ne.jp/nowokay/20080516/1210944893
http://d.hatena.ne.jp/nowokay/20080516/1210944893
http://okwave.jp/qa/q5935311.html
http://okwave.jp/qa/q4942769.html
http://workingout.jugem.jp/?cid=4