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

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

• There are examples in the docs. If those approaches do not work for you, then the details of your problem are needed to figure out a better approach, if in fact it is possible. Dec 4, 2022 at 0:46
• Sorry, SE was struggling with the latex in the post and only the first paragraph posted. The full post is now up. Dec 4, 2022 at 0:52
• @NolanKing Please, pay attention, that $\theta 1=\pm \pi/2, \theta 1'=0$ are also equilibrium states. Dec 4, 2022 at 6:01
• Please show us definition of eqs written with Mathematica code rather than $\LaTeX$. If you're having difficulty in making it look nice, press Ctrl+Shift+I to convert it to InputForm before copying. For more info, read this: mathematica.meta.stackexchange.com/a/1585/1871 Dec 4, 2022 at 8:21
• The code is attached now. Dec 4, 2022 at 13:52

## 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 SimplifyPWToUnitStep and introduce a smooth approximate UnitStep:

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

sol =
NDSolveValue[
SimplifyPWToUnitStep@
SimplifyPWToUnitStep@{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]
`