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.

  • $\begingroup$ You might want to look into using WhenEvent[] to handle your discontinuity. $\endgroup$ – J. M. is away Apr 12 '16 at 18:06
  • $\begingroup$ 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. $\endgroup$ – Chris K Apr 12 '16 at 19:51
  • $\begingroup$ I am having a look into the WhenEvent[] command right now. Thank you for the suggestion @J.M. $\endgroup$ – Lloyd Apr 14 '16 at 8:16
  • $\begingroup$ @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). $\endgroup$ – Lloyd Apr 14 '16 at 8:22
  • $\begingroup$ 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. $\endgroup$ – Chris K Apr 14 '16 at 12:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.