Mathematica代码绘制这个微分方程的图形？

Question

有谁知道将跟踪下图的Mathematica代码？

这是图的等式,具有常系数的二阶线性微分方程:

以下是此等式跟踪的图表:

从"时间序列分析和预测示例"一书中引用:

......其中δ(t)是脉冲(δ)函数,就像豌豆击中一样,在时间t = 0时迫使摆远离其平衡,a是豌豆撞击的大小.很容易想象,这个二阶微分方程所追踪的曲线是时间的阻尼正弦函数,但是,如果摩擦或粘度足够大,(过阻尼)摆可能会在指数曲线之后逐渐停止而没有交叉中心线.

Answer 1

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中的控制系统函数,如下所示

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]

更新

严格地说,上述结果DSolve不是手工推导出这个问题的结果.正确的解决方案应如下所示

(见这也供参考)

正确的分析解决方案由

我在这里得到了这个问题和类似案例(第一章).

使用上面的解决方案,正确的响应将如下所示:

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]

如果您将上面的图表与上面显示的第一个图表进行比较,DSolve您可以看到附近的差异t=0.