4
$\begingroup$

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.

enter image description here

$\endgroup$
  • $\begingroup$ See Assumptions $\endgroup$ – Sascha Dec 28 '15 at 17:22
  • 1
    $\begingroup$ See this question on Mathematica.SE and especially the answer from Michael E2 $\endgroup$ – Sascha Dec 28 '15 at 17:30
  • $\begingroup$ As to the documentation on Assumptions and why DSolve is not listed supporting Assumptions: To my knowledge DSolve internally may use e.g. Simplify or Integrate which support Assumptions and this may solve your issue $\endgroup$ – Sascha Dec 28 '15 at 17:34
  • $\begingroup$ @Sascha: I used assumptions and it seems to have gotten further, with a new error message. Please see update. Thanks for that tip! $\endgroup$ – Moo Dec 28 '15 at 17:37
  • 1
    $\begingroup$ I added the message name. It should help others search for solutions when they get the same error. $\endgroup$ – Michael E2 Dec 28 '15 at 19:05
4
$\begingroup$

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}]}

Mathematica graphics

$\endgroup$
  • $\begingroup$ @MichealE2: Thanks for that wonderful answer. For some reason, your image is not showing. Also, is it possible to change eqns = {x'[t] == y[t], y'[t] == Piecewise[{{-x[t] + c Sign[x[t]], Abs[x[t]] > c}, {0, Abs[x[t]] <= c}}]}; ics = {x[0] == 1, y[0] == 1}; and do it numerically with various initial conditions? The numerical solution captured the essence of the phase plot. Lastly, is it possible to get StreamPlot to plot the solution using various values of $c$ in the previous update in my comment? I updated what the phase portrait would look like under these changes with $a$ being $c$ $\endgroup$ – Moo Dec 28 '15 at 19:42
  • 1
    $\begingroup$ @Variable (1) The image didn't work because I had to use HTML instead of Markdown. (2) You can use ParametricNDSolve. (3) StreamPlot does a lot of stuff automatically that sometimes drives me up the wall. How is this: StreamPlot[{y, ny[x, y]}, {x, -2, 2}, {y, -1, 1}, AspectRatio -> Automatic, StreamScale -> {1.5, All, 0.05}, StreamPoints -> Join[Table[{1, i/5}, {i, 5}], -Table[{1, i/5}, {i, 5}]]]? Sometimes, though, it's easier to do the kind of thing you want directly using NDSolve and ParametricPlot. See mathematica.stackexchange.com/q/60838 for plotting arrows. $\endgroup$ – Michael E2 Dec 28 '15 at 20:08
  • $\begingroup$ Right clicking my Magic-mouse opens up your ugly image in a new window or tab on my OS X system. No need for Ctrl clicking. $\endgroup$ – m_goldberg Dec 29 '15 at 2:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.