I am finding some troubles when using NDSolve with this piece of code:

    (*  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 

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

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

    (* 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_] = -Abs[1/(H0*t)]   
    fini[t_] = Abs[f0*(gini[t]^l)]
    finider[t_] = -Abs[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]]

    (* Initial and final time *)

    t0 = 1 
    tf = 1*10^9 

    (* NDSolving *)

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

As I think, my problems are related with the time-dependent function V[f[t]], defined as

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

where $l\ll 1$, so that both pieces of V[f[t]] (approximately) coincide at $f_0$. (This is the case for both V[f[t]] and its first derivative, as I checked)

Let me recall that at the beginning I needed to introduce absolute values in the functions fini, finider and gini in order to avoid the appearance of complex numbers in NDSolve and thus the corresponding mistakes. After that, I got the mistake related to Filippov continuation, which I solved unplugging the "DiscontinuityProcessing" method as shown in the last command of NDSolve (I am working with Mathematica 10).

However, the code is still failing at the discontinuity point, i.e., NDSolve detects the singularity at $f_0$. I do not know if Mathematica is not able of handling V[f[t]] as I wrote down (I am trying with $l\sim 10^{-20}$, even I tried with much smaller values, but perhaps Mathematica is still handling both pieces of V[f[t]] as discontinuous) or if I have to add another instruction of Mathematica which I do not know.

What could I do? Thank you very much in advance.


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.