Here's the problem:

For one of my classes, we're supposed to use Mathematica to solve the equation y''=y-b*y' for a variety of b values and specified boundary conditions and plot the outputs. There I have no problem and Mathematica runs well. However, the next problem is to replace the -b*y' with -b*sign(y') and repeat part one. When I do this, Mathematica runs for a long time, then gives up and outputs the input. I tried removing the boundary conditions and just doing the general solution, but it can't figure that out, either. It isn't giving any errors, it just doesn't work. Any ideas? This is what my input looks like:

In[5]:=  beta = 0.2;
         Solution = DSolve[{y''[x] == - y[x] - beta * Sign[y'[x]], y[0] == 1, y'[0] == -1}, y[x], x]
Out[6]= DSolve[{y''[x] == -0.2 Sign[y'[x]] - y[x], y[0] == 1, y'[0] == -1}, y[x], x]
  • $\begingroup$ Try NDSolve for {x, 0, 20}? $\endgroup$ – Michael E2 Sep 25 '14 at 18:52
  • $\begingroup$ Here's a question: does the homework specifically ask for an analytic solution to the boundary value problem, or does it only ask for the visual plots? If you're only interested in plotting the results, then replacing the symbolic DSolve with the numerical NDSolve will work better. DSolve is symbolic, and thus it may take a long time (or forever) to find a symbolic solution, whereas NDSolve is numerical, and thus is not (usually) subject to such complications. $\endgroup$ – DumpsterDoofus Sep 25 '14 at 19:13

Try to solve the problem numerically

beta = 0.2;

As you have observed, the problem is not solved in symbolic form:

Solution = 
 DSolve[{y''[x] == -y[x] - beta*Sign[y'[x]], y[0] == 1, y'[0] == -1}, y[x], 

(* -> DSolve[{(y^\[Prime]\[Prime])[x] == -0.2 Sign[Derivative[1][y][x]] - y[x], 
  y[0] == 1, Derivative[1][y][0] == -1}, y[x], x] *)

Now using NDSolve for x in the interval 0 to 10 (say):

yy[x_] = y[x] /. NDSolve[{y''[x] == -y[x] - beta*Sign[y'[x]], y[0] == 1, y'[0] == -1}, y[x], {x, 0, 10}] [[1]]

(* InterpolatingFunction[{{0.`,10.`}},"<>"][x] *)

This works well, and the plot looks good

Plot[yy[x], {x, 0, 10}]

(* 140925_plot_yy.jpg *)

enter image description here

The following part has been edited because of an error (26.09.14 00:52):

The solution lookes like a damped oscillation. We can try to derive the behaviour of the solution for large x as follows:

Multiplying the differential equation by y' and observing y' Sign(y') = Abs(y') we have

$y' y'' = - y y' - beta |y'|$

which can be written as

$\frac{1}{2} \frac{d}{dx} (y'^2 + y^2) = - beta |y'|$

Up to this point the reasoning is ok but in the following there is a flaw. Now, assuming $beta > 0$ the right hand side is always $<0$ which means that the expression $\frac{1}{2} (y'^2 + y^2)$ is only decreasing with x, which means that it will go to 0 which means in turn that $y\to 0$

Sorry, I shall clarify the asymptotic behaviour later.

Hope this helps, Regards Wolfgang


Edit 2: I think I've figured out how NDSolve constructs the solution, which has been incorporated in the revision below.


The term -beta Sign[y'[x]] is discontinuous when beta != 0, so the IVP is not guaranteed to have a solution. NDSolve does not indicate that, but it is obvious. In such cases NDSolve automatically processes the discontinuity (references given below). This confused me for a while, but now I think I've figured it out.

Second, the differential equation alternates between $$(y - \beta)'' = -(y - \beta) \quad\hbox{and}\quad (y + \beta)'' = -(y + \beta)$$ depending on the sign of $y'$. The solution to each is a sine curve with an axis of symmetry $y = \pm \beta$, respectively. When the differential equation has a solution, it will consist of a continuous sequence of half periods of sine curves. Each half period will cover an interval (of length $\pi$) where $y'$ does not change sign; therefore the curve itself will stretch between two consecutive extrema. Over consecutive segments, $y'$ changes sign and the axis of symmetry shifts by $2\beta$. This has the effect of subtracting $2\beta$ from the amplitude. This can continue until the amplitude is less than $2\beta$, which can happen for only finitely many steps. When the amplitude would be less than $2\beta$, there is no solution, that is to say, a solution can exist only over interval of the form $(-\infty, x_0)$. Edit 1: I should add that there is a special case in which a solution can be extended to the whole real line. When the amplitude of the sine curve is an odd multiple of $\beta$, the amplitude will step down until it equals $\beta$; at the end of the interval the solution will satisfy y(x_0) = y'(x_0) = 0 and the solution can be extended to $y(x) = 0$ for all $x > x_0$. (For an example, change the initial conditions below to ics = {y[0] == 6/5, y'[0] == 0} for beta = 0.2. The initial amplitude will be 5 beta because the axis of symmetry will be y == beta.)

On the other hand, NDSolve, which I suggested in a comment, creates a solution that becomes a nonzero constant after the amplitude falls below 2 beta. That seemed curious, especially since there is no message indicating that the solution might be wrong.

ode[beta_] := y''[x] == -y[x] - SetPrecision[beta, 100]*Sign[y'[x]];
ics = {y[0] == 1, y'[0] == -1};
residual[beta_] := ode[beta] /. Equal -> Subtract;

 Solution = NDSolve[{ode[beta], ics}, y, {x, 0, 20}, WorkingPrecision -> 20];
 Plot[{y[x], y'[x], -beta Sign[y'[x]]} /. First[Solution] // Evaluate,
  {x, 0, 20}, PlotRange -> All, 
  GridLines -> {ArcTan[-1 + beta, 1] + Pi Range[0, 6], beta Range[-5, 5]}, 
  PlotLegends -> {y[x], y'[x], -beta Sign[y'[x]]},
  Prolog -> {ColorData[97, 3], Dashed, 
    Line[Thread /@ {{#, y[First@#]}, {#, y[Last@#]}}] & /@ 
      Partition[ArcTan[-1 + beta, 1] + Pi Range[0, 5], 2, 1] /. First@Solution}],
 {{beta, 0.2}, 0, 2},
 TrackedSymbols :> Manipulate

Mathematica graphics

With WorkingPrecision -> MachinePrecision instead of arbitrary precision, there is some jumping around of the sign of y'[x], which is what I would expect if numerical integration is continued beyond the end point of the domain of the true solution. (When the sign of y'[x] changes and Abs[beta] > Abs[y[x]], then sign of y''[x] changes, and drives y'[x] back toward zero.)

Mathematica graphics

The residual from the ode tracks well up until ArcTan[-1 + 0.2, 1] + 3 Pi or 11.6703, at which point it is about 0.12, the same as the value of y[x] for x > 11.6703.

Block[{beta = 0.2},
 LogPlot[residual[beta] /. First@Solution // Abs,
   {x, 0, 20}, GridLines -> {ArcTan[-1 + beta, 1] + Pi Range[0, 6], None}]

Mathematica graphics

At any point there are three possibilities for Sign[y'[x]], namely -1, 0, 1. If the solution reaches a point where y'[x] == 0 and Abs[y[x]] < beta, then y'[x] could decrease, remain the same, or increase. Given that y''[x] has a value, there is only one choice, but it won't satisfy the ode. But the other possibilities won't satisfy the ode either, as shown by the residual for each choice not equaling zero:

residual[beta] /. Sign[y'[x]] -> {-1, 0, 1} /. First@Solution /. x -> 12.
(* {-0.319375, -0.119375, 0.0806248} *)

Why does NDSolve fail to give a warning?

Answer: NDSolve uses Filippov continuation automatically when discontinuous functions such as Sign appear explicitly in the differential equation. See the tutorial Events and Discontinuities in Differential Equations. This automatic processing of discontinuous functions can be turned off with the option Method -> {"DiscontinuityProcessing" -> False}.

NDSolve[{ode[0.2], ics}, y, {x, 0, 20}, 
  Method -> {"DiscontinuityProcessing" -> False}, 
  "ExtrapolationHandler" -> {Indeterminate &}];

NDSolve::mxst: Maximum number of 887537 steps reached at the point x == 11.677314840499118`. >>

In this case the phase vector field points in exactly opposite directions along y' == 0 for -beta < y < beta, so that the solution (and its continuation) is stationary:

Block[{beta = 0.2},
 StreamPlot[{p, -y - beta*Sign[p]}, {y, -2 beta, 2 beta}, {p, -2 beta, 2 beta},
  Axes -> True, Epilog -> {Red, PointSize[Large], Point[{{beta, 0}, {-beta, 0}}]}]

Mathematica graphics

From the first answer -- The following detects the problem faster. I suspect this is because the discontinuity occurs where y'[x] is zero, and therefore high Accuracy accuracy would be difficult to achieve with MachinePrecision.

One can get NDSolve to integrate up to ArcTan[-1 + 0.2, 1] + 3 Pi and give one of two warnings if AccuracyGoal is set to Infinity or to a high enough value:

NDSolve[{ode[0.2], ics}, y, {x, 0, 20}, AccuracyGoal -> Infinity]

NDSolve::ndtola: A component of the solution at x == 11.670315306734453` is essentially zero and the tolerance specified by the AccuracyGoal option could not be achieved. >>

NDSolve[{ode[0.2], ics}, y, {x, 0, 20}, AccuracyGoal -> 40]

NDSolve::ndsz: At x == 11.670315306734453`, step size is effectively zero; singularity or stiff system suspected. >>


It's a very interesting homework problem. If it's for an introductory undergraduate course in differential equations, I suspect a complete solution would not be expected. One of the problems with using computers to solve problems is that there is a tendency to take the answer and run with it. One sometimes forgets to think about the problem and the answer and to think of a way to check them. Perhaps that will be the "moral of story," so to speak, when your instructor reviews the assignment.

  • $\begingroup$ (at) Michael E2: I have outlined before in my solution how to proceed in general. I was stuck because I couldn't assess the asymptotic behaviour. I thought to have proved that y->0 for x->oo. But the numeric results looked different, as you have also confirmed afterwards. To be continued ... $\endgroup$ – Dr. Wolfgang Hintze Sep 26 '14 at 18:29
  • $\begingroup$ @Dr.WolfgangHintze Yes I suppose we both worked on it away from the site. I thought the conclusion that $y \rightarrow \infty$ was perhaps wrong, so I was waiting for you to finish your second answer before reading it. There are still issues to address, but I'm doubtful that it is worth the time (well, my time, at least). For instance, the solution above does not have a second derivative when the sign of y'[x] changes. If the solution is not twice-differentiable, it is not clear we can call it a solution. $\endgroup$ – Michael E2 Sep 26 '14 at 18:45
  • $\begingroup$ (at) Michael E2: I had some trouble understanding the contant asymptotic != 0 exhibited by the numeric solutions. First I got rid of the "explosion" you have observed by smoothing out the jump. Then I had the idea of the phase space consideration which led me to the calculation of the exact asymptotic value of y. My own argument of the never increasing distance of the system point from the origin in phase space turned out to be decisive in an unexpected way. Please read my answer #3 for details. $\endgroup$ – Dr. Wolfgang Hintze Sep 26 '14 at 22:19

Wolfgang's answer #3

It seems apropriate to create a new answer because the idea of my second answer is not taken up here.


We present here the complete analytic solution of the problem. First we consider the problem with the sharp step replaced by a smooth one. Then a study of the phase space trajectories of the system will give the key insight into the solution. Here we conclude that the system stops after a finite sumber of oscillations. This also solves the tricky problem of the asymptotic behaviour of y. Finally, the "time" dependence is calculated which completes the solution. All this is done for the case b=0.2 but the method is general.

Smoothing out the step

Let us take instead of Sign the smooth function ArcTan so that the diffenrential equation becomes

eq = y''[x] + y[x] == -b 2/\[Pi] ArcTan[p y'[x]]

(* y[x] + (y^\[Prime]\[Prime])[x] == -((2 b ArcTan[p Derivative[1][y][x]])/\[Pi]) *)

Here the parameter p controls the steepness of the jump.

With[{p = 5}, Plot[2/\[Pi]  ArcTan[p t], {t, -5, 5}]]
(* 140926_Smooth _Step.jpg *)

enter image description here

Solving the equation with NDSolve gives for increasing steepness

yy[x_] = With[{b = 0.2, p = 5}, 
   y[x] /. First[
     NDSolve[y''[x] + y[x] == -b 2/\[Pi] ArcTan[p y'[x]] && y[0] == 1 && 
       y'[0] == -1, y[x], {x, 0, 25}]]];

Plot[yy[x], {x, 0, 25}, PlotRange -> {-1.3, 1}, 
 PlotLabel -> 
  "y(x) calculated with NDSolve for smoothed-out step\nParameters: b = 0.2, p \
= 5", AxesLabel -> {"x", "y(x)"}]

(* 140926_y(x)_smooth _p _ 5.jpg *)

enter image description here

This lookes like a normal damped oscillation.

But now p = 50 leads to a different behaviour:

yy[x_] = With[{b = 0.2, p = 50}, 
   y[x] /. First[
     NDSolve[y''[x] + y[x] == -b 2/\[Pi] ArcTan[p y'[x]] && y[0] == 1 && 
       y'[0] == -1, y[x], {x, 0, 25}]]];

Plot[yy[x], {x, 0, 25}, PlotRange -> {-1.3, 1}, 
 PlotLabel -> 
  "y(x) calculated with NDSolve for smoothed-out step\nParameters: b = 0.2, p \
= 50", AxesLabel -> {"x", "y(x)"}]
(* 140926_y(x)_smooth _p _ 50.jpg *)

enter image description here

Now there are only three passes through zero. Afterwards y stays negative and is slowly approaching y = 0.

Letting now p=5*10^3

yy[x_] = With[{b = 0.2, p = 5*10^3}, 
   y[x] /. First[
     NDSolve[y''[x] + y[x] == -b 2/\[Pi] ArcTan[p y'[x]] && y[0] == 1 && 
       y'[0] == -1, y[x], {x, 0, 25}]]];

Plot[yy[x], {x, 0, 25}, PlotRange -> {-1.3, 1}, 
 PlotLabel -> 
  "y(x) calculated with NDSolve for smoothed-out step\nParameters: b = 0.2, p \
= 5,000", AxesLabel -> {"x", "y(x)"}]
(* 140926_y(x)_smooth _p _ 5k.jpg *)

enter image description here

Leads to a constant asymptotic values <0 for y.


(* -0.115909 *)

Comparing this with the sharp step

yy[x_] = With[{b = 0.2}, 
   y[x] /. First[
     NDSolve[y''[x] + y[x] == -b Sign[y'[x]] && y[0] == 1 && y'[0] == -1, 
      y[x], {x, 0, 25}]]];

Plot[yy[x], {x, 0, 25}, PlotRange -> {-1.3, 1}, 
 PlotLabel -> 
  "y(x) calculated with NDSolve for sharp step\nParameter: b = 0.2", 
 AxesLabel -> {"x", "y(x)"}]
(* 140926_y(x)_sharp.jpg *)

During evaluation of In[438]:= NDSolve::mxst: Maximum number of 10000 steps reached at the point x == 11.670368725029757`. >>

enter image description here

shows a breakdown of NDSolve as mentioned by Michael E2 but this happens at the same constant value of y which in the smoothed-out case is extended by NDSolve to large values of x.

Now this numeric result for the asymptotic behaviour needs to be verified independently because I have shown analytically earlier for the original Sign-function that

1/2 d/dx( y'^2 + y^2 ) = - b | y' |

Which means that the quantity y'^2 + y^2 cannot grow with increasing x. At first I thought it should go to zero together with |y'|, but this is not conclusive.

Can the equation of motion have a constant solution ya? The answer is yes: we have then y'' -> 0. This would ya -> - b Sign[y']. Now assume that for some x > x0, y' has a small positive value, i.e. Sign[y'] = +1 this leads to ya-> -b even if now y'->+0.

Now there is no contradiction: the left hand side 1/2 d/dx( y'^2 + y^2 ) -> 1/2 d/dx( 0 + y^2 ) = 1/2 d/dx( y1^2 ) = 0 and the right hand side -b |y'| -> 0 as well.

We shall make this more exact in the next section.

Phase space and asymptotics

Instead of pursuing my proposal with the piecewise trig functions (as in my soluton #2) let us look at the phase space (y, y'=v) of the problem.

The equation of motion is replaced by the system

dy/dx = v dv/dx = - y - b Sign(v)

The trajectories in phase space are given by dividing these two equations which leads to

dv/dy = - (y + b Sign(v))/v

For v<0 we have

dv/dy = - (y - b)/v

this can be integrated to give

v^2 + (y-b)^2 = const

This means the trajectories in the lower v-plane are circles centered at (y,v) = (b,0). The radius r1 is determined by the initial condition at x = 0: v(0) = -1, y(0) = 1, i.e.

r1 = Sqrt(1 + (1-b)^2)

Similarly, in the upper v-plane the trajectories are circles centered at (y,v) = (-b,0).

Let's draw a picture and look at the development of our example in phase space.

With[{b = .2},
 r1 = Sqrt[1 + (1 - b)^2]; w = ArcCos[(1 - b)/r1];
 pE = {r1 - 7 b, 0};
 pD = {-r1 + 5 b, 0};
 pC = {r1 - 3 b, 0};
 pB = {-r1 + b, 0};
 pA = {1, -1}; 
 Show[Graphics[{Line[{{-1.2, 0}, {1.2, 0}}], 
    Line[{{0, -1.5}, {0, 1}}], 
    Point[{0, 0}], {PointSize[Medium], Point[{b, 0}], 
     Point[{-b, 0}]}, {PointSize[Medium], RGBColor[1, 0, 0], 
     Point[pA], Point[pB], Point[pC], Point[pD], Point[pE]}, 
    Circle[{b, 0}, r1, {\[Pi], 2 \[Pi] - w}], 
      Circle[{-b, 0}, Abs[r1 - (4 k - 2) b], {0, \[Pi]}], {k, 1, 2}], 
     1], Flatten[
     Table[Circle[{b, 0}, Abs[r1 - 4 k b], {\[Pi], 2 \[Pi]}], {k, 1, 
       1}], 1], {Dashed, 
     Circle[{b, 0}, Abs[r1 - 8 b], {\[Pi], 2 \[Pi]}]}, 
    Text["A = {1,-1}", {1, -0.9}], Text["B", {-r1 + b - 0.1, 0.1}], 
    Text["C", {r1 - 3 b + 0.1, 0.1}], 
    Text["D", {-r1 + 5 b - 0.05, 0.1}], 
    Text["E", {r1 - 7 b + 0.05, 0.1}], Text["-b", {-b, 0.2}], 
    Text["b", {b, 0.2}], Text["v", {0, 1.1}], Text["y", {1.4, 0}]}]]]
(* 140926_Graph _phase-space.jpg *)

enter image description here

The system point P = (y,v) moves along successive circles of given center and radius as shown in the following table

$\left( \begin{array}{cccc} \text{segment} & \text{center} & \text{radius} & \text{y of end point} \\ A\to B & \{b,0\} & \text{r1=Sqrt[1+(1-b)${}^{\wedge}$2]} & \text{-r1+b} \\ B\to C & \{-b,0\} & \text{r2=r1-2b} & \text{r1-3b} \\ C\to D & \{b,0\} & \text{r3=r1-4b} & \text{-r1+5b} \\ D\to E & \{-b,0\} & \text{r4=r1-6b} & \text{r1-7b} \\ \end{array} \right)$

The journey stops at E = {r1-7b,0}. Why? Well, as we have shown, the distance of our system point P from the origin must never increase for increasing x. This condition is indeed met along the whole path from A to E, but would be violated if P would enter the dashed circle.

Hence the asymptotic value ya of y is reached after a finite time and is given by ya = r1 - 7b = Sqrt[1+(1-b)^2] - 7b /.b->0.2 = -0.119375. This value is now very close to the one observed with NDSolve.

I think this detailed study of the example showed all the ingredients necessary to find the asymptotic value of the solution. It would require just more labour to explore the behaviour for other values of b, especially the number of transits and hence the sign of the asymptotic value.

Also the "time" dependence of y, i.e. y(x), can be obtained easily from the phase space trajectory integrating the equation dy/dx = v giving x(y) = Integrate[1/v[y],y]. But this integral with an integrand of the type 1/Sqrt[ r^2 - (y \[PlusMinus] b)^2] is elementary and gives an ArcSin which leads to the time dependent solution which consists of segments:

y[k](x) = \[PlusMinus] b + r[k] Sin[x + c[k]]    
v[k](x) = r[k] Cos[x + c[k]]


x[k-1] <= x < x[k]

k = 1, 2, ..., kmax

Here r[k] are given in the table,

x[0] = 0, x[k] = Pi(k-1/2)+ArcSin[(1-b)/r[1]],

c[k] are constants determined by yk = y-value of the end point in the table.

kmax = 4 in our case.

I leave it for the interested reader to write down the explicit expressions and compare them with the previous results.

PS: I should remark that the procedure presented here improves and even replaces my solution attempt #2.

A surprisingly fascinating problem !

Regards, Wolfgang

  • $\begingroup$ 1. I like the way your phase diagram explains why the amplitude steps down by $2\beta$ over an interval of length $\pi$. It's easier to see than my explanation. 2. Somehow I'm not seeing how the smoothed equation illuminates the OP's problem. The solution for p = 5000 approaches zero, not a negative constant. (The eigenvalues of the linearization are $(-p \pm \sqrt{p^2-25\pi^2})/(5\pi)$, whose real components are negative.) 3. By "equation of motion", do you mean the OP's DE, or the one multiplied by $y'$? Note that if a DE is multiplied by $y'$, any constant trivially becomes a solution. $\endgroup$ – Michael E2 Sep 27 '14 at 19:20
  • $\begingroup$ ad 1. thanks. But the main result is that for the asymptotic value ya of y which is given by ya = Sqrt[1+(1-b)^2] - 7b /.b->0.2 = -0.119375. ad 2. first of all it seems a valid extension of the problem to be considered. Then, from the sequence of my graphs you can see that the asymptotic value is different from zero over an increasingly longer time as the steepness increases. This observation convinced me that ya !=0 which I then showed and calculated by means of the Phase space diagrams. ctd. in next comment. ad 3. I mean the original equation y''+ y = - b Sign[y']. $\endgroup$ – Dr. Wolfgang Hintze Sep 29 '14 at 6:42
  • $\begingroup$ ad 2. ctd: On the other hand, you are right, that for every finite p the asymptotic value is zero as can easily be seen from the approximation ArcTan[z] ~= z. Then y''+y = - 2/Pi b p y', and y = Exp[I w x] leads to w1 = I pb/Pi (1+Sqrt[pi^2/(pb)^2 -1], w2 = I pb/Pi (1-Sqrt[pi^2/(pb)^2 -1], and for p->oo: w1 ~= I 2pb/pi, w2 ~= I pi/(2pb). This means that from w2 y has an additive term Exp[-x pi/(2bp)] which decays exponentially the slower the larger p becomes. $\endgroup$ – Dr. Wolfgang Hintze Sep 29 '14 at 6:59
  • $\begingroup$ ad 2, ctd 2: Let me add that there is no difficulty with the approximation ArcTan(y' p) ~= y'p and p>>1. In fact, the approximation is valid for |y'p|<<1. But for the slowly decaying part, y~Exp[-x pi/(2bp)], we have |y'p| = |y p pi/(2bp)| = |y pi/(2b)| << 1 i.e. |y|<<(2/pi) b while p cancels out. $\endgroup$ – Dr. Wolfgang Hintze Sep 29 '14 at 16:47
  • $\begingroup$ I take it, however, that p >> 1 no longer "leads to a constant asymptotic values <0 for y," since we seem to agree that y -> 0, right? And if the equation of motion is the OP's $y''+y=-\beta \,\hbox{Sign}[y']$, then there are no nontrivial constant solutions, since a constant $y$ does not satisfy the differential equation ($\hbox{Sign}[y']=0$ for $y'=0$ in Mathematica). Do we agree about that? Now the smoothed DE is the explanation for Filippov's continuation method mentioned in the NDSolve documentation. Is that what you're talking about? $\endgroup$ – Michael E2 Sep 29 '14 at 17:07

Wolfgang's answer #4, 03.10.14

Reformulation of problem and solution

The following brief answer #4 presents finally the complete exact solution of the problem. It is based on all the intermediate results obtained so far, and therefore presents only the final results without derivations.

First of all, I have simplified the problem without loss of generality in two respects

a) put $b = 1$
b) "normalize" the initital conditions, see below

The exact solution of the DE (without losing generality we have put b = 1 in the original DE)

$y''(x) + y(x) = - sgn(y'(x))$

with the initial conditions (I call this "normalization", it can be done too without losing generality)

$y(0) = a, y'(0) = -0$ (so that $sgn(y'(0)) = -1)$

is given by the following function (it gives $y(x)$ as well as $y'(x)$)

y[x_, a_] := Module[{k, m, rv},
  k = 1 + Floor[x/\[Pi]];
  m = Ceiling[(a - 1)/2];
  rv = {If[
     k <= m, (-1)^(k + 1) + (a - 2 k + 1) Cos[
        x], (-1)^(m + 1) + (a - 2 m + 1) Cos[m \[Pi]]],
    If[k <= m, -(a - 2 k + 1) Sin[x], 0]}]

Plotting it with a = 4.9

With[{a = 4.9}, 
 Plot[{y[x, a][[1]], y[x, a][[2]]}, {x, 0, 4 \[Pi]}, 
  PlotRange -> {-1.1 a, 1.1 a}, 
  PlotLabel -> 
   "Exact solution of the DE (with Sign, normalized to b = 1)\nPosition y(x) \
(blue curve) and derivative y'(x) (red curve)\nInitial values y(0) = " <> 
    ToString[a] <> ", y'(0) = -0", AxesLabel -> {"x", "y(x),y'(x)"}]]
(* 141003_plot _y _yp _a _ 4_ 9.jpg *)

enter image description here

and with a = 6.9

With[{a = 6.9}, 
 Plot[{y[x, a][[1]], y[x, a][[2]]}, {x, 0, 4 \[Pi]}, 
  PlotRange -> {-1.1 a, 1.1 a}, 
  PlotLabel -> 
   "Exact solution of the DE (with Sign, normalized to b = 1)\nposition y(x) \
- blue curve, derivative y'(x) - red curve\nInitial values y(0)=" <> 
    ToString[a] <> ", y'(0)=-0", AxesLabel -> {"x", "y(x),y'(x)"}]]
(* 141003_plot _y _yp _a=_ 6_ 9.jpg *)

enter image description here

Comparison with the numerical solution whows excellent agreement.
In fact,

sol = With[{a = 6.9, \[Epsilon] = 10^(-4)}, 
    NDSolve[y''[x] + y[x] == -Sign[y'[x]] && (y[0] == 
        a) && (y'[0] == -\[Epsilon]), y[x], {x, 0, 9.4}]][[1]];
yn[x_] := y[x] /. sol

With[{a = 6.9}, 
 Plot[{y[x, a][[1]], yn[x] + 0.4}, {x, 0, 9.4}, PlotRange -> {-1.1 a, 1.1 a}, 
  PlotLabel -> 
   "Comparision of exact and numerical solutions\nof the DE (with Sign, \
normalized to b = 1)\nExact y(x) (blue curve) and numerical y(x) shifted by \
0.4 (red curve)\nInitial values y(0) = " <> ToString[a] <> 
    ", y'(0) = -0, -10^(-4) resp.", AxesLabel -> {"x", "y(x),y'(x)"}]]
(* 141003_plot _compare.jpg *)

enter image description here

The number of periods $m(a)$ performed by the system and the final point $yE(a)$ are given as a function of the parameter a by the following expressions

m[a_] := Ceiling[(a - 1)/2]

yE[a_] := Module[{m, rv},
  m = Ceiling[(a - 1)/2];
  rv = (-1)^(m + 1) + (a - 2 m + 1) Cos[m \[Pi]]]

Plot[{yE[a], m[a]}, {a, 0, 10}, 
 PlotLabel -> 
  "Exact solution of the DE (with Sign)\nNumber of periods (red curve) and \
final value of y(x) (blue curve)\nas a function of the initial value y=a \
(note that y'(0) = -0)", AxesLabel -> {"a", "yE(a), m(a)"}]
(* 141003_plot _m _yE.jpg *)

enter image description here

Concluding remarks

1) We have given the exact solution to the problem, i.e. give the dependent variable $y$(x) as a function of the independent variable $x$.
The System carries out a finite number of oscillations before coming to rest at some point $yE$ between -1 and +1.
2) The numerical solution breaks down close to the final point, but it is in excellent agreement with the exact solution up to that point.
3) It is paradoxial that the final value of the solution (which is reached in finite time) is not a solution to the DE. In fact $y''=0$ leads to $y = \[PlusMinus] 1$ rather than $yE$ which is smaller in absolute value.
This problem is still open and needs further clarification. I admit, this type of problem was new to me.
4) I very much appreciate the intense cooperation with Michael E2, and hope we can continue it to resolve the remaining open question.

Best regards, Wolfgang


Wolfgang's answer #2

The problem can be solved exactly up to the numerical solution of a goniometric equation (i.e. one involving Sin and Cos) for the zeroes of y'.

1. General theory

Indeed, the differential equation is either (letting $beta \to b$)

$y'' + y = +b \to y'<0 $


$y'' + y = -b \to y'>0 $

The solution is of the type

$y(x) = U + R*cos(x) + S*sin(x)$

with three constants $U$, $R$, and $S$ which will be defined below soon.

In the first interval the solution $y[1]$ is valid for $0<= x <= x1$ where $x1$ is the smallest positive root of $y1'(x) = 0$ and has the form

$y1(x) = b + R[1] cos(x) + S[1]sin(x)$

From the initial conditions we have explicitly

$y[1](0) = 1 = b + R[1]$
$y[1]'(0) = -1 = S[1]$


$R[1] = 1-b, S[1] = -1$

$x1$ is determined from $y[1]'(x1) = 0$ which leads to

$x1 = 2 arctan[1-b+Sqrt[b^2-2b+2]]$

In the second interval $x1 < x < x2$ we have $y'>0$ and therefore $y'' + y = -b$ with the solution

$y[2](x) = - b + R[2] cos(x) + S[2] sin(x)$

The requirement of continuity for $y(x)$ and $y'(x)$ at $x = x1$ gives

$y[2](x1) = y[1](x1) \to - b + R[2] cos(x1) + S[2] sin(x1) = + b + R[1] cos(x1) + S[1] sin(x1)$ $y[2]'(x1) = y[1]'(x1) \to - R[2] sin(x1) + S[2] cos(x1) = - R[1] sin(x1) + S[1] cos(x1)$

which gives

$R[2] = R[1] + 2 b cos(x1)$
$S[2] = S[1] + 2 b sin(x1)$

In the third interval $x2< x < x3$ we have $y' < 0$ and therefore $y'' + y = +b$ with the solution

$y[3](x) = + b + R[3] cos(x) + S[3] sin(x)$

The requirement of continuity for y(x) and y'(x) at x = x2 gives

$y[3](x2) = y[2](x2) \to + b + R[3] cos(x2) + S[3] sin(x2) = - b + R[2] cos(x2) + S[2] sin(x2)$
$y3'(x1) = y2'(x1) \to - R[3] sin(x2) + S[3] cos(x2) = - R[2] sin(x2) + S[2] cos(x2)$

which gives

$R[3] = R[2] - 2 b cos(x2) = R[1] + 2 b cos(x1) - 2 b cos(x2)$
$S[3] = S[2] - 2 b sin(x2) = S[1] + 2 b sin(x1) - 2 b sin(x2)$

And so on, giving for general k the solution

$y[k](x) = -(-1)^k b + R[k] cos(x) + S[k] sin(x)$

and the coefficients

$R[k] = R[1] - 2 b \sum _{m=1}^{k-1} (-1)^m cos(x[m])$
$S[k] = S[1] - 2 b \sum _{m=1}^{k-1} (-1)^m sin(x[m])$

Here $x[k]$ is the smallest root $> x[k-1]$ of $y[k]'(x) = 0$.

2. Implementation in Mathematica

to be done, later

3. Example

to be done, later

Regards, Wolfgang

  • $\begingroup$ I beg your pardon. At the moment I'm being lost in the phase space of this interesting problem. $\endgroup$ – Dr. Wolfgang Hintze Sep 26 '14 at 10:42

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.