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

このQ&Aのポイント
  • 二階微分を数値計算する際に問題が生じることがあります。
  • 一階微分と同様の計算方法を用いると計算結果が爆発的に大きくなることがあります。
  • 無限小の次数の影響が計算結果に大きく関与している可能性があります。
回答を見る
  • ベストアンサー

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

関数(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

>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) が正しい。

Kokorochaniuna
質問者

お礼

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

Kokorochaniuna
質問者

補足

中心差分、前進差分、後退差分を含むExcelVBAのプログラムを作成しました。遅くなりもうしわけありません。 【Excelのワークシート上の設定】 ExcelのCellには、 B4セルにtの値が、 C4には関数が、例えば=sin(B4)のように入力されているものとします。 D4セルには、dtの値が、例えば0.00001のように入力されているものとする。 【表示結果】 下記マクロを実行すると、 E4セルに中心差分による一階微分の値が表示されます。 F4セルに中心差分による二階微分の値が表示されます。 G4セルに前進差分による一階微分の値が表示されます。 H4セルに前進差分による一階微分の値が表示されます。 I4セルに前進差分による一階微分の値が表示されます。 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

その他の回答 (1)

  • Tacosan
  • ベストアンサー率23% (3656/15482)
回答No.2

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

Kokorochaniuna
質問者

お礼

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

Kokorochaniuna
質問者

補足

下に、中心差分、後退差分も含んだソースコードを書きました。 遅くなりもうしわけありません。 後の人のために、参考にしたものを、以下に補足的に記載します。 ■参考にした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

関連するQ&A

  • 全微分で2階導関数を求めるについて

    Z=f(X,Y),X=φ(t)Y=Ψ(t)がC2級のときtの関数Z=f(φ(t),Ψ(t))の2階導関数を求める問題の解き方が分かりません。 1階導関数は、全微分の公式を用いてすぐ求めることができるのですが、2階導関数の場合、d/dtを両辺に掛けたとき1階導関数を求めたときに式に現れる∂z/∂xに対してどのように式を変形していけばいいのかわかりません。回答よろしくお願いします。

  • 二階の全微分について

    物理でxyの座標を極座標に変換し加速度を計算するなかで、2階の全微分に困っています。あまり、微分積分は慣れていないので、丁寧に教えていただけると助かります。 http://okwave.jp/qa/q2707943.html でも、同じような質問があります。 一階の全微分はわかりますが、2階の全微分で項が増えるのがわかりません。 具体的には、 Z=f(X,Y), X=g(t) Y=h(t)で、 dZ/dt=(∂Z/∂x)dx/dt+(∂Z/∂y)dy/dt まではよくわかり、これを二階にするときはまず、第1項目(∂Z/∂x)dx/dtが {∂/∂x(∂Z/∂x)dx/dt}dx/dt+{∂/∂y(∂Z/∂x)dx/dt}dy/dt となるだと思うのですが、(∂Z/∂x)d/dt(dx/dt)という項も加わるようです。詳しくその考え方を教えていただけますでしょうか?

  • 媒介変数表示での二階微分に関して

    媒介変数表示の二階微分を計算していたときに疑問がでたので質問します。 x,yがともにtの関数とします。 このときなぜyのxでの二階微分を、次のように計算してはならないのか教えてください。 d^2 y /dt^2 ---------- d^2 x /dt^2

  • 3階微分って何がわかるの??

    微分法の応用にて、関数f(x)の1階微分の正負では関数f(x)の増減が分かり、2階微分では関数f(x)の凹凸がわかるところまでは理解したのですが・・・3階微分すると何が分かるのか分からないので教えてください。

  • 数(3)の微分についてです。

    媒介変数で表された関数の微分法についてなのですが、教科書に下のような説明が書いてあります。 x=f(t),y=g(t)と表され、x,yがtについて微分可能のとき 合成関数の微分法により dy/dx=dy/dt*dt/dx ・・・(1) したがって dy/dx=dy/dt*1/dx/dy=dy/dt/dx/dt=g`(t)/f`(t) (1)の合成関数の微分っていうのはyがtで微分できて、tがxで微分できるときに使えるんですよね?てことはyがtの関数で、tはxの関数で無ければならないと思うのですが、最初に与えられているのはyはtの関数、xはtの関数ってことだけで、tはxの関数であるとは限らないと思うのです。なので上の証明はx=f(t)の逆関数が存在する時しか成り立たないのではないのでしょうか?何故いつも成り立つのかがわかりません。 初歩的な質問ですみませんm(__)m

  • 微小値の微分

    数学の微分に関する質問です。 微小な値を微分すると、その結果もまた微小値となるのでしょうか? 例えば、時間tの関数θ(t)という関数があるとして、 θ(t)を微小値とするとき(|θ(t)|<<1)、θ(t)を微分したd{θ(t)}/dtもまた微小値と言えるのでしょうか? 僕は微小値(この場合θ(t))は微小であり続けるため、あまり大きく変化することはできない、つまり微小な値を微分すると、その結果もまた微小値になると思っています。 この問題はどのように考えればよいのでしょうか? どなたかご教示願います

  • 偏微分

    偏微分を用いて、全微分をするとき 例えばx,y,zの時間に依存する変数からなる関数f(x,y,z)を時間で全微分する時、 df/dt=(df/dx)(dx/dt)+(df/dy)(dy/dt)+(df/dz)(dz/dt) となると思うのですが、 仮に、x,を時間だけでなく、もう一つ時間に依存する関数n(t)を与えるとします、 つまり X=x+n(t) f(x) => f(X)=f(x+n(t)) になるとします。 その時、時間の全微分はどうなるのでしょうか? f(x+n(t))はxとn(t)に依存しているので、f(x,n(t))と書いて f(x+n(t))=f(x,n(t)) df(x+n(t))/dt=(df(x,nt)/dt)=(df/dx)(dx/dt)+(df/dn)(dn/dt) としてもいいんでしょうか? 後どのような時、偏微分しても可能なのか教えて頂ければ幸いです。 どなたか分かる方よろしくお願いします。

  • 2つの関数の積を数値計算で解く

    時間tで微分すると、共にtの関数であるS(t)、I(t)の積の形で表せる関数F(t)を 数値計算で解きたいのですが、やり方がわかりません。 関数が1つだけの時のやり方はわかるのですが、関数が2つ、しかも積となるとお手上げで…。 式の形は dF(t)/dt=-αS(t)I(t)   (αは係数) となります。このような形の式を数値計算で解くやり方を解説しているページ、 または書籍をご存知の方がいらっしゃいましたらご教授願います。

  • 非斉次な2階の線形微分方程式

    fを実数として、非斉次な2階の線形微分方程式を積分定数を用いて解け。(xはtの関数) (d^2x/dt^2)+9x=fe^(2it) x=Ae^(2it)とすると、 A=f/5 ここからどうすればよいかわかりません。 詳しい解説お願いします。

  • 非斉次な2階の線形微分方程式

    fを実数として、非斉次な2階の線形微分方程式を積分定数を用いて解け。(xはtの関数) (d^2x/dt^2)+9x=fe^(3it) x=Ae^(3it)とすると、 0A=f ここからどうすればよいかわかりません。 詳しい解説お願いします。