Discontinuous ODE results in NDSolve::ndinid warning, initial condition is not in the range specified by the discrete variable (internal variable)

Question

The equations and code are below.

The short version is:
I know the error is because NDSolve is doing something in the background to handle the discontinuous ODE. Adding the following to NDSolve produces the correct result but it's incredibly ineffecient.

    Method -> "DiscontinuityProcessing" -> False, AccuracyGoal -> 6,
    PrecisionGoal -> 6, MaxSteps -> 10^7

Solving it this was is incredibly slow for a second order system. What's the right way to handle discontinuous ODEs?

The ODE: $$ \theta 1''(t) = \frac{-\theta 1'(t)^2 \sin(\theta 1(t)) \cos(\theta 1(t))+2 \sin(\theta 1(t))-u(t) \cos(\theta 1(t))}{2-\cos^2(\theta 1(t))} $$ where $$ u(t) = 2 \tan(\theta 1(t))-\theta 1'(t)^2 \sin(\theta 1(t)) + u_1(\theta1(t),\theta1'(t))(2 - \cos^2(\theta1(t))) $$ and $$ u_1(z_1,z_2) = \begin{cases}\begin{cases}1 & z_2 \leq -\sqrt{2z_1} \newline-1 & z_2 >-\sqrt{2z_1} \end{cases} & z_1 > 0 \newline\begin{cases}1 & z_2 < \sqrt{-2z_1}\newline-1 & z_2 \geq \sqrt{-2z_1} \end{cases} & z_1 \leq 0 \newline\end{cases} $$ For reference, the first two terms of the control eliminate the right hand side of $\theta1''$, and the last term is a time optimal control that drives the state to $(\theta1 = 0, \theta1'=0)$, which comes from the Pontryagin Maximum Principle.

I tried to use the following code to solve it

    sol = NDSolve[{eqs, θ1[0] == -Pi/6, θ1'[0] == -1},θ1[t],{t, 0, 10}];

which returns the warning

The reason the control has such an ugly piecewise definition is to check that we're never taking the square root of a negative number. So, it's not clear to me what is happening that's skipping all these cases and evaluating the square root anyway.

I see in the reference page for NDSolve there is a WhenEvent option that can be implemented, but it's not clear how to specify that the action should only be taken once (when the state reaches the switching curve, and not continue to cycle through the discrete variables for each step that it remains on the switching curve).

ClearAll[x1, x, y, s, u, eqs, u1, sol, uop];
l1 = 1;
m1 = 1;
mCart = 1;
g = 1;

x = s[t] + l1 Sin[θ1[t]];
y = l1 Cos[θ1[t]];

T = (1/2) m1 (D[x, t]^2 + D[y, t]^2) + (1/2) mCart (s'[t])^2;
V = m1 g ( y - 1);
L = T - V; 
ELE = Flatten[
  Solve[
    FullSimplify[D[D[L, {{θ1'[t], s'[t]}}], t]] - 
      FullSimplify[D[L, {{θ1[t], s[t]}}]] == {0, u[t]},
     {θ1''[t], s''[t]}
    ] /. {Rule -> Equal}
  ]
u1[t] = (u[t] /. Solve[ELE[[1]] /. {θ1''[t] -> 0}, u[t]] // 
    FullSimplify)[[1]]
optimalcontrol[z1_, z2_] =
 Piecewise[{
   {Piecewise[{
      {1, z2 <= -Sqrt[2 z1]},
      {-1, z2 > -Sqrt[2 z1]}}], z1 > 0},
   {Piecewise[{
      {1, z2 < Sqrt[-2 z1]},
      {-1, z2 >= Sqrt[-2 z1]},
      {0, z1 == 0 && z2 == 0}}], z1 <= 0}
   }]
eqs = ELE[[1]] /. {u[t] -> 
     u1[t] + optimalcontrol[θ1[t], θ1'[
         t]] (-2 + Cos[θ1[t]]^2)/Cos[θ1[t]]} // 
  FullSimplify
sol = NDSolve[
   {eqs, θ1[0] == -Pi/6, θ1'[0] == -1},
   θ1[t],
   {t, 0, 10}
   ];

Answer 1

Solution 1

Simplifying definition of optimalcontrol to

optimalcontrol[z1_, z2_] = 
 With[{pm1 = Piecewise[{{1, z1 > 0}}, -1]}, 
  Piecewise[{{-pm1, pm1 z2 + Sqrt[2 z1 // Abs] > 0}}, pm1]]

and make use of Simplify`PWToUnitStep and introduce a smooth approximate UnitStep:

appro = With[{k = 10^5}, ArcTan[k #]/Pi + 1/2 &];

sol = 
   NDSolveValue[
    Simplify`PWToUnitStep@
      Simplify`PWToUnitStep@{eqs, θ1[0] == -π/6, 
        θ1'[0] == -1} /. UnitStep -> appro, θ1, {t, 0, 10}]; // AbsoluteTiming
(* {0.0235535, Null} *)

Plot[sol[t], {t, 0, 10}]

Solution 2

It turns out that, after turning off DiscontinuityProcessing, discard the step control of NDSolve by setting FixedStep as the Method of NDSolve works wonders:

sol = 
   NDSolveValue[{eqs, θ1[0] == -Pi/6, θ1'[0] == -1}, θ1, {t, 0, 10}, 
     Method -> {FixedStep, Method -> Automatic, DiscontinuityProcessing -> False}, 
     StartingStepSize -> 0.001]; // AbsoluteTiming
(* {0.10727, Null} *)

ListPlot[sol, Frame -> True]

I've used the hidden syntax of ListPlot, see this post for more info.

This is surprising, I should say. It's rare to see example in which adaptive ODE solver is inferior to fixed-step ODE solver.

Solution 3

Though the default adaptive ODE solver (it's LSODA, I believe) doesn't work well on the problem, StiffnessSwitching does if we adjust MaxStepSize a bit:

sol = 
   NDSolveValue[{eqs, θ1[0] == -Pi/6, θ1'[0] == -1}, θ1, {t, 0, 10}, 
     Method -> {StiffnessSwitching, DiscontinuityProcessing -> False}, 
     MaxStepSize -> 0.1]; // AbsoluteTiming
(* {0.0140824, Null} *)

ListPlot[sol, Frame -> True, PlotRange -> All]