WhenEvent&NDSolve: How to detect saddle point?

Question

I'm considering a simple ode. With the help of WhenEvent and Discrete Variable it is possible to detect the sign of the solution x'[t] and switch the ode accordingly(avoiding Sign-function).

sol = NDSolveValue[{x''[t] == 1 - x[t] - .05 sgn[t] (1 + x[t]) ,
x[0] == 0, x'[0] == 1, sgn[0] == 1 ,WhenEvent[x'[t] == 0  , sgn[t] -> -sgn[t] ]}, {x, sgn}, {t, 0, 25}, DiscreteVariables -> {sgn[t] }] 

Plot[{sol[[1]][t], sol[[2]][t]}, {t, 0, 25}]

The sign-changes of x'[t]are detected very well until t~=21. Here the solution shows a saddle point.

My question: How to define two WhenEvents, which distinguish extrema(x'[t]==,x''[t]!=0) and saddle points(x'[t]==,x''[t]==0)?

Thanks!

Answer 1

It's a bit tricky: First, in WhenEvent[], only the "state" variables are substituted. This means you can't use the highest order derivative (maybe WRI should rethink that). One trick is to differentiate the ODE to raise the order, but that doesn't work here because the ODE is discontinuous. I resorted to using the "right-hand side" of the ode as a substitute for x''[t]. Second, getting the change of sign in x''[t] is tricky because it is discontinuous at the event x'[t] == 0. However, you can examine x''[t] on both sides of the event with rhs[x[t], ±sgn[t]]. From these two values, you can deduce whether the second derivative changes sign or not.

vars = {x, sgn, extrema, saddle};
rhs[x_, s_] := 1 - x - .05 s (1 + x);
sol = NDSolve[{x''[t] == rhs[x[t], sgn[t]],
    x[0] == 0, x'[0] == 1,
    sgn[0] == 1, extrema[0] == 0, saddle[0] == 0,
    WhenEvent[
     x'[t] == 0 && rhs[x[t], sgn[t]] rhs[x[t], -sgn[t]] > 0,
     {extrema[t] -> 1 + extrema[t], sgn[t] -> -sgn[t]}],
    WhenEvent[
     x'[t] == 0 && rhs[x[t], sgn[t]] rhs[x[t], -sgn[t]] <= 0,
     {saddle[t] -> 1 + saddle[t], sgn[t] -> -sgn[t]}]},
   vars, {t, 0, 25}, DiscreteVariables -> Rest@vars];

Plot[Through[vars[t]] /. First@sol // Evaluate, {t, 0, 25}, 
 PlotLegends -> vars]

Answer 2

Here are some thoughts building on @MichaelE2's answer.

I'm visually oriented, so I converted your system into two first-order ODEs (with v[t]==x'[t]) to look at this in the phase-plane.

sol = NDSolve[{
  v'[t] == 1 - x[t] - 0.05 sgn[t] (1 + x[t]),
  x'[t] == v[t],
  x[0] == 0, v[0] == 1, sgn[0] == 1, 
  WhenEvent[v[t] == 0, sgn[t] -> -sgn[t]]},
  {x, v, sgn}, {t, 0, 25}, DiscreteVariables -> {sgn[t]}][[1]];

First make sure nothing changed:

Plot[Evaluate[{x[t], sgn[t]} /. sol], {t, 0, 25}]

Now we can plot the x[t] vs v[t]=x'[t] phase plane. The tricky part is there are two sets of streams depending on sgn[t]. I made the sgn[t]==1 streams blue and the sgn[t]==-1 ones gray.

{xmin, xmax} = {-0.2, 2.4}; {vmin, vmax} = {-1.2, 1.4};
sp1 = StreamPlot[{v, 1 - x - 0.05 (1 + x)},
  {x, xmin, xmax}, {v, vmin, vmax},
  StreamStyle -> Lighter[Blue, 0.5], PlotRangePadding -> 0];
sp2 = StreamPlot[{v, 1 - x + 0.05 (1 + x)},
  {x, xmin, xmax}, {v, vmin, vmax}, StreamStyle -> Gray];
Show[sp1, sp2, Graphics[Line[{{xmin, 0}, {xmax, 0}}]], FrameLabel -> {"x", "v=x'"}]

Well, that's an eye-full, but look along the v==0 line and notice that the blue and gray arrows agree in direction on the outside but in the middle they are pointed at each other.

Now plot the numerical solution and color-code it based on sgn[t]:

pp = ParametricPlot[{x[t], v[t]} /. sol, {t, 0, 25}, 
  ColorFunction -> Function[{x, y, t}, If[(sgn[t] /. sol) > 0, Blue, Black]],
  ColorFunctionScaling -> False, PlotStyle -> Thick];
Show[sp1, sp2, pp, Graphics[Line[{{xmin, 0}, {xmax, 0}}]], FrameLabel -> {"x", "v=x'"}]

Zooming in on the interesting part:

{xmin, xmax} = {0.75, 1.2}; {vmin, vmax} = {-0.2, 0.2};
...

So I'd say that NDSolve is doing exactly what you asked: switch sgn[t] when x'[t]==0 and then carry on.

If you want one equation for v'[t] when Sign[v[t]]>0 and the other when Sign[v[t]]<0, why not just use Sign?

sol = NDSolve[{
     v'[t] == 1 - x[t] - .05 Sign[v[t]] (1 + x[t]),
     x'[t] == v[t],
     x[0] == 0, v[0] == 1}, {x, v}, {t, 0, 25}][[1]];

Plot[Evaluate[x[t] /. sol], {t, 0, 25}]

Here the solution gets trapped once it hits the region where the streams point at each other. We can find that region by solving where each vector field changes directions:

Solve[1 - x - 0.05 (1 + x) == 0, x]
Solve[1 - x + 0.05 (1 + x) == 0, x]
(* {{x -> 0.904762}} {{x -> 1.10526}} *)

Putting it together:

{xmin, xmax} = {0.7, 1.2}; {vmin, vmax} = {-0.2, 0.2};
Show[
  StreamPlot[{v, 1 - x - Sign[v] 0.05 (1 + x)},
  {x, xmin, xmax}, {v, vmin, vmax},
  StreamStyle -> Gray, StreamPoints -> Fine],
  ParametricPlot[{x[t], v[t]} /. sol, {t, 0, 25}, PlotStyle -> Red],
  Graphics[{Blue, Opacity[0.2], Rectangle[{xmin, 0}, {xmax, vmax}]}],
  Graphics[{Black, Thick, Line[{{0.904762, 0}, {1.10526, 0}}]}],
  FrameLabel -> {"x", "x'=v"}, 
  PlotRange -> {{xmin, xmax}, {vmin, vmax}}, PlotRangePadding -> 0
]