Closed-form solutions to $x''+\frac{k}{m}\ x+\mu\ g\ \text{sgn}(x')=0$


I am looking for simple mechanics models that could have closed-form solutions that achieves finite extinction times where it becomes zero for their own system dynamics and stays there forever after.

Looking for simple systems I found in Wikipedia that a mass sliding in a horizontal plane under Coulomb friction is modeled by the differential equation of the Newton's 2nd law as: $$ m\ x'' = -F - \mu\ m\ g\ \text{sgn}(x')$$ where the mass $m$, the earth's gravity acceleration $g$, and the kinetic coefficient of friction $\mu$ are positive constants.


I found here: brick sliding in an horizontal plane after an initial push that the mass sliding after an initial push indeed slides until it stops moving, being their behavior described through a piecewise polynomial closed-form solution, representing the scenario where the force $F = 0$. As example is pretty obvious, which is good since can be easily found the procedure is right since results matches the classic answers found through energy analysis of the system.

Now, to move into the next step of difficulty, I want to model the exact example of the Wikipedia page, the case where the mass is attached to a spring, where the only additional force present besides the Coulomb damping. So modeling the spring as the classical force $F = k x$ with $k$ the string constant ($k>0$), the previous equation becomes: $$ m\ x'' = -k\ x - \mu\ m\ g\ \text{sgn}(x')$$

The question

I have tried unsuccessful to solve the equation as I did for the case $F=0$, and so far I don't find any papers with closed-form solutions to the equation: $$ x'' = -\frac{k}{m}\ x - \mu\ g\ \text{sgn}(x')$$

If I made every constant equal to one, then Wolfram-Alpha shows the following:


As expected for a non-linear equation there are multiple solutions. I am specially interested in the solutions were the mass stop moving (which is impossible to represent accurately through a linear ODE or non-piecewise power series, as explained here - otherwise it will violate the Identity Theorem), but differently from the mentioned example where I used a self-named endiness constraint: there exists a time $T>0$ such as $x(t) = 0,\ \forall t>T$, but for what I have found on the papers, in this mass-spring system also exists the possibility that the mass stops moving in a different position than the original rest position, so it could have a final position $x_f$ constant, such as $x(t) = x_f,\ \forall t>T$.

I hope you could find the closed-form solution that stops moving, showing the equations that determines the finite extinction time $T$ and the final position $x_f$, showing how you found them.

Added later____________

Thinking in how the system could stop moving in a position different from equilibrium, I am expecting to find a solution of the form: $$f(t) = x(t)\theta(T-t)+x(T)\theta(t-T)$$ for some function $x(t)$ such as $f'(t)$ is a function finite duration (it achieve by self-dynamics the value zero "continuously" and stay there forever after), and with $\theta(t)$ is the Heaviside step function.

Here I got messed with the integration constant by going backwards from $f'(t)$ to $f(t)$ since I don't know how you make appear the term $x(T)\theta(t-T)$ (this is the problem with these special functions that are hidden distributions - if you could explain this also, it would be great). This was finally answered here, but I still don't figure out how the solution will become a constant: for example think in the following numerical aproximation $x''+x+\text{sgn}(x')=0$ with arbitrary initial conditions $x(0)=-\frac{3\pi}{2}$ and $x'(0)=\frac{\pi}{2}$ shown numerically here:

numerical simple example

I don't have any clue of how the solution will satisfy the differential equation after it stops moving, neither how it will match the "discrete-fix-value but variable term" made by $\text{sgn}(x')$. If you could explain it theoretically to have insight I will appreciate it, since I have intuition that the solution will be within the theory of distributions, at least in part (thinking in thing like $xf(x)=xg(x) \Rightarrow f(x)=g(x)+c\delta(x)$).

Also, the solution should have some decaying term as is shown numerically for $x''+x+\text{sgn}(x')=0,\ \ x(0)=3\pi/2,\ \ x'(0)=10$ in here: example of decay

This last point is important: as example, if I use the ansatz $x(t)=c_1\sin(t)+c_2\cos(t)-\text{sgn}\left(c_1\cos(t)-c_2\sin(t)\right)$ for the equation $x''+x+\text{sgn}(x')=0$ and I ignore every distribution-alike term like $\delta(f(t))\equiv 0$ I can see in Wolfram-Alpha that for every value of $t$ the differential equation is fulfilled, but as you could see from the solution the decaying term is missing, so it is not the solution for the problem - here maybe the distribution theory have something to say.

What I did: $$\begin{array}{r l} x''+x+\text{sgn}(x') = 0 & \Biggr| \frac{\partial}{\partial t} \\ \Rightarrow x'''+x'+2 x''\delta(x') =0 & \Biggr| \cdot x'\\ \Rightarrow x'x'''+(x')^2+\require{cancel}\cancel{2x''\cdot\underbrace{x'\delta(x')}_{\text{since }x\delta(x)\ =\ 0}} = 0 & \Biggr| y=x'\\ \Rightarrow y(y''+y)=0 & \Biggr| y\neq 0\ \forall t\\ \Rightarrow y(t) = c_3\sin(t)+c_4\cos(t) & \Biggr| \int \, dt \\ \Rightarrow x(t) = c_1\sin(t)+c_2\cos(t) + c & \end{array}$$ but for inspection one could notice that the integration constant $c$ cannot really be constant since when using the solution in the equation $x''+x+\text{sgn}(x')=0$ I need something to "kill" the term $\text{sgn}(\cdot)$, so the "constant" must not be constant, which it is also weird. And also the decaying term is still missing.

My guess is that the decaying term should be a piece-wise polynomial of the form $(T-t)^q\theta(T-t)$ such as every derivative which rise delta function got canceled by $x^n\delta(x) = 0$, similar to what happened in brick sliding in an horizontal plane after an initial push, but here since the velocity profile becomes zero, but the solution becomes a constant that can be different from zero, it is not easy to figure out how the answer will become a constant value - I hope that maybe this could help someone else to make a/the closed-form solution.

I have checked in online Octave if the decaying behavior were just a mistake of Wolfram-Alpha (since the answer founded looks it solves at least the procedure I use to find it), but it looks is right, both could happen: (i) a decaying behavior, and (ii) and ending position different from equilibrium - both behaviors unattainable by a simple harmonic oscillator.

f = @(t,y) [y(2);-sign(y(2))-y(1)]; 
t0 = 0; y0 = [10;-4];
[ts,ys] = ode45(f,[t0,20],y0,opt);

Plot of the numerical solution to $x''+x+\text{sgn}(x')=0,\,x(0)=10,\,x'(0)=-4$: Numerical simulation in Octave

2nd Added later____________

Thinking in the following form: $$\begin{array}{r c l} \text{let }y & = & x-\text{sgn}(x') \\ \Rightarrow y' & = & x'-\delta(x')x'' \\ \Rightarrow y'' & = & x'' -\delta'(x')(x'')^2-\delta(x')x''' \end{array}$$

Now, since $$\delta(f(x)) = \sum\limits_{n}|f'(x_n)|^{-1}\delta(x-x_n)\text{ with }f(x_n)=0,\,f'(x_n)\neq 0 $$

I will work from now on assuming that every term $\delta(\cdot)\equiv 0$ since they will affect only a zero-measure points, this hoping to find the solution outside this problematic points and later figure out how to solve in this places:

So with this, the previous analysis becomes: $$\begin{array}{r c l} \text{let }y & = & x-\text{sgn}(x') \\ \Rightarrow y' & = & x' \\ \Rightarrow y'' & = & x'' \end{array}$$ So now I have $y''+y+\text{sgn}(y') = 0 \iff x''+ x -\text{sgn}(x')+\text{sgn}(x')=0 \iff x''+x = 0$

But now I don't know how to solve $x''+x=0$ without a only-trigonometric solution. I tried unsuccessfully something of the form: $$ x = (T-t)\theta(T-t)(c_1\sin(t)+c_2\cos(t)) $$ which at least cancel every polynomial term letting only trigonometric functions to be cancelled.

If I plot the following solution in online Octave:

f = @(t,y) [y(2);-sign(y(2))-y(1)]; 
t0 = 0; y0 = [30;-4];
[ts,ys] = ode45(f,[t0,50],y0,opt);

Is not hard to see that the decay is indeed linear, so I don't believe I am too lost about the solution. Hope you could share some ideas. Linear decay of the solution

Finally I have tested the linear decay and it isn't purely lineal (more clear at the end), as can be seen contrasted in the following graph:

Not purely lineal at the end.

3rd attempt________________

Following the change of variables that @eyeballfrog have used in their answer, lets think in the differential equation as: $$z(\tau)''+\frac{k}{m}\ z(\tau)+\mu\ g\ \text{sgn}(z'(\tau))=0$$

Now let $\tau = w t$ with $w = \sqrt{\frac{k}{m}}$, so $\frac{k}{m} = w^2$ and $z(\tau) = \frac{\mu g}{w^2}x(w\ \tau)$. Now I would like to know which happens in $ z(\tau)''+ w^2 z(\tau)+\mu\ g\ \text{sgn}(z'(\tau))=0$.

Then, I have that: $$\begin{array}{r c l} z'(\tau) = \frac{\mu g}{w^2} \frac{\partial}{\partial \tau} x(w\tau) = \frac{\mu g}{w^2} x'(w\tau) w & = & \frac{\mu g}{w}x'(w\tau) \\ \Rightarrow z''(\tau) = \frac{\partial}{\partial \tau}\left( \frac{\mu g}{w}x'(w\tau)\right) = \frac{\mu g}{w} x''(w \tau) w & = & \mu g\ x''(w\tau)\\ \end{array}$$

Replacing I will have now that: $$\begin{array}{c} \mu g\ x''(w\tau)+ w^2 \frac{\mu g}{w^2}x(w\ \tau)+\mu\ g\ \text{sgn}(\frac{\mu g}{w}x'(w\tau))=0 \\ \iff \mu g \left(x''(w\tau) +x(w\tau)+\text{sgn}(\frac{\mu g}{w}x'(w\tau))\right) = 0 \\ \end{array}$$

Now using that $\mu g \neq 0$ and also that $\frac{\mu g}{w}>0$ such its true that

$$\text{sgn}(\frac{\mu g}{w}x'(w\tau)) = \frac{\frac{\mu g}{w}x'(w\tau)}{|\frac{\mu g}{w}x'(w\tau)|} = \frac{\frac{\mu g}{w}x'(w\tau)}{|\frac{\mu g}{w} | \cdot|x'(w\tau)|} = \text{sgn}(x'(w\tau))$$

then I will have I only need to find the solutions to: $$x''(t) +x(t)+\text{sgn}(x'(t)) = 0$$

Which is why I am working in this equation now: the solution will work only if $\mu>0$ and $g>0$ and $k>0$ and $m>0$, keep this in mind.

Now, I have found to make appear a decaying behavior in the solution but I don't think is a $100\%$ rigorous:

$$\begin{array}{r l} x''+x+\text{sgn}(x') = 0 & \Biggr| \frac{\partial}{\partial t} \\ \Rightarrow x'''+x'+2 x''\delta(x') =0 & \Biggr| \cdot \delta(x')\\ \Rightarrow \delta (x')\left(x'''+ x' \right) + 2x''(\delta(x'))^2= 0 & \Biggr| \text{assuming arbitrarily that }(\delta(x))^2\equiv\delta(x)\\ \Rightarrow ?\quad \delta (x')\left(x'''+ x' + 2x''\right) = 0 & \Biggr| y = x' \\ \Rightarrow y''+y+2y'=0 & \\ \Rightarrow y(t) = c_3\exp(-t)+c_4\ t \exp(-t) & \Biggr| \int \, dt \\ \Rightarrow x(t) = c_1\exp(-t)+c_2\ t \exp(-t) + c & \end{array}$$

Which is interesting since indeed got a polynomial term, also jointly with an exponential, but in a way that the derivative keep the same structure except for a constant: this thinking in that the speed profile needs to end at zero at time $T$, but the position needs to become a constant, kind of fit with this result.

Since the decay looks quite linear, I tried for the exponential term something of form $\frac{e^{T}}{e^{T}-1}\left(1-e^{T-t}\right)$ in order it affects little at the beginning, but that ends at zero so it keeps the same decay as the polynomial term, again of the form $(T-t)$ so it have a zero at time $T$, and since the trigonometric function is already tested as solving the solution, I incorporated for now arbitrarily as a linear addition that will be predominant at the end of the movement (should be the same function as the one attached to the decay, at least in principle, in order to keep the number of integration constants that could be determined by the initial conditions):

exponential attempt

Which looks quite promising for being just a not-fitted mix of the solutions I have being found (not rigorously), so I think maybe the answer have the form: $$x(t) = M\theta(T-t)(t-T)(1-e^{-(T-t)})\left(a\sin(t)+b\cos(t)\right)+N\left(a\sin(t)+b\cos(t)\right)+C$$

or like $$x(t) = M\theta(T-t)|(t-T)(1-e^{-(T-t)})|\left(a\sin(t)+b\cos(t)\right)+N\left(a\sin(t)+b\cos(t)\right)+C$$

with some constants $M,\ N,\ T,\ a,\ b,\ C$ to be determined by initial conditions, but I have no clue if the $\text{sgn}(x')$ should be added, while looking good in theory, it will introduce unobserved "jumps" on the solution.

Best numerical attempt so far_________

I uploaded because is quite good, even when I just search by hand some values for the constants, without making any "goodness-of-fit" approach for estimating them:

f = @(t,y) [y(2);-sign(y(2))-y(1)]; 
t0 = 0; y0 = [30;-4];
[ts,ys] = ode45(f,[t0,50],y0,opt); 
x = (30+3/2)*((50-ts)/50)*(exp(50)/(exp(50)-1)).*(1-exp(ts-50)).*cos(ts)-3/2*cos(ts)+1/2*sin(ts);

Best attempt so far

Added later: The solution is not a pure/piecewise trigonometric function

As some of the answers I have got solve the equation by cases, solutions displayed, looking "reasonable" in the math shown, they are made by pure trigonometric functions, or by piecewise constant amplitude decays by half-cycle, which from the numerical solutions could be seen they aren't right, since the plot shows that every lobe is asymmetric proving that the solution is under a gradual decay, as make sense from the physical point of view of the problem, since friction is always present and not piecewise during the movement of the object through time.

As a fast example in online Octave, here the following plot counting approximately the grid squares between zeros and the point of min/max:

f = @(t,y) [y(2);-sign(y(2))-y(1)]; t0 = 0; y0 = [10;-4]; opt=odeset('RelTol',1e-3,'AbsTol',1e-4); [ts,ys] = ode45(f,[t0,20],y0,opt); plot(ts,ys(:,1),'b'), set(gca,'xtick',[0:0.5:20]), set(gca,'ytick',[-10:0.5:10]), grid on;

not pure trigonometric neither piecewise trigonometric

  • $\begingroup$ How are you modelling the sign when $x'=0$? $\endgroup$
    – FShrike
    Commented Feb 8, 2023 at 19:36
  • $\begingroup$ @FShrike as $0$ same as its classic definition Sign function $\endgroup$
    – Joako
    Commented Feb 8, 2023 at 20:25
  • 1
    $\begingroup$ Ok. Well, since $x’$ must be continuous (for it to be second differentiable) so $x’^{-1}(0)$ is closed and the complement is open. We can then safely conceive of differentiation on $\Bbb R\setminus\{x’^{-1}(0)\}$ and on each (open) component you can solve this with the standard simple harmonic motion formula. All that remains to do is piece together the disparate oscillators together on $x’^{-1}(0)$. There can’t really be a single equation closed form since (I think) an arbitrary number of pieces is possible $\endgroup$
    – FShrike
    Commented Feb 8, 2023 at 20:56
  • $\begingroup$ @FShrike From what I have read here, I believe that within the interval given by the time of the "initial condition" $t_0$ and the finite extinction time $T>t_0$ where it stops moving, the solution still hold uniqueness. $\endgroup$
    – Joako
    Commented Mar 20, 2023 at 19:19
  • $\begingroup$ I suggest to investigate a perturbative solution for weak Coulomb friction. As the harmonic oscillator case has been solved, this shouldn't be too difficult. This should lead to a damped oscillation. This is likely a singular perturbation problem since the large time limit and zero coulomb friction limit will not commute $\endgroup$ Commented Apr 6, 2023 at 17:02

5 Answers 5


Let's define the natural frequency $\omega = \sqrt{k/m}$ and the frictional length scale $L = \mu g/\omega^2$, then normalize to $\bar{x} = x/L$ and $\bar{t} = \omega t$. The differential equation takes the considerably simpler form $$ \bar{x}'' + \bar{x} = -\operatorname{sgn}(\bar{x}'). $$ First, suppose the initial condition is $\bar{x}(0) = \bar{x}_1 >0$, $\bar{x}'(0) = 0$. Now if $\bar{x}_1 \le 1$, the spring can't overcome the force of friction and the solution will be $x(t) = x_1$. Otherwise, the object will move towards zero, so we have $\operatorname{sgn}(x') = -1$. That gives the inhomogeneous differential equation $\bar{x}'' + \bar{x} = 1$, which given the intial conditions solves to $$ \bar{x}(t)= 1+\left(\bar{x}_1 - 1\right)\cos(\bar{t}), $$ and in particular, the object will come to a stop at $\bar{x}_2 = 2-\bar{x}_1$ after time $\bar{t} = \pi$. Now, if $|\bar{x}_2| \le 1$, once again the object stops. Since we assumed that $\bar{x}_1 > 1$, this condition can only be satisfied if $\bar{x}_2$ is negative, so we will have $|\bar{x}_2| = \bar{x}_1 - 2$. Repeating this process until $|\bar{x}| < 1$, we have $$ |\bar{x}_{n+1}| = |\bar{x}_{n}|-2\Longrightarrow |\bar{x}_{n+1}| = \bar{x}_1 - 2n, $$ and $|\bar{x}| < 1$ will occur after $\lfloor \bar{x}_1/ 2\rceil$ iterations, where $\lfloor x\rceil$ is the nearest integer function. Thus, the time $T$ and location $X_f$ where an object at rest at initial position $x_1$ comes to a stop will be $$ \omega T(x_1) = \pi\left\lfloor \frac{ x_1}{ 2L}\right\rceil\;\;,\;\;X_f(x_1) = (-1)^{\left\lfloor x_1/(2L)\right\rceil}\left(x_1 - 2L\left\lfloor \frac{x_1}{ 2L}\right\rceil\right) $$ If the initial condition has nonzero velocity, we need to find where and how long it takes to turn around. Calling these initial conditions $x_0$ and $v_0$ and again assuming $x_0 > 0$, the solution is $$ \operatorname{sgn}(v_0)\bar{x}(\bar{t}) = 1 - [\operatorname{sgn}(v_0)\bar{x}_0 - 1]\cos(\bar{t}) + |\bar{v}_0|\sin(\bar{t}), $$ where $\bar{v}_0 = v_0/(L\omega)$. Finding the turning point of this thing involves some fun algebra, but I'll skip the details and jump to $$ T_0(x_0,v_0) = \frac{1}{\omega}\tan^{-1}\left[\frac{\bar{v}_0}{\bar{x}_0\pm 1}\right]\;\;\;,\;\;\; x_1(x_0,v_0) = L\left[\sqrt{\bar{v}_0^2+(\bar{x}_0 \pm 1)^2} \mp 1\right] $$ Putting this all together gives \begin{eqnarray} T(x_0,v_0) &=& \frac{1}{\omega}\tan^{-1}\left(\frac{\bar{v}_0}{\bar{x}_0\pm 1}\right)+\frac{\pi}{\omega}\left\lfloor \frac{\bar{x}_1}{2}\right\rceil\\ X_f(x_0,v_0) &=& (-1)^{\lfloor x_1/2\rceil}\left(\bar{x}_1 - 2\left\lfloor \frac{\bar{x}_1}{2}\right\rceil\right)L \end{eqnarray} where again, $\omega = \sqrt{k/m}$, $L = \mu g/\omega^2$, $\bar{x}_0 = x_0/L$, $\bar{v}_0 = v_0/(L\omega)$, $v_0 = \pm |v_0|$, and $$ \bar{x}_1 = \sqrt{\bar{v}_0^2+(\bar{x}_0 \pm 1)^2} \mp 1. $$ The numerical checks I did seem to agree with this. There's another way to do this using the total mechanical energy $\bar{x}'^2 + \bar{x}^2$ which gave me the same answer, but keeping track of all the signs in that gets annoying fast.

  • $\begingroup$ Thanks for taking the time to answer. Unfortunately I have to travel for work and I cannot review it yet in detail. Which I notice so far is that some edition have being made and the first definitions of the variables with a bar have being lost, It was intentional? Since in the current form it doesn't show from where you are making appear those variables. Also, since you jump some steps, Did you solve it or you are using some source? If any source is used I ask you to share it to look for details, since I highly interested in figure out why it works, but my background in math is a bit limited $\endgroup$
    – Joako
    Commented Feb 12, 2023 at 2:48
  • 2
    $\begingroup$ @Joako Definitions now restored. This was all my own work. $\endgroup$ Commented Feb 13, 2023 at 17:20
  • $\begingroup$ Hi. I have started recently to trying to reproduce your result to understand what you did, and I have a question about the change of variables: You have assume that is true something of the form $$\text{sgn}\left(\frac{\mu g \sqrt{km}}{k}\bar{x}'\right) = \text{sgn}(\bar{x}')$$ due $\mu>0,\,g>0,\,k>0,\,m>0$, right? $\endgroup$
    – Joako
    Commented Mar 23, 2023 at 11:46
  • $\begingroup$ I have another issue with your analisis: as in other ansatz I have been trying using the solution to $y''+y+1=0$ like $$x(t)=c_1\sin(t)+c_2\cos(t)-\text{sgn}(c_1\cos(t)-c_2\sin(t))$$ for solving $x''+x+\text{sgn}(x')=0$, is that they are pure trigonometric functions, without the decay expected for a spring experiment which can be seen in numerical examples like this one which use initial conditions $x(0)=3\pi/2$ and $x'(0)=10$... any thoughts? $\endgroup$
    – Joako
    Commented Mar 24, 2023 at 21:25
  • $\begingroup$ Could you add a picture of a numerical example of your formula? (for $x''+x+\text{sgn}(x')=0$)... I want to see if it shows a decaying behavior in time (I still reviewing it, and I got lost within symbols). Thanks beforehand. $\endgroup$
    – Joako
    Commented Mar 24, 2023 at 22:53

The following phase-space analysis confirms many of the results of eyeball frog using simple geometrical principles. (Note however that in contrast to eyeball frog's analysis, I do not assume that ALL physical parameters in the problem can be scaled to unity. The model below still depends on one parameter $c$. The qualitative properties of the system are essentially the same for all choices of $c$ however, so eyeball frog's analysis in which $c=1$ is qualitatively correct and representative of all other cases. )

After scaling variables one can write the governing equation as $$\ddot x+ x = -c \ sgn(\dot x)$$ where $c$ is a friction constant.

Below we illustrate the system behavior in a particular numerical case ($c=.8$). The diagram below is called a phase-space diagram. The phase space coordinates are $x$= position, $y =\dot x $=velocity. The diagram suggests that some solutions get stuck quickly, and others oscillate a bit and then get stuck. We will now confirm this conjecture with a quantitative analysis that fully describes the behavior of the system throughout each oscillation.

(i) Consider an oscillation that begins at the far left axis with initial data $y=0$ and $x_0<0$. (The trajectory in phase space with these initial conditions travels clockwise with $x$ increasing.) Note that in any time interval during which $\dot x>0$, we may make the substitution $u= x+c$ into the governing equation and deduce that $u$ satisfies the simple harmonic oscillator equation $\ddot u + u =0$. This portion of the trajectory is motion on a semi-circular arc in phase space whose center is $x=-c$. This tells us that $u^2+ \dot u^2$ is conserved during this portion of the trajectory. Moreover, the analytical form of the solution is $ u= u_0 \sin t, \dot u= u_0\cos t$ where $u_0= x_0+c$. Returning to $x$ coordinates, we deduce that the analytical solution is $x(t)= (x_0+c) \sin t -c$ and $ \dot x= (x_0+c) \cos t$. This solution is valid until the trajectory crosses the axis defined by $y=0$. Thus we proceed to the next stage.

(ii) Once the trajectory crosses the positive $x$ axis at some new value $x_1$, a similar calculation can be performed. When $x$ is decreasing, the trajectory travels on a semi-circular arc in phase space centered at $x=c$. Note that this center differs from the center described in (i).

P.S. It is useful to note that the system will dissipate energy with every oscillation. Integrating e.g. over a time period where $x$ is increasing from say $x_1$ to $x_2$, we deduce that $$(2) \qquad \frac{1}{2}[ (\dot x_2)^2 + x_2^2)-(\dot x_1^2+ x_1^2)] = -c (x_2- x_1)<0$$ (The same argument works for the next half-cycle.)

Thus the total energy ( =potential + kinetic) $=E$= $\frac{ 1}{2} ( x^2+ \dot x^2)$ declines because of friction. That is, the trajectories move closer to the origin in phase space as time increases.

Summary. If $x$ is initially in the far fourth quadrant in phase space it can travel on a wide arc into the first quadrant, then the second, third, fourth, then complete a tight arc in the fourth quadrant and stop. On the other hand, with less energetic initial conditions, fewer oscillations occur before the solution hits the horizontal axis.

The motion on the $x$ axis is simply a sequence of to-and fro oscillations that diminish in amplitude and then halt.

In effect, the line segment in phase space between $ -c<x<c$ is a wall that the trajectories swirl around and eventually hit.

P.P.S. As a mechanical analogy, one can imagine a spring that has one end pinned and one end free, and with each half-cycle, as the spring pauses,the location of the pin jumps instantaneously. The jumps made by the pin get smaller in equal increments.

enter image description here

  • $\begingroup$ Thanks for taking the time to answer. But I don't really get your point with this answer, Are you wanting to pinpoint that a non-decaying solution is mistaken? (as the purely trigonometric solution founded by @eyeballfrog ? - I have now this opinion that it must be decaying) $\endgroup$
    – Joako
    Commented Mar 30, 2023 at 21:32
  • $\begingroup$ Exactly. The solutions lose energy (due to friction) with each half-oscillation until they come to a halt on the x axis, where y= velocity =0. You can vary the friction parameter and make other plots using WolfAlpha Streamplot $\endgroup$
    – MathFont
    Commented Mar 30, 2023 at 23:04
  • $\begingroup$ Thanks. I got interested exactly because of this, that it stops moving: an ODE must be nonlinear with at least one Non-Lipschitz point where it becomes singular in order to be possible to hold solutions with a finite extinction time. Differently from traditional ODEs which solutions are Power Series, which due the Identity Theorem can at best have solutions that vanishes at infinity, so, never-ending in time. Meaning this, that at least in what accuracy could say, there are 2 different kind of phenomena: those who stop moving and those who not, but today we model everything as the last ones. $\endgroup$
    – Joako
    Commented Mar 31, 2023 at 2:05
  1. Let us for convenience go to dimensionless barred coordinates for position, velocity and time, $$\begin{align} x ~=~& \frac{\mu g}{\omega^2} \bar{x}, \cr v ~=~& \frac{\mu g}{\omega} \bar{v},\cr t ~=~& \omega^{-1} \bar{t},\end{align}\tag{1}$$ cf. eyeballfrog's answer. Here we have introduced the undamped oscillator frequency $$ \omega ~\equiv~ \sqrt{\frac{k}{m}}. \tag{2}$$ We will not write the bars explicitly from now on. This has the effect of scaling the parameters $$ \mu g~=~1~=~\omega.\tag{3}$$ (The dimensionful parameters $\mu g$ and $\omega$ can easily be recovered in the end by dimensional analysis.)

  2. Then OP's 2nd-order ODE can be written as a system of two 1st-order ODEs $$\begin{align} \dot{x}~=~&v,\cr \dot{v}~=~&-x-{\rm sgn}(v), \end{align}\tag{4}$$ in an $(x,v)$ phase space.

  3. OP's system has a standard Hamiltonian formulation for $v\neq 0$. The Hamiltonian is (up to an arbitrary unphysical choice of energy-level) $$\begin{align} H(x,v)~=~&\frac{v^2}{2}+\frac{x^2}{2}+x{\rm sgn}(v)\cr ~=~&\frac{v^2}{2}+\frac{(x+{\rm sgn}(v))^2}{2}-\frac{1}{4}, \qquad v~\neq~ 0. \end{align}\tag{5}$$ In each half-plane of the $(x,v)$ phase space, the energy (5) is conserved, so the trajectories are semicircles with center $(-{\rm sgn}(v),0)$, cf. MathWonk's answer.

  4. The notion of energy (5) is not conserved when a trajectory crosses the $x$-axis in the the $(x,v)$ phase space. As long as the amplitude $|x|\geq 1$, we can continuously extend the solution from one semicircle to the next, and make the amplitude $|x|$ smaller and smaller. The dissipation of energy happens because every time we cross the $x$-axis we change the notion of energy (5).

  5. For an amplitude $|x|<1$, we can not extend the semicircle solution to the $x$-axis. Mathematically, the solution does not exist for initial values $|x|<1$ and $v=0$. (It is tempting to suggest that physically the particle is at rest at a constant position, but mathematically this is ill-defined.) The problem is that the ODE (4) is not Lipschitz continuous, and hence we cannot rely on the Picard-Lindelöf theorem to guarantee the existence and uniqueness of solutions.

  6. The discontinous ${\rm sgn}$ step function in OP's ODE is physically unrealistic$^1$. Let us regularize$^2$ OP's OPE (4) to instead $$\begin{align} \dot{x}~=~&v,\cr \dot{v}~=~&-x-f(v), \end{align}\tag{6}$$ where we have introduced a regularization function $$ f(v)~=~\begin{Bmatrix} 1 &{\rm for}& v \geq b^{-1} \cr bv &{\rm for}& |v|\leq b^{-1} \cr -1 &{\rm for}& v \leq -b^{-1} \end{Bmatrix}~=~\max(\min(bv,1),-1).\tag{7} $$ Here the regularization parameter $b\gg 1$ is the slope in Fig. 1.

    $\uparrow$ Fig. 1: The regularization function $v\mapsto f(v)$ with regularization parameter $b=3$.

  7. The ODE (6) is Lipschitz continuous, so existence and uniqueness of solutions hold. Effectively, we have introduced a linear drag in a narrow horizontal strip $|v|\leq b^{-1}$ in the $(x,v)$ phase space. This is just an overdamped oscillator$^3$ with well-known exponentially decaying solutions, cf. Fig. 2.

    $\uparrow$ Fig. 2: The $(x,v)$ phase space for the ODE (6) with regularization parameter $b=3$.

    The regularized model (6) is unaltered outside the strip, so the trajectories outside are still semicircles with center $(-{\rm sgn}(v),0)$. The amplitude $|x|$ decays towards the origin linearly (exponentially) in time $t$ for $|x|\geq 1$ ($|x|<1$), respectively.

    The regularized model (6) may still not be physically realistic, but at least it is mathematically well-defined.


$^1$ Ideally from a physical perspective, one should try to formulate a more realistic model of friction than just Coulomb friction. This is beyond the scope of this answer.

$^2$ From a mathematical perspective, there is a huge freedom in the choice of regularization. Our choice was partly influenced by a wish to obtain well-behaved closed-form piecewise analytic solutions.

$^3$ The 1D damped oscillator is dissipative but remarkably it still has a Hamiltonian formulation, cf. section 4 in my Phys.SE answer here.

  • $\begingroup$ Thanks for your answer. I would like to know Why you think this solution is somehow non-physical? I am really interested into this. I believe they are actually an accurate model than classic ones since without a Non-Lipschitz point in time an ODE can only have solutions that are never-ending in time, and series expansions cannot do the trick neither since cannot match a constant value. (...) $\endgroup$
    – Joako
    Commented May 15, 2023 at 20:42
  • $\begingroup$ (...) Instead, after the comments in this question, I think now that a better model for 2nd order systems (at least near their stopping time), will be something like: $$x''+a\ x'+b\ x+c+d\ \text{sgn}(x')=0$$ since it could stop moving, even its possible to stop in a place outside equilibrium, and have closed form solutions in the same way @eyeballfrog solve the main question in his answer. My intuition tells me that having never-ending solutions for a pendulum should be the non-physical model. Hope you could elaborate into this. $\endgroup$
    – Joako
    Commented May 15, 2023 at 20:45

$\color{green}{\textbf{The task standing.}}$

Let $\;P=\dot x = \dfrac{\text dx}{\text dt}= P(x(t)),\;$ then $$\ddot x=\dot P=\dfrac{\text dP}{\text dt}=\dfrac{\text dP}{\text dx} \dfrac{\text dx}{\text dt} = P(x)P'(x).$$ Then the given equation $$\ddot x+ mg \operatorname{sgn}\dot x +\dfrac km x=0\tag{$\Diamond$}$$ takes a form of $$PP'+ m g\operatorname{sgn}P=-\dfrac km x.\tag1$$

$\color{green}{\textbf{Solution for the negative derivative.}}$

If $P\le 0,\;$ then $$\dfrac 12(P^2)'=mg-\dfrac kmx,\quad P^2=2mgx-\dfrac km x^2 +W_0,$$ where $\,W_0\;$ is the starting energetic constant, which defines the amplitude of the oscilations. $$P=\dfrac{dx}{dt}=-\sqrt{W_0+2mgx-\dfrac km x^2},\tag2$$ Denote $\;a_0=\dfrac{mW_0}k, \quad b=\dfrac{m^2g}k>0,\quad \omega=\sqrt{\dfrac mk},\quad$ then $$\dfrac{\text dx}{\text dt}=-\dfrac1{\omega}\;\sqrt{a_0+2bx-x^2},\quad \omega\,\text dt= -\dfrac{\text dx}{\sqrt{a_0+2b x-x^2}},\quad \omega(t-t_0)=\arctan \dfrac{b-x}{\sqrt{a_0+2bx-x^2}}\ge 0, $$ $$\sec^2(t)=1+\tan^2\omega(t-t_0) = 1+\dfrac{(b-x)^2}{a_0+b^2-(b-x)^2} =\dfrac{a_0+b^2}{a_0+b^2-(b-x)^2},$$ $$(b-x)^2 = (a_0+b^2)-(a_0+b^2)\cos^2\omega(t-t_0)=(a_0+b^2)\sin^2\omega(t-t_0),$$ $$x(t)=b+s\sqrt{a_0+b^2}\sin\omega(t-t_0),\quad s=\pm1,$$ $$\dot x(t)=s\sqrt{a_0+b^2}\cos\omega(t-t_0)\le 0,\quad s=-\operatorname{sgn}(\cos\omega(t-t_0))$$ $$x(t)=b-\sqrt{a_0+b^2}\sin\omega(t-t_0)\operatorname{sgn}(\cos\omega(t-t_0)).\tag3$$

Solution $(3)$ is valid, until $$\operatorname{sgn}(\cos\omega(t-t_0))=\,\text{constant}.\tag4$$

$\color{green}{\textbf{Solution for the positive derivative.}}$

If $P\ge0,\;$ then from $(1)$ $$\dfrac 12(P^2)'=-mg-\dfrac kmx,\quad P^2=-2mgx-\dfrac km x^2 +W_1,$$ where $\,W_1\;$ is the starting energetic component, which defines the amplitude of the finite oscilations. If $\;W_1=0,\;$ then moving is not possible.

$$P=\sqrt{W_1-2mgx-\dfrac km x^2},\tag5$$ Denote $\;a_1=\dfrac{mW_1}k, \quad b=\dfrac{m^2g}k\ge 0,\quad \omega=\sqrt{\dfrac mk},\;$ then $$\dfrac{\text dx}{\text dt}=\dfrac1{\omega}\;\sqrt{a_1-2bx-x^2},\quad \omega\,\text dt= \dfrac{\text dx}{\sqrt{a_1-2b x-x^2}},\quad \omega(t-t_1)=\arctan \dfrac{b+x}{\sqrt{a_1-2bx-x^2}}, $$ $$1+\tan^2\omega(t-t_1) = 1+\dfrac{(b+x)^2}{a_1-b^2-(b+x)^2} =\dfrac{a_1-b^2}{a_1-b^2-(b+x)^2},$$ $$(b+x)^2 =a_1-b^2-(a_1-b^2)\cos^2\omega(t-t_1)=(a_1-b^2)\sin^2\omega(t-t_1),$$ $$x=-b+s\sqrt{a_1-b^2}\sin\omega(t-t_1),$$ $$\dot x=s\sqrt{a_1-b^2}\cos\omega(t-t_1)\ge 0,\quad s=\operatorname{sgn}(\cos\omega(t-t_0)),$$ $$x=-b+\sqrt{a_1-b^2}\sin\omega(t-t_1)\operatorname{sgn}(\cos\omega(t-t_0)).\tag6$$

Solution $(6)$ is valid under condition $(4).$

$\color{green}{\textbf{Example of trajectory.}}$

Taking in account $(3),(4),(6),$ easily to get the common solution in the form of $$x(t)=\left\{\begin{align} &b-\sqrt{a_0+b^2}\sin\omega(t-t_0)\operatorname{sgn}(\cos\omega(t-t_0)), \quad\text{if}\quad \dot x(t)\le 0\\ &-b+\sqrt{a_1-b^2}\sin\omega(t-t_1)\operatorname{sgn}\cos(\omega(t-t_0)),\quad\text{if}\quad \dot x(t)\ge 0\\[2pt] &a>0,\quad b>0. \end{align}\right. \tag7$$

From the OP task, based on the equality $\ddot x+\operatorname{sgn}\dot x+1=0,$ with the starting point $$x(0)=\dfrac32\pi,\quad \dot x(0)=10.$$ should $\omega=\sqrt{\dfrac km}=1,\quad b=mg=1,\quad \dot x(0)=-\sqrt{a-1} \,|\cos\omega(t-t_0))|=10.$

i.e. this data set looks not physically.

Let us consider the task $$\mathbf{\ddot x+6\operatorname{sgn} \dot x+1=0,\quad x(0)=6,\quad \dot x(0)= -10<0}.\tag 8$$ Assume $\quad m=k=1,\;$ then $\quad b=g=6,\quad p=\sqrt{a+36},\quad q=\sqrt{a-36},$

$$x(t)=\left\{\begin{align} &6-p\sin\omega(t-t_0)\operatorname{sgn}(\cos\omega(t-t_0)), \quad\text{if}\quad \dot x(t)\le 0\\[4pt] &-6+q\sin\omega(t-t_1)\operatorname{sgn}\cos(\omega(t-t_0)),\quad\text{if}\quad \dot x(t)\ge 0,\quad (p>6,\;p>q>0)). \end{align}\right. \tag{7a}$$

$$\dot x(t)=\left\{\begin{align} &-p\;|\cos\omega(t-t_0))|, \quad\text{if}\quad \dot x(t)\le 0\\[4pt] &q|\cos(\omega(t-t_0))|,\quad\text{if}\quad \dot x(t)\ge 0. \end{align}\right. \tag{7b}$$

If $\;t\in\left[0, \dfrac\pi2\right],\quad $ then $$\dot x_1(t)\le 0,\quad x_1(t)=6-p_1\sin(t),\quad \dot x_1(0)=-p_1=-10,\quad p_1=10.$$

Therefore, $\quad x_1(t)=6-10\sin(t),\quad t\in\left[0,\dfrac\pi2\right],$

Conjugational conditions for the next piece of the trajectory should be $$x_2\left(\dfrac\pi2\right)=x_1\left(\dfrac\pi2\right)=-4,\quad \dot x_2\left(\dfrac\pi2\right) = \dot x_1\left(\dfrac\pi2\right) = 0.$$

If $\;t\in\left[\dfrac\pi2,\dfrac32\pi\right],$ then $$\quad \dot x_2(t)\ge 0,\quad x_2(t)=-6-q_2\sin t,\quad x_2\left(\dfrac\pi2\right)= -6+q=-4,\quad q_2=2.$$

Therefore, $x_2(t)=-6+2\sin(t),\quad t\in\left[\dfrac\pi2,\dfrac32\pi\right].$

Conjugational conditions for the next piece of the trajectory should be $$x_3\left(\dfrac32\pi\right)=x_2\left(\dfrac32\pi\right)=-8.$$

If $\;t\in\left[\dfrac32\pi, \dfrac52\pi\right]\quad $ then $$\dot x_3(t)\le 0,\quad x_3(t)=6+p_3\sin(t),\quad \dot x_3\left(\dfrac32\pi\right)=6-p_3=--8,\quad p_3=14.$$

Therefore, $\quad x_3(t)=6+14\sin(t),\quad t\in\left[\dfrac32\pi,\dfrac52\pi\right],$

Conjugational conditions for the next piece of the trajectory should be $$x_4\left(\dfrac52\pi\right)=x_3\left(\dfrac52\pi\right)=20.$$

If $\;t\in\left[\dfrac52\pi, \dfrac72\pi\right],$ then $$\quad \dot x_4(t)\ge 0,\quad x_4(t)=-6+q_4\sin t,\quad x_4\left(\dfrac32\pi\right)= -6+q_4=20,\quad q_4=26.$$

Therefore, $x_4(t)=-6+26\sin(t),\quad t\in\left[\dfrac52\pi,\dfrac72\pi\right].$

Conjugational conditions for the next piece of the trajectory should be $$x_5\left(\dfrac72\pi\right)=x_4\left(\dfrac72\pi\right)=-32.$$

If $\;t\in\left[\dfrac72\pi, \dfrac92\pi\right],$ then $$\quad \dot x_5(t)\le 0,\quad x_5(t)=6-p_5\sin t,\quad x_5\left(\dfrac72\pi\right)= 6+p_5=-32,\quad p_5=38.$$

Therefore, $x_5(t)=6+38\sin(t),\quad t\in\left[\dfrac72\pi,\dfrac92\pi\right].$

Starting pary of the obtained trajectory is shown at the picture below.

Trajectory Plot


Let's ignore the extraneous constants, since the other answer shows how to eliminate them through rescaling. You have two different equations of motion: $$ \begin{eqnarray} x''(t) &=& -x(t) - 1 &\qquad& {\text{when }} x'(t) > 0;\\ x''(t) &=& -x(t) + 1 &\qquad& {\text{when }} x'(t) < 0.\\ \end{eqnarray} $$ Each has simple trigonometric solutions of the same form: $$x(t)=c_1\cos t + c_2\sin t \pm 1 = x(0)+c_1(\cos t - 1)+c_2\sin t.$$ Switching from one to the other occurs when $x'(t)=0$, so we only need to consider initial conditions where $x'(0)=0$. Since this implies $c_2=0$, in fact $$ x''(t) = -c_1\cos t=-x(t)+x(0)-c_1=-x(t)\mp 1, $$ from which we deduce that $c_1 = x(0) \pm 1$. We find $$ x(t) = x(0) + (x(0) \pm 1)(\cos t - 1) =x(0)\cos t \mp (1-\cos t) $$ for each interval over which the velocity is nonzero (taking $t=0$ to be the start of the interval). The velocity goes to zero again when $t=\pi$, at which time $x(\pi)=-x(0)\mp 2.$ And then iterate. So there's no decay, as such; there are just half-sinusoids with different amplitudes (and the same period) stitched together. The motion is smooth except at the turning points, where the second derivative has a jump discontinuity of magnitude $2$. For an example, suppose $x(0)=10$ (so the initial force is in the negative direction). Then the first piece of the solution is $x(t)=9\cos t + 1$. At $t=\pi$, the velocity is zero again, at $x=-8$. The next piece of the solution is therefore $x(t+\pi)=-7\cos t -1$. At $t=2\pi$, the velocity is zero again, at $x=6$. Etc. The pendulum swings from $10$ to $-8$ to $6$ to $-4$ to $2$ to $0$, at which point it stops. (It will stop whenever the velocity hits zero in the range $|x|\le 1$, since friction is then greater than the restoring force.)

Note that we can compare this to the numerical solution plotted in the question, which started at $x=10$ with $x'=-4$. Let's say that the true start was something like $x \approx 11$ with $x'=0$ at $t \approx -1/2$. (We could determine this exactly if required.) Then we expect turning points at $$x\approx -9,7,-5,3,-1,$$ at times $$t\approx \pi-1/2, 2\pi-1/2,3\pi-1/2,4\pi-1/2,5\pi-1/2 \\ \approx 2.6, 5.8, 8.9, 12.1, 15.2,$$ and a stop at the final turning point (that is, at $x \approx -1$). This describes the figure quite well.

  • 1
    $\begingroup$ Don't rely on your eyes :). Looking at the arrays generated by this command, the first minimum is $-8.8477$ at $t=2.7437$. The next relative maximum is $6.8479$ at $t=5.8663$. The next relative minimum is $-4.8477$ at $t=9.0050$. The next relative maximum is $2.8475$ at $t=12.150$. The final stopping point is at $-0.8474$ for $t \ge 15.289$. The amplitudes clearly behave exactly as predicted (decrease by $2$ each swing). The quarter-period durations are $3.1226,3.1387,3.1450,3.1390$, which are all $\pi$ to a fraction of a percent. $\endgroup$
    – mjqxxxx
    Commented Apr 6, 2023 at 18:45
  • 1
    $\begingroup$ I think what's counterintuitive here is that you might expect some "gradual reduction in the amplitude", because certainly if you look at the frictionless spring energy $(y')^2/2 + y^2/2$, that decreases continuously. But with a constant friction force, the individual lobes are not distorted from the original sinusoidal shape; the ODE solution over each quarter-period is still a simple sine wave (with a shifted center). You can consider using a nonlinear friction term like $-c \sqrt{|y'|}$ to get a result more like what you're expecting. $\endgroup$
    – mjqxxxx
    Commented Apr 6, 2023 at 18:53
  • 1
    $\begingroup$ Final note: I just looked at your most recent image and see that you have a key misunderstanding. The segments that are pure sinusoids are those between adjacent minima and maxima, not those between adjacent zeroes. The solution has a jump discontinuity in its second derivative at each minimum and maximum. $\endgroup$
    – mjqxxxx
    Commented Apr 6, 2023 at 18:58
  • 1
    $\begingroup$ No... because they are sinusoids centered at $\pm 1$, not at $0$. $\endgroup$
    – mjqxxxx
    Commented Apr 7, 2023 at 0:14
  • 1
    $\begingroup$ I believe the formulas in @eyeballfrog 's answer are correct... I haven't checked the derivation of the initial turning point, but it looks right. (In my answer I assumed $v_0=0$ for simplicity.) The envelope function will differ very slightly from being linear... a straight line passes through all the relative maxima or minima, but those lines aren't quite tangent to the curve. $\endgroup$
    – mjqxxxx
    Commented Apr 7, 2023 at 17:03

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .