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

解決済

質問者：Kokorochaniuna
質問日時：2011/07/28 13:37
回答数：2件

関数(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://oshiete.goo.ne.jp/qa/1903228.html

通報する

この質問への回答は締め切られました。

質問の本文を隠す

回答 (2件)

ベストアンサー優先
最新から表示
回答順に表示

No.1ベストアンサー

回答者： tknakamuri
回答日時：2011/07/28 14:20

>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)

が正しい。

この回答への補足

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

【Ｅｘｃｅｌのワークシート上の設定】
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