Closed-form solutions to $x''+\frac{k}{m}\ x+\mu\ g\ \text{sgn}(x')=0$
Introduction______________________
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:
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:
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];
opt=odeset('RelTol',1e-3,'AbsTol',1e-4);
[ts,ys] = ode45(f,[t0,20],y0,opt);
plot(ts,ys(:,1),'b');
Plot of the numerical solution to $x''+x+\text{sgn}(x')=0,\,x(0)=10,\,x'(0)=-4$:
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];
opt=odeset('RelTol',1e-3,'AbsTol',1e-4);
[ts,ys] = ode45(f,[t0,50],y0,opt);
plot(ts,ys(:,1),'b');
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.
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:
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):
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];
opt=odeset('RelTol',1e-3,'AbsTol',1e-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);
plot(ts,ys(:,1),'b',ts,x,'r');
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;
Best Answer
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.