# Problem with NDSolve and piecewise functions: Failure computing Filippov continuation

I am trying to solve a first order differential equation using NDSolve. The differential equation is

g[x[t]] x’[t] + K1 x[t] + f[x[t]] = K2


with x’==0. Function f is a continuous linear piecewise function and g is df/dx, and it is not a continuous function. It is easy to obtain the exact solution for the different regions, and it is a periodic oscillation with a period of about 1.555 (with the data below). My problem is about the numerical solution because I can’t obtain it using NDSolve. I have tried different options, even using ListInterpolation for the points in f and g, but the results I obtain is

NDSolve::smpf: Failure to project onto the discontinuity surface when computing Filippov continuation at time 2.031744647155892.


The code is:

dataOrigIni = {a -> 100.2, b -> 0.99, x1 -> 1.99/100.2, x2 -> 1/0.99, K1 -> 13, K2 -> 2.5};
f[x_] := Piecewise[{{b x, x <= -x2}, {-(((a - b) x1 x2)/(x2 - x1)) + (b x2 - a x1)/(x2 - x1) x, -x2 <= x <= -x1}, {a x, -x1 <= x <= x1}, {((a - b) x1 x2)/(x2 - x1) + (b x2 - a x1)/(x2 - x1) x, x1 <= x <= x2}, {b x, x >= x2}}];
g[x_] := Piecewise[{{b , x <= -x2}, {(b x2 - a x1)/(x2 - x1), -x2 < x < -x1}, {a , -x1 <= x <= x1}, {(b x2 - a x1)/(x2 - x1), x1 < x < x2}, {b , x >= x2}}];
eqn0 = {g[x[t]] x'[t] + K1 x[t] + f[x[t]] == K2, x == 0} /.dataOrigIni
sol0 = NDSolve[eqn0, x[t], {t, 0, 10}]


Any suggest?

Thank you

Images from the closed solution:  I have put together Alex' solution using NDSolve and the obtained by DSolve through theoretical considerations. The comparison is below They are very similar, but jumping in x[t] is not observed and the periodicity is missing because an increasing drift can be observed. My experience with Mathematica is limited, so I don't know if this is a normal behaviour.

• How did you get exact solution with non exact parameters? Feb 28, 2021 at 11:15
• I mean I can exactly solve the differential equation in each region given by the piecewise function (I have closed expressions for x[t]). The solutions are exponential functions, first increasing and then decreasing, giving the periodic oscillation I have previously mentioned. Function f is an odd function (its shape is similar to a deformed "N" letter, more or less)
– art
Feb 28, 2021 at 11:37
• But all your parameters except K1 are Real, not exact. Therefore, your solution is not exact, but some kind of numerical solution. Even with exact parameters this kind of solution is not differentiable in some points, and therefore it can't be solved with NDSolve directly. Probably we need to add WhenEvent to this model. Feb 28, 2021 at 11:52
• Images added to the original question
– art
Feb 28, 2021 at 23:20
• The nonlinear characteristics f[x] comes from the driving point of a non-linear device (V=f[I]) inside a circuit with a linear resistor (K1), a linear capacitor and a voltage source (K2). I wanted to analyze its behaviour depending on current. I have the expressions that control the current through the circuit, the stable oscillation condition (which meet in this case), but I wanted to know if I could obtain a numerical solution in order to analyze more complex circuits.
– art
Mar 2, 2021 at 8:13

This is not solution just some remarks to the NDSolve methods. First, we can pass all interval {t,0,10} with explicit Euler, but this solution not converge to that proposed by art with number of steps increasing

dataOrigIni ={a -> 100.2, b -> 0.99, x1 -> 1.99/100.2, x2 -> 1/0.99, K1 -> 13, K2 -> 2.5} ;
f[x_] :=
Piecewise[{{b x,
x <= -x2}, {-(((a - b) x1 x2)/(x2 - x1)) + (b x2 - a x1)/(x2 -
x1) x, -x2 <= x <= -x1}, {a x, -x1 <= x <=
x1}, {((a - b) x1 x2)/(x2 - x1) + (b x2 - a x1)/(x2 - x1) x,
x1 <= x <= x2}, {b x, x >= x2}}];
g[x_] := Piecewise[{{b,
x <= -x2}, {(b x2 - a x1)/(x2 - x1), -x2 < x < -x1}, {a, -x1 <=
x <= x1}, {(b x2 - a x1)/(x2 - x1), x1 < x < x2}, {b, x >= x2}}];
eqn0 = { x'[t] + K1 x[t]/g[x[t]] + f[x[t]]/g[x[t]] == K2/g[x[t]],
x == 0} /. dataOrigIni;
sol0 = NDSolve[eqn0, x[t], {t, 0, 10},
Method -> {"DiscontinuityProcessing" -> False,
"TimeIntegration" -> "ExplicitEuler"}, MaxSteps -> 10^6,
StartingStepSize -> 4.5 10^-2]


Visualization

{Plot[x[t] /. sol0[], {t, 0, 10}, PlotRange -> All,
AxesLabel -> {"t", "x"}],
Plot[f[x[t]] /. dataOrigIni /. sol0[], {t, 0, 10},
PlotRange -> All, AxesLabel -> {"t", "f"}]} It looks very similar to the art's picture, but with StartingStepSize ->5 10^-3 it looks quiet different If we varying StartingStepSizethen we get new solution for every new step size. Therefore, there is no unique solution for this problem at t>2.031744647155892 it why we have message and NDSolve` stops.