# Harmonic oscillator with damping attrition

Considering the following physical situation:

and writing the following code:

m := 1.52
g := 9.81
us := 0.15
uk := 0.10
k := 2.12
xi := 4.00
vi := 0.00
tmax := 10

P := m g
Fs := us P
Fk := uk P
Fe[t_] := -k x[t]

sol = NDSolve[{
Fe[t] - Sign[x'[t]] Fk == m x''[t],
x[0] == xi,
x'[0] == vi},
x, {t, 0, tmax}];

Plot[Evaluate[{
Sign[x'[t]] Fs,
Fe[t],
x[t]} /. sol],
{t, 0, tmax},
AxesLabel -> {"t", "fct[t]"},
PlotLegends -> {"-Fs", "Fe", "x"}]


you get the following graph:

which shows that oscillations are over for $t \approx 8\,s$ causes kinematic friction.

On the other hand, putting us = 0.40 you get this other graph:

which shows that oscillations are over for $t \approx 3\,s$ causes static friction.

Question: is it possible to automate all this by making the x(t) graph plot until the motion stops?

• Ok, there are two situations with different amplitudes of you resulting functions, and what? Could you clarify, what namely do you want to automate? Aug 14 '17 at 11:42
• @Rom38: I would like to write a code that, in reference to the two examples mentioned above, plots for 0 <= t <= 8 and 0 <= t <= 3, i.e. only when the motion is actually established.
– TeM
Aug 14 '17 at 12:05
• Two questions: You use Fk in NDSolve but plot Fs. Is this intended? Secondly, what exactly happens at 8s resp. 3s? I don't see anything significant at those times in the plots you provided. (Also, x[t] seems nearly identical between the two plots, is this intended?) Aug 14 '17 at 12:13
• @Mathe172: The answers are in the physical phenomenon. When the body is in motion it is subject to the elastic force Fe and the kinematic friction Fk, at points where the velocity is zero (motion reversal) the body is subject to the elastic force Fe and the static friction Fs. For this reason, the body can stop for two reasons: either in a generic point causes the kinematic friction Fk or in the inversion points to cause static friction Fs.
– TeM
Aug 14 '17 at 12:24
• Maybe you added the wrong image, because even in your second graph the motion continues till 8 s.
– user484
Aug 14 '17 at 12:34

This is a perfect use case for WhenEvent:

sol = NDSolve[
{
Fe[t] - Sign[x'[t]] Fk == m x''[t],
x[0] == xi,
x'[0] == vi,
WhenEvent[x'[t] == 0 && Fs > Abs[k x[t]], tmax = t; "StopIntegration"]
}
, x, {t, 0, Infinity}]


Note that this automatically sets tmax, so you don't need to specify anything before. The only thing to note is that you can't replace k x[t] with Fe[t], since WithEvent doesn't see the x[t] in that case. You could write (note the Evaluate wrapped around the condition)

WhenEvent[Evaluate[x'[t] == 0 && Fs > Abs[Fe[t]]], tmax = t; "StopIntegration"]


if you really want to write Fe[t].

For us=0.15:

For us=0.4:

• Fantastic! I had thought of something more artisan to establish tmax before the plot, but this solution is definitely more elegant, exactly what I was looking for! I also thank all the others, unfortunately I do not know English well, I find it difficult to explain the basic physical phenomenon!
– TeM
Aug 14 '17 at 13:03
• Glad it works for you. Please don't forget to mark your question as answered. Aug 14 '17 at 13:09
• Sorry if I break still but I still have a request. If us = 0 I would like WhenEvent to "disappear" and tmax = 10. How could I do it?
– TeM
Aug 14 '17 at 20:47
• I can't test it right now, but something like WhenEvent[If[us != 0, x'[t] == 0 && Fs > Abs[k x[t]], t >= 10], tmax = t; "StopIntegration"] should do the trick Aug 14 '17 at 20:52
• Incredible, you know how to answer any question, you are too strong and generous too. Thank you!
– TeM
Aug 14 '17 at 21:00