CODEWITHSUNDEEP
wolfram-mathematica
Contango
Contango

Reputation: 80222

Mathematica code to draw a graph of this differential equation?

Does anyone know the Mathematica code that will trace the graph below?

Here is the equation for the graph, a second order linear differential equation with constant coefficients:

enter image description here

Here is the graph traced by this equation:

enter image description here

Quote from the book "Times Series Analysis and Forecasting By Example":

... where δ(t ) is an impulse (delta) function that, like a pea shot, at time t = 0 forces the pendulum away from its equilibrium and a is the size of the impact by the pea. It is easy to imagine that the curve traced by this second order differential equation is a damped sinusoidal function of time although, if the friction or viscosity is sufficiently large, the (overdamped) pendulum may gradually come to rest following an exponential curve without ever crossing the centerline.

Upvotes: 1

Views: 1961

Answers (1)

Nasser
Nasser

Reputation: 13131

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 graphics

Notice that, for zero initial conditions, another option is to use the Control system functions in Mathematica as follows

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 graphics

Update

Strictly speaking, the result of DSolve above is not what can be found by hand derivation of this problem. The correct solution should come out as follows

(see this also for reference)

The correct analytical solution is given by

Mathematica graphics

which I derived for this problem and similar cases in here (first chapter).

Using the above solution, the correct response will look like this:

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

Plotting it:

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 graphics

If you compare the above plot with the first one shown above using DSolve you can see the difference near t=0.

Upvotes: 6

Related Questions