# Kinetic Friction in Mathematica, weird behaviour

I found the behavoiur of Sign function weird in the code below.

when $T=10$

T = 10;
A = 1;
Fk0 = 2;
m1 = 1;
x0[t_] = A Sin[2 Pi t];
Fk[t_] = -Fk0  Sign[x1'[t]];
sol = NDSolve[{
m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0},
{x0, x1}, {t, 0, T}, MaxSteps -> Infinity]
Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True]


but when $T=40$ , I got

which is how it is supposed to be. I wonder what I did wrong here.

-
no errors/warnings? – george2079 Feb 2 '14 at 13:53
It returns no error. – Dadan Ari Wibowo Feb 2 '14 at 13:57

Using WhenEvent[] to get a typical kinetic friction behavior:

A = 1;
m1 = 1;
x0[t_] := A Sin[2 Pi t];
Fk0 = m1;
Fk[t_] := -Fk0 a[t];
sol[T_?NumericQ] :=
NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0, a[0] == 1,
WhenEvent[x1'[t] == 0, a[t] -> Sign[x1'[t]]]}, {a, x1}, {t, 0, T},
MaxSteps -> Infinity, DiscreteVariables -> a, MaxStepSize -> 10^-3]
GraphicsRow@{
Plot[{x0''[t], x1[t] /. sol[10], x1'[t] /. sol[10]}, {t, 0, 10}, Frame -> True,
Evaluated -> True,  PlotStyle -> {Red, Blue, Green}],
Plot[{x0''[t], x1[t] /. sol[40], x1'[t] /. sol[40]}, {t, 0, 40},
Frame -> True, Evaluated -> True, PlotStyle -> {Red, Blue, Green}]}


-
It could be interesting to add a static friction component at those instants when x1'[t]==0 – Dr. belisarius Feb 2 '14 at 16:55
does WhenEvent detect a sign change, or only trigger if the numerical value is precisely zero? – george2079 Feb 2 '14 at 20:02
@george2079 You can adjust the detection method by using options, but NDSolve being a numerical method implementation can't rely in exact values – Dr. belisarius Feb 2 '14 at 20:04
@george2079 I believe it detects a sign change. – Daniel Lichtblau Feb 2 '14 at 21:43
@george2079: it is documented to detect sign changes: "f==0: the function f of the solution variables crosses zero" (from ref/WhenEvent). – Albert Retey Feb 3 '14 at 9:34

As it turns out the only thing wrong with the original solution is a stroke of bad luck.. Once you hit the zero point at 5, the default time step is one and the solution marches out values at 6,7,8,9,10.. where the forcing function is zero..

Simply limiting the step size fixes things..

 T = 10;
A = 1;
Fk0 = 2;
m1 = 1;
x0[t_] = A Sin[2 Pi t];
Fk[t_] = -Fk0 Sign[x1'[t]];
sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10,
x1[0] == 0}, {x0, x1}, {t, 0, T}, MaxSteps -> Infinity,
MaxStepSize -> .01]
Plot[{x0''[t], x1'[t], 10 Sign[x1'[t]]} /. sol, {t, 0, T},
Evaluated -> True, Frame -> True]


Here is the original showing the actual solution points..

Show[{Plot[x1'[t] /. sol, {t, 0, T}, Evaluated -> True, Frame -> True],
Graphics[{PointSize[.01],
Point@{#, x1'[t] /. First@sol /. t -> # } & /@
First@Cases[sol, f_InterpolatingFunction, Infinity][[1, 3]]}]}]


By the way there is a cleaner way to display the solution points using EvaluationMonitor described in the documentation.

-
In v8.0.4, MaxStepSize -> .01 doesn't work, but MaxStepSize -> .03 does. I think that a more strict control for the step is needed, though I don't know how to achieve that… – xzczd Feb 4 '14 at 10:30
yes 9.0.1. I find it frankly a bit disturbing that the default time step algorithm has made such a poor choice here abruptly jumping by a factor of ~1000. There should be a way to limit the change in time step instead of just setting an absolute limit but I can't readily see a way. – george2079 Feb 4 '14 at 19:20
I found it, see the 3rd approach in my answer. – xzczd Feb 5 '14 at 8:28

Dadan Ari Wibowo, george2079 and belisarius are all in v9, right? With this question, I noticed some slight changes are made for NDSolve in v9.

The syntax of NDSolve in OP's code is incorrect for v8.0.4. x0, as a known function, isn't allowed to be the 2nd argument of NDSolve.

Then the behavior of NDSolve at $t=5$ is different:

T = 6;
A = 1;
Fk0 = 2;
m1 = 1;
x0[t_] = A Sin[2 Pi t];
Fk[t_] = -Fk0 Sign[x1'[t]];
sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0}, x1, {t, 0, T}]
Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True]


Notice that the option MaxSteps -> Infinity has been removed or it'll eat up all my computer memory.

Finally, I've found 3 solutions for this problem. The 1st one is to use option WorkingPrecision -> 16:

sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0},
x1, {t, 0, T}, WorkingPrecision -> 16]
Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True]


I have to mention that only WorkingPrecision -> 16 will work in this case i.e. even WorkingPrecision -> 15 or WorkingPrecision -> 17 will generate the same error. This solution is probably a coincidence caused by numerical error so it's not that instructive.

The 2nd approach is to use Method -> "StiffnessSwitching":

sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0},
x1, {t, 0, T}, Method -> "StiffnessSwitching"]
Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True]
(* The plot is same as the former so I'd like not to "repost" it. *)


I found this approach mainly by… er… intuition so I can't give a theoretical explanation for it. (From the failing plot we see that something seems to happen at zero point so I guess it might be a stiffness issue and it seems to be! ) Maybe someone more experienced in NDSolve can give a detailed instruction?

Inspired by this post, I've found the 3rd solution:

sol = NDSolve[{m1 x1''[t] == -m1 x0''[t] + Fk[t], x1'[0] == 10, x1[0] == 0}, x1, {t, 0, T},
Method -> {"FixedStep", Method -> "ExplicitRungeKutta"},
StartingStepSize -> 1, MaxStepFraction -> Infinity]
(* Plot[{x0''[t], x1'[t] /. sol // Evaluate}, {t, 0, T}, Frame -> True] *)
<< "DifferentialEquationsNDSolveUtilities"
StepDataPlot[sol]


Here I deliberately chose 1 as the step size so all of the zero point should have been hit, but correct result is still generated, so it's hard (at least for me) to tell what actually happened inside NDSolve, after all, it's one of the most complex functions in Mathematica.

Not sure if my 3 approaches will work for v9.

-