# NDSolve for an (slightly) discontinuous Piecewise input

I define a Piecewise discontinuous function V[f[t]], which depends on a function f[t] I want to get numerically using NDSolve later, as

$$V(f(t))=\left\{ \begin{array}[cc]\\ V_0 + \frac{1}{2}m^2\ f^2(t) + l\left( f^n(t) - f_0^n \right) & , f(t)>f_0 \\ V_0 + \frac{1}{2}m^2\ f^2(t) & , f(t)\leq f_0 \\ \end{array} \right.$$

where $V_0, m, f_0$ and $l$ are constant parameters. Here one can see that, taking $l\ll 1$, the step between both sides of V[f[t]] (approximately) disappears at $f_0$ and thus V[f[t]] is (almost) continuous. (This is the case for both V[f[t]] and its first derivative derV[f[t]], as I previously checked). Writing down this this as Mathematica code (I take $l\sim 10^{-18}$ as an example), and explicitely defining V[f[t]] and derV[f[t]] in order to avoid Indeterminate pieces in the definitions),

    (*  Definition of function V[f[t]] and its first derivative derV[f[t]],
and related constants   *)

l = 2.4*10^-18
n = 2/3
f0 = 10
V0 = 10
m = 6*10^-9

V[f[t_]] = Piecewise[{{V0 + (1/2)*m^2*f[t]^2 + l*((f[t])^n - (f0)^n), f[t] <= f0},
{V0 + (1/2)*m^2*f[t]^2, f[t] >= f0},
{V0 + (1/2)*m^2*f0^2, True} }]

derV[f[t_]] = Piecewise[{
{D[V0 + (1/2)*m^2*f[t]^2 + l*((f[t])^n - (f0)^n), f[t]], f[t] <= f0},
{D[V0 + (1/2)*m^2*f[t]^2, f[t]], f[t] >= f0} ,
{m^2*f0, True} }]


I want to solve two equations, EQ1 and EQ2, which depend on both V[f[t]]and derV[f[t]]. More concretely, these equations depend on two undetermined functions f[t] and g[t] which I want to NDSolve. Let me just put the corresponding code segment with the initial conditions I want to take for NDSolve later: for fand g, fini[t],\ gini[t] and for the first derivative of f, finider[t], as well as the parameter H0,

    (* Numerical parameter for initial conditions for equations EQ1 and EQ2  *)

H0 = 3.3*10^-7

(* Initial conditions for NDSolving functions f and g *)

gini[t_] = -1/(H0*t)
fini[t_] = f0*(gini[t]^l)
finider[t_] = -10*gini[t]*H0*fini[t]

(* Equations to solve, EQ1 and EQ2 *)

EQ1 = (Derivative[1][g][t]/g[t]^2) - Sqrt[(1/3)*((1/2)*(Derivative[1][f][t]/g[t])^2 + V[f[t]] )]

EQ2 = Derivative[2][f][t] + 2*(Derivative[1][g][t]/g[t])*Derivative[1][f][t] + g[t]^2*derV[f[t]]


I want to NDSolve both equations EQ1and EQ2 in order to numerically compute the functions f[t] and g[t] in an interval of times $[t_0, t_f]$, where $t_0,\ t_f\leq 0$. Taking into account all of the previous stuff I have written this part as

    (* Initial and final time *)

t0 = -10000
tf = -1000

(* NDSolve *)

y = NDSolve[{EQ1 == 0, EQ2 == 0,
g[t0] == gini[t0], f[t0] == fini[t0], Derivative[1][f][t0] == finider[t0]},
{g, f}, {t, t0, tf},
MaxSteps -> 10000000, PrecisionGoal -> 10, AccuracyGoal -> 90,
Method -> If[$VersionNumber > 8, {"DiscontinuityProcessing" -> False}, Automatic]]  where the command Method in NDSolve allows me to avoid one mistake about the Filippov continuation arising for Mathematica versions > 8.0. - as mine (10.0.). Here comes the point I got stucked at. Executing this code, Mathematica gives the mistake NDSolve::ndsz: At t == -10000., step size is effectively zero; singularity or stiff system suspected. My code works for non-Piecewise V[f[t]], so that it seems that NDSolve fails at the discontinuity point$f_0$. NDSolve stucks the discontinuity, even though$l$is defined at such a small scale as$10^{-18}$. (The surprising fact is that this mistake also happens if we define$l=0$explicitely in the code...). As a result, I think that I am not defining NDSolve in an appropiate way for handling those discontinuities. I would like to NDSolve for a determinate range of values which includes the discontinuity point f0 instead of solving separately for the two ranges of V[f[t]], i.e.,$f<f_0$and$f\geq f_0$and joining both sides of f[t] and g[t]. Is there a way or some kind of special Mathematica commands for NDSolving for such kind of step discontinuity, without failing at the discontinuity point$f_0\$? Thanks in advance.

• You might want to look into using WhenEvent[] to handle your discontinuity. – J. M. will be back soon Apr 12 '16 at 18:06
• Your code ran as-is for me, on Mathematica 10.4 on a Mac. f[t] looked constant and g[t] like a negatively sloping line. – Chris K Apr 12 '16 at 19:51
• I am having a look into the WhenEvent[] command right now. Thank you for the suggestion @J.M. – Lloyd Apr 14 '16 at 8:16
• @ChrisK Right, I obtained the same results as you. Actually that was wrong, as both EQ1and EQ2 appear like 'splitted' when copying my code into Mathematica. This happens because of how I wrote my code in my question (just for reader's convenience). Now you will notice my trouble after copying and testing the code again (or hopefully not, let us see :-) ). (I also forgot to add my value for the parameter m in the code, this is now fixed). – Lloyd Apr 14 '16 at 8:22
• I tried your code without the Piecewise by setting V[f[t_]] = V0 + (1/2)*m^2*f[t]^2 and derV[f[t_]] = D[V0 + (1/2)*m^2*f[t]^2, f[t]] and get the same NDSolve::ndsz error. So I don't think it's necessarily the Piecewise causing trouble. – Chris K Apr 14 '16 at 12:32