# Ndsolve problem with Sign (friction)

In a very simple model I consider the motion of a body with frictional force.

X = ParametricNDSolveValue[{x''[t] ==1 -    x[t] - \[Mu] Sign[x'[t]]  Abs[1 + x[t]], x[0] == 0,x'[0] == 0}, x, {t, 0, 50}, \[Mu]]


Setting x'[t]==0, x''[t]]==0 I would expect the stationary solution x[t]==1 but Mathematica evaluates

Plot[Table[X[\[Mu]][t], {\[Mu], { .1, .2, .5 }}], {t, 0, 50}, GridLines -> {None, {1}},PlotRange->All]


different asymptotic solutions depending on \[Mu]

What's going wrong here?

I know that NDSolve sometimes shows problems with Sign[..]-function and substituted the Sign-function

X = ParametricNDSolveValue[{x''[t] ==1 -    x[t] - \[Mu]  Tanh[x'[t]]  Abs[1 + x[t]], x[0] == 0,x'[0] == 0}, x, {t, 0, 50}, \[Mu]]
Plot[Table[X[\[Mu]][t], {\[Mu], {.1, .2, .5}}], {t, 0, 50},GridLines -> {None, {1}}, PlotRange -> All]


Now the result shows, as expected, a unique stationary solution for different \[Mu]

My questions:

• Isn't NDSolveable to solve non-smooth odes?
• Can I use special methods to avoid this problem?
• Shouldn't I use Sign and NDSolve together?

Thanks!

• You may want to look at WhenEvent>Applications>Friction Models. It uses Sign. Maybe it can be adapted to your purpose. – Tim Laska Jun 19 at 13:33
• Thanks, I'll take a look. – Ulrich Neumann Jun 19 at 13:35
• You have an oscillatory normal force. At some point the friction force will be greater than the restoring force and then the mass will be stuck. There is no reason why this should be at the same asymptotic displacement. If the normal force was constant the point of sticking will vary with the initial conditions. – Hugh Jun 19 at 16:49
• @Hugh Thanks. I didn't model stip slick effect. My argument expecting a unique asymptotic solution is simple substituting x'==x''==0into the ode. – Ulrich Neumann Jun 19 at 18:40
• What signs would you use for 0 and what about Abs? I don't think you can just set velocity and acceleration to zero. This answer is the correct approach. – Hugh Jun 19 at 18:44

## 2 Answers

As @Hugh points out, the particle will remain at rest if it comes to rest at a point where the restoring force is less than or equal to the frictional force. The interval where this occurs may be solved for (cps below).

ics = {x[0] == 0, x'[0] == 0};
ode = {x''[t] == 1 - x[t] - μ Sign[x'[t]] Abs[1 + x[t]]};

acc = x''[t] /. First@Solve[ode, x''[t]] /. {x'[t] -> v, x[t] -> x};
{accm, accp} = Simplify[acc /. Sign[_] -> {-1, 1}, x > 0];
cps = x /. Solve[# == 0, x] & /@ {accm, accp} // Flatten
(*  {(-1 - μ)/(-1 + μ), (1 - μ)/(1 + μ)}  *)


The physical intuition about the forces shows up in the NDSolve code in the phase portrait of the ODE. Along the interval defined by cps (red line), the vector field points in opposite directions and is normal to the red line. The Filippov sliding mode continuation results in no movement (that is, x[t] remains constant once the system is in a state along the red line).

Block[{μ = 0.1},
With[{a = x''[t] /. First@Solve[ode, x''[t]] /. {x'[t] -> v, x[t] -> x}},
StreamPlot[{v, a}, {x, 0, 2}, {v, -1, 1},
FrameLabel -> {HoldForm[x], HoldForm[v]},
Epilog -> {Red, Line[Transpose@{cps, {0, 0}}]}]
]]


In short, the solution with Sign is correct.

Reply to comment, which is too long for a comment:

@UlrichNeumann writes: One point I want to mention. If a stationary solution on the "red line" is calculated , obviously x'[t]==0. Still this solution must fullfill the ode, in the stationary case 0== 1-x[t] which implies xRED[t]==1 ???

To handle the sort of paradox or contradiction in discontinuous differential equations that you note, one has to extend one's notion of what such an DE and its solutions are. (Aside: The Green's function is perhaps the earliest example, but not particularly relevant here.) Here, Mma is applying Filippov's idea of replacing the differential equation by a differential inclusion. This is a little involved and I'm not really an expert. I'll skip some of the technical details. I hope it's clear and accurate enough.

Consider this formulation of your ODE: \eqalign{ \dot x &= v\cr \dot v &= f(x,v) = 1 - x - \mu \mathop{\text{sgn}}(\dot x) \,|{1 + x}| \cr }\tag{1} Let $$X = (x,v) \equiv (x, \dot x)$$. Then in the Filippov theory, the phase vector $$\dot X = (\dot x, \dot v)$$ should belong to a set $$F(X)$$, and $$\dot X \in F(X)$$ is called a differential inclusion. The set $$F(X)$$ is defined according to values of the vector field $$\dot X$$ in a neighborhood of $$X=(x, v)$$. In this case, Filippov's construction (details omitted) leads to the following definition of $$F$$: $${\rm(a)}\ F(x,0)=\{(0,a) \colon a_1 \le a \le a_2\}\quad \text{and}\quad {\rm(b)}\ F(x,v)=\{(v,f(x,v))\}\ \text{for}\ v\ne0\,,\tag{2}$$ where the interval $$[a_1,a_2]$$ is discussed in the next paragraph. Note $$\dot X \equiv (\dot x, \dot v) \in F(x,v)$$ means that when $$v\ne0$$, the differential equation (1) is satisfied, which always happens at points where the vector field is continuous.

The case that interests us, of course, is the case when $$v=0$$, where the vector field is discontinuous. In that case, when $$(x,v)$$ is on the red line, we have $$a_1 \le 0 \le a_2$$, so that by (2a) above $$a_1 \le \dot v \equiv \ddot x \le a_2$$ and $$\dot v = 0$$ is permitted by the differential inclusion $$\dot X \in F(x,0)$$. (In fact $$a_1 = \mu +(\mu -1)\, x+1$$ and $$a_2 = -\mu -(\mu +1)\, x+1$$.) Outside the red line, still along $$v=0$$, $$a_1$$ and $$a_2$$ have the same sign. They are both positive to the left of the line, which is reflected by the upward flow in the phase portrait; and the situation is the opposite to the right.

I'll leave aside why, out of all the choices in $$[a_1,a_2]$$ on the red line when $$v=0$$, the acceleration $$\ddot x$$ should be zero in this case according to the Filippov theory. It is obvious from the intuitive physical setup that it should be. I hope that an outline of a theory of discontinuous DEs that allows $$\ddot x$$ to be zero, and thus gets around the apparent contradiction, is sufficiently helpful.

For more, see the tutorial linked above and search for "sliding mode"; also Filippov, A., Differential Equations with Discontinuous Right Hand Sides, Kluwer Academic Publishers, 1988. Somewhere I maybe should have remarked that the Filippov solution implemented by NDSolve is part symbolic processing and part (mostly) numerical integration.

• @ MichaelE2 Thank you for your detailled explanation. Especially the phasediagram is a very helpful modell which I frequently use. I will analyse your approach in detail. One point I want to mention. If a stationary solution on the "red line" is calculated , obviously x'[t]==0. Still this solution must fullfill the ode, in the stationary case 0== 1-x[t] which implies xRED[t]==1 ??? – Ulrich Neumann Jun 19 at 18:51
• @UlrichNeumann Its been a long time since I've done ODEs, so this is just a guess! but imagine as x'[t] approaches arbitrarily close to zero, but Sign[x'[t]] is still positive. The ODE in this "stationary" case is then 0=0 and doesn't constrain the solution in the same way. Someone who knows better than me will need to weigh in, but I can imagine that a numerical solution could be affected by this. – mbrig Jun 19 at 23:10
• @UlrichNeumann I replied, briefly?, in an update to my answer. – Michael E2 Jun 20 at 1:32
• @ MichaelE2 Thank you very very much. I have to elaborate your hints... What I'm missing, perhaps not understanding, is the fact Sign[0]==0 in the friction modell. – Ulrich Neumann Jun 20 at 9:32
• @UlrichNeumann In the Filippov method, the value of Sign[0] does not really matter since the function is discontinuous there. If you want to stick to the classical formulation of ODEs as equations to be satisfied at all times, then you run into all sorts of trouble. For instance, there can be no solution x[t] such that x''[t] satisfies the ODE at the time when x'[t] == 0, because of the type of discontinuity of Sign[x] at x == 0. (If there's a finite jump in the value of a derivative when crossing t == a, then the derivative must undefined at t == a by basic analysis/calculus.) – Michael E2 Jun 20 at 17:58

Not an answer, but Tanh also approaches the Sign solution as one sharpens the transition so maybe Sign is behaving appropriately.

X = ParametricNDSolveValue[{x''[t] ==
1 - x[t] - μ Tanh[ m x'[t]] Abs[1 + x[t]], x[0] == 0,
x'[0] == 0}, x, {t, 0, 50}, {μ, m}]
plts = Table[
Plot[Evaluate@Table[X[μ, 10^lm][t], {μ, {.1, .2, .5}}], {t,
0, 50}, GridLines -> {None, {1}}, PlotRange -> {0, 2},
PlotLegends -> Automatic], {lm, -3, 4, 0.1}];
ListAnimate[plts]


• Good point, thank you. I tried too to find the optimal factor inside Tanh[]. For me it seems to be a numerical problem, don't know how to tell it NDSolve – Ulrich Neumann Jun 19 at 18:28