WhenEvent&NDSolve: How to detect saddle point?

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!

• Can't tell if the result is the expected one, but what about WhenEvent[x'[t] == 0 && x''[t] == 0, sgn[t] -> -sgn[t]]? Jul 5, 2019 at 15:57
• Consider Plot[x''[t] /. sol, {t, 0, 25}] and the discrete change of sign in x''[t] (i.e. such that x''[t] != 0) that occurs around t~=21. It does not quite meet your definition of saddle. Jul 5, 2019 at 16:23
• @ C.E. Thanks, I tried the proposed Event, but NDSolve doesn't find a solution (MMA v 11.0.1) Jul 5, 2019 at 16:34
• @ Michael E2 Thanks, so it is perhaps an inflection point ! I tried an additional event WhenEvent[x''[t] == 0 && x'[t] == 0 , sgn[t] -> sgn[t] ] without success. Jul 5, 2019 at 16:39
• @UlrichNeumann I guess I was wondering if the discontinuous change of sign of x''[t] disqualified it as a saddle or not. I was guessing that you would include it, but I wasn't sure. Jul 5, 2019 at 16:42

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,
vars, {t, 0, 25}, DiscreteVariables -> Rest@vars];

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


• @ MichaelE2 Thank you for your answer, I have to analyse in detail! The saddle is well detected but sgn shouldn't change!. I changed the saddle -event to s[t]-> s[t] but doesn't work . Jul 5, 2019 at 17:04
• @UlrichNeumann The saddle event presumes that sgn[t] changes sign. If you do not change its sign, then x''[t] will not change sign and there will be an extremum, not a saddle. I'm not sure why you want sgn[t] not to change at a potential saddle, but if it's to prevent saddles in the solution, then it seems to me it's working (except my variables extrema and saddle no longer count extremas and saddles, which seems unimportant and/or fixable). Jul 5, 2019 at 17:24
• @ Michael E2 : sgn[t]~Sign[x'[t]] is what I'm looking for. I agree that x''[t] must change sign in the saddle point. In the plot of your answer I see a positive slope of the blue curve near the saddle. Jul 5, 2019 at 17:31
• @UlrichNeumann OK, I'll try to think about sgn[t]~Sign[x'[t]] -- discontinuities can be tricky. I don't think in every case there is a solution to an ODE with discontinuous coefficients. Plot x'[t] over a small range including the saddle: it decreases to zero (approximately) and then increases immediately with an angle like absolute value. A discrete plot should show a line segment with a positive slope, since an instantaneously zero slope is impossible to draw via line segments (which is what Plot[] does). Jul 5, 2019 at 17:38

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
]


• Thank you very much for your helpful comprehensive answer. With my attempt I tried to verify the solution of NDSolve of the ode with Sign (expecting a asymptotic solution x[t]==1). I need some time to elaborate your answer... Jul 6, 2019 at 13:09
• No prob. I think it could end up anywhere on that line segment from x=0.905 to 1.105. Jul 6, 2019 at 13:19
• That's the point where I'm struckling: If you check the residuum of the solution you'll see a value near to zero for t<21, but for t>21 you see nonvanishing switching constant values. Does that mean the numerical solution doesn't approximate the ode for t>21??? Jul 6, 2019 at 13:43
• Could you explain more about what you meant by the residuum of the solution? Jul 6, 2019 at 14:09
• @UlrichNeumann I believe it's the same as the problem here. I don't see how to use sgn to get around the issue with Sign. (I personally am convinced that Sign is the right way to go, and the linked problem gives the correct solution.) Jul 6, 2019 at 21:43