DSolve with Piecewise Function in System of DEQs

Question

I have been messing around with this problem from MSE, which is given by:

$$ \ddot{x} = \begin{cases} -x + c\cdot \operatorname{sgn}(x)& |x| > c\\ 0 & |x|\leq c \end{cases} $$ where $c > 0$ are some given constants. I can solve the entire problem by hand, but was curious how to do it in MMA.

For this discussion, let $c = 1$, but would like to do this for a general $c$ as it makes a difference in the solution. After writing the second-order DEQ as a system of first-order equations, we get:

  eqns = {x'[t] == y[t], 
          y'[t] == Piecewise[{
            {-x[t] + 1 Sign[x[t]], Abs[x[t]] > 1},
            {0, Abs[x[t]] <= 1}
           }],
          x[0] == 1, y[0] == 1};

  sol = Assuming[x[t] ∈ Reals && y[t] ∈ Reals, DSolve[eqns, {x, y}, t]]

I get the following error message and it hangs:

DSolve::bvnul For some branches of the general solution, the given boundary conditions lead to an empty solution.

Suggestions?

Note: One can manually find the solutions for the three cases as:

$$y^2 + (x \mp c)^2 = C_{1,2}, y = C_3$$

The phase plane would be semicircles that have straight lines (elongated by the setting of $c$).

Update

The phase portrait should look like the following (note the $a$ in the picture is $c$). They are not nicely connected like that.

Answer 1

Here's some headway. One should note at the start that second-order, nonlinear, discontinuous differential equations are hard to deal with symbolically. I'm not exactly sure where the current theory is at present, but Mathematica's solution to this equation is a beast to compute with. You may be much better off trying to handle the DE numerically with NDSolve.

First, your error message indicates that Mathematica can find the "general" solution. Whether this contains the particular solution is not guaranteed; see, for instance, DSolve misses a solution of a differential equation and DSolve not finding solution I expected.

So we should first get the general solution and try to solve for the initial conditions. We could try to set generic initial conditions ics = {x[0] == c1, y[0] == c2}, but that failed for me. So lets split the system like this:

eqns = {x'[t] == y[t], 
   y'[t] == 
    Piecewise[{{-x[t] + 1 Sign[x[t]], Abs[x[t]] > 1}, {0, Abs[x[t]] <= 1}}]};
ics = {x[0] == 1, y[0] == 1};

Now DSolve will rewrite the variables x[t] and y[t] as just x and y at some point in its computation. We can glimpse it here:

sol = DSolve[eqns, {x, y}, t];

Integrate::pwrl: Unable to prove that integration limits {1,x[t]} are real. Adding assumptions may help. >>

Integrate::pwrl: Unable to prove that integration limits {1,x} are real. Adding assumptions may help. >>

Integrate::pwrl: Unable to prove that integration limits {1,x} are real. Adding assumptions may help. >>

General::stop: Further output of Integrate::pwrl will be suppressed during this calculation. >>

(Note it does return a solution, but it is in terms of unevaluated integrals. We can do better.)

It looks like we should assume all forms of the variables to be real, just to be safe. There are two solutions, which agrees with the OP's implicit solution.

sol = Assuming[
   x[t] ∈ Reals && x ∈ Reals && y[t] ∈ Reals && y ∈ Reals, 
   DSolve[eqns, {x, y}, t]];
Length@sol
(*  2  *)

Now it turns out Solve balks at trying to solve the initial conditions. This should be no surprise, since that is basically what happened inside DSolve when it returned the DSolve::bvnul error. So let's try a numeric solution:

FindRoot[ics /. Last@sol /. {C[1] -> c1, C[2] -> c2}, {{c1, 1}, {c2, 1}}]
(*  {c1 -> -1.20162*10^-16, c2 -> 2.82843}  *)

They look like 0 and Sqrt[8] -- how lucky!

2.8284271247461907`^2
(*  8.  *)

Check:

ics /. Last@sol /. {C[1] -> 0, C[2] -> Sqrt[8]}
(*  {True, True}  *)

Whoopee...To see what you're up against, here's an image of the solution (CTRL-click to open image in a new window on a MAC; I suppose right-click might work on Windows):

---

Addendum: Example numeric approach

It does appear to be periodic, so we can add a WhenEvent to compute the period as well:

solN = NDSolve[{eqns, ics, 
    WhenEvent[x[t] == ics[[1, 2]] && Abs[y[t] - ics[[2, 2]]] < 10^-6, "StopIntegration"]},
   {x, y}, {t, 0, 20}];

The period appears to be 4 + 2 Pi:

x["Domain"] /. solN
(*  {{{0., 10.2832}}}  *)

Phase curve:

ParametricPlot @@ {{x[t], y[t]} /. solN, Flatten[{t, x["Domain"] /. solN}]}