1

下のグラフをトレースする Mathematica コードを知っている人はいますか?

グラフの方程式は次のとおりです。これは、係数が一定の 2 次線形微分方程式です。

ここに画像の説明を入力

この方程式によってトレースされたグラフは次のとおりです。

ここに画像の説明を入力

本「Times Series Analysis and Forecasting By Example」からの引用:

... ここで δ(t ) はインパルス (デルタ) 関数で、エンドウ豆のショットのように、時間 t = 0 で振り子をその平衡から遠ざけ、a はエンドウ豆による衝撃の大きさです。この 2 次微分方程式によってトレースされる曲線が時間の減衰正弦波関数であることは容易に想像できますが、摩擦または粘性が十分に大きい場合、(過減衰) 振り子は、交差することなく指数曲線に従って徐々に停止する可能性があります。中心線。

4

1 に答える 1

6 
eq = m z''[t] + c z'[t] + k z[t] == a DiracDelta[t];
parms = {m -> 1, c -> .1, k -> 1, a -> 1};
sol = First@DSolve[{eq /. parms, z[0] == 1, z'[0] == 0}, z[t], t];
Plot[z[t] /. sol, {t, 0, 70}, PlotRange -> All, Frame -> True, 
 FrameLabel -> {{z[t], None}, {Row[{t, " (sec)"}], eq}}, 
 GridLines -> Automatic]

Mathematicaグラフィックス

初期条件がゼロの場合、別のオプションは、次のようにMathematicaの制御システム関数を使用することです。

parms = {m -> 10, c -> 1.2, k -> 4.3, a -> 1};
tf = TransferFunctionModel[a/(m s^2 + c s + k) /. parms, s]
sol = OutputResponse[tf, DiracDelta[t], t];

Plot[sol, {t, 0, 60}, PlotRange -> All, Frame -> True, 
 FrameLabel -> {{z[t], None}, {Row[{t, " (sec)"}], eq}}, 
 GridLines -> Automatic]

Mathematicaグラフィックス

アップデート

厳密に言えば、DSolve上記の結果は、この問題の手作業による導出で見つけることができるものではありません。正しい解決策は次のようになります

（これも参照してください）

正しい解析解は次の式で与えられます。

Mathematicaグラフィックス

これは、この問題と同様のケースについてここ（最初の章）で導き出しました。

上記のソリューションを使用すると、正しい応答は次のようになります。

parms = {m -> 1, c -> .1, k -> 1, a -> 1};
w = Sqrt[k/m];
z = c/(2 m w);
wd = w Sqrt[1 - z^2];
analytical = 
  Exp[-z w t] (u0 Cos[wd t] + (v0 + (u0 z w))/wd Sin[wd t] + 
     a/(m wd) Sin[wd t]);
analytical /. parms /. {u0 -> 1, v0 -> 0}

 (* E^(-0.05 t) (Cos[0.998749 t] + 1.05131 Sin[0.998749 t]) *)

それをプロットする：

Plot[analytical /. parms /. {u0 -> 1, v0 -> 0}, {t, 0, 70}, 
 PlotRange -> All, Frame -> True, 
 FrameLabel -> {{y[t], None}, {Row[{t, " (sec)"}], 
    "analytical solution"}}, GridLines -> Automatic, ImageSize -> 300]

Mathematicaグラフィックス

を使用して上記のプロットを上記の最初のプロットと比較すると、のDSolve近くの違いを確認できますt=0。

于 2013-01-06T02:13:52.560 に答える