I have these differential equations and I wanted to solve its solutions then plug it into an expression for temperature (which I have listed in the plot below), but the expression for temperature is a piece-wise function one for $t<t_b$ and the other is $t>t_b$.

Mp = 2.4353*10^18;
m = 1.8*10^13;
Rm = (1.8*10^13)/(2.4353*10^18);
v = 13.2 Mp;
Rv = 13.2;
\[Lambda] = 10^-13;
Tb = 10 (m/Mp)^4;
tf = 10^7; 
sol = NDSolve[{dy'[t] + 3 H[t] (1 + 10) dy[t] + \[Lambda] (y[t]^2 - Rv^2) y[t] == 0, H[t] == Sqrt[(0.5 dy[t]^2 + 0.25 \[Lambda] (y[t]^2 - Rv^2)^2 + \[Rho][t])/3], \[Rho]'[t] + 4 H[t] \[Rho][t] == 3 H[t] 10 dy[t]^2, y'[t] == dy[t], y[0] == -30, dy[0] == 9*10^-6, \[Rho][0] == Rm^4}, {y, dy, H, \[Rho]}, {t, 0, tfin},MaxStepFraction -> 0.001, MaxSteps -> 5*10^7]
ysol = y[10^7] /. sol;
ypsol = dy[10^7] /. sol;
Hsol = H[10^7] /. sol
\[Rho][10^7] /. sol;
tb = ((3 Hsol (1 + 10))/(2 \[Lambda] Rv^2)) (Log[(40 Pi^2 Hsol^2 (1 + 10)^2)/(\[Lambda]^2 Rv^6 10)])
Plot[Piecewise[{{Tb Exp[Hsol (x - tb)],x < tb}, {Tb Exp[(\[Lambda] Rv^2 (x-tb))/(6 Hsol (1 + 10))], x > tb}}], {x, 0, 10^7}, PlotRange -> All]

It gave me an error that tb is not a machine sized number, etc. Can anyone help me with this situation? Also, I'm fairly new with mathematica so I've got few more problems, does my code look correct? For example, I substituted "sol" to a function of x, but every quantity needs to be evaluated at some t, so I'm not sure. Just a note, the number 10 you see is something that I can change later so I just left it like that, example (1+10), I did not add it so that it is easier to change later.

Basically, I need to plot $T$ vs. $t$

ADDITIONAL QUESTION: Suppose I have these codes, can anyone help me determine the proper $tf$ and $y[0]$ to get a plot like below, note that in my code $T = \frac{3 \rho[t]}{10 \pi^2}$ where $T$ is the temperature. Also, dy[0]=b shouldn't be too sensitive, $10^{-6}$ would be ok, it is just $y[0]$ and $tf$ that I somehow have some difficulties, although $0 \lt y[0] \lt 50$ (maybe it can go beyond 50 but 500 is too much I think). For some combination of $y[0]$ and $tf$, I'm getting an error saying ""For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions". I'm not sure how to determine the correct combination to get the plot, it's not necessary to get the exact same plot, but $T$ should somehow cut the line $H$ and go below it, like a slanted V.


Mp = 2.4353*10^18;
m = 1.8*10^13;
Rm = (1.8*10^13 )/(2.4353*10^18);
\[Sigma] = 2.24*10^19;
\[Sigma]m = (2.24*10^19)/(2.4353*10^18);
\[Lambda] = ((1.8*10^13)/(2.4353*10^18))^4;
tf = 10^9;
sol[a_, b_] := NDSolve[{dy'[t] + 303 H[t] dy[t] + 4 \[Lambda] (y[t]^2 - \[Sigma]m^2) y[t] == 0, H[t] == Sqrt[(0.5 dy[t]^2 + \[Lambda] (y[t]^2 - \[Sigma]m^2)^2 + \[Rho][t])/3], \[Rho]'[t] + 4 H[t] \[Rho][t] == 300 H[t] dy[t]^2, y'[t] == dy[t], y[0] == -a, dy[0] == b, \[Rho][0] == Rm^4}, {y, dy, H, \[Rho]}, {t, 0, tf}]
sol[1, 1/100000]
Manipulate[Plot[{Evaluate[H[t] /. First@sol[a, b]], Evaluate[(3 \[Rho][t]/10 Pi^2)^(1/4) /. First@sol[a, b]]}, {t, 0, tf}, PlotRange -> Automatic], {{a, 1}, 0, 50, Appearance -> "Labeled"}, {{b, 1/100000}, 0, 0.00009, Appearance -> "Labeled"}]

I have resolved the previous problem but suppose I wanted to plot $H$ and $T = (3\rho(t)/10 \pi^2)^{1/4}$ vs. $N$. I have written the code below but when I plot them it just shows the "$N$" axis, no axis for $H$ and $T$. What seems to be wrong in the code? I tried figuring out if I have some syntax error etc. but so far, I don't see anything wrong.

Mp = 2.4353 10^18; 
m = 1.8 10^15; 
Rm = (1.8 10^15 )/(2.4353 10^13); 
\[Sigma] = 2.24 10^19; 
R\[Sigma] = (2.24 10^19)/(2.4353 10^18); 
\[Lambda] = ((1.8 10^15)/(2.4353 10^18))^4; 
tf = 10^8; 
H[t_] = Sqrt[1/3 (0.5 y'[t]^2 + \[Lambda] (y[t]^2 - R\[Sigma]^2)^2 + \[Rho][t])]; 

sol = ParametricNDSolveValue[Rationalize[{y''[t] + 3 (1 + 0.5) H[t] y'[t] + 4 \[Lambda] (y[t]^2 - R\[Sigma]^2) y[t] == 0, \[Rho]'[t] + 4 H[t] \[Rho][t] == 1.5 H[t] y'[t]^2, D[Nx[t], t] == H[t], y[0] == -a, y'[0] == b, \[Rho][0] == Rm^4, Nx[0] == 0}, 0], {y, \[Rho], Nx}, {t, 0, tf}, {a, b}, AccuracyGoal -> 16, MaxSteps -> Infinity] 

Manipulate[ParametricPlot[Evaluate[{{Nx[t], ((3 \[Rho][t])/(10 \[Pi]^2))^(1/4)}, {Nx[t], H[t]}} /. Thread[{y, \[Rho], Nx} -> sol[a, b]]], {t, 0, tf}, PlotRange -> Automatic], {{a, 1}, 0, 2, Appearance -> "Labeled"}, {{b, 1/100000}, 0, 0.00009, Appearance -> "Labeled"}]
  • 1
    $\begingroup$ The value for tfin seems to be missing and Hsol "seems funny" $\endgroup$
    – Bill
    Commented Apr 25, 2017 at 18:51
  • $\begingroup$ @Bill I have edited my post, another problem is that I don't know how to determine the proper time $tf$ to calculate the differential equation, I'm just guessing, maybe $10^7$, and the proper initial conditions for the $dy[0]=10^{-6}$ is ok but for $y[0]$ I'm not sure. I should get a plot of $T$ vs $t$ such that the plot looks like a slanted "V", the function for t<tb should die but will be connected by the function for t>tb which should be increasing, but as I know the duration for the function t>tb takes longer to increase therefore the whole plot should be slanted. $\endgroup$
    – mathemania
    Commented Apr 25, 2017 at 19:11
  • 1
    $\begingroup$ Still, tfin is missing, you only have a tf defined in your code, if you modify tfin in NDSolve to tf, and all the /.sol to /.sol[[1]] (please think about why {2} > 1 doesn't "work") then everything works well, the resulting graph isn't interesting though, clearly something is wrong with tb or tf. (tf is much smaller than tb! ) $\endgroup$
    – xzczd
    Commented Apr 26, 2017 at 2:30
  • $\begingroup$ @mathemania To have an output not just an empty plot, do this with tb=Flatten[....][[1]]. $\endgroup$
    – zhk
    Commented Apr 26, 2017 at 2:35
  • $\begingroup$ @xzczd Thanks for that info. Is there a way to better determine the proper tf, I mean in general? Or is it just a guessing game? $\endgroup$
    – mathemania
    Commented Apr 26, 2017 at 5:17

1 Answer 1


Not a answer, but too long for a comment.

Well, as said above, I'm afraid it's hard to determine a and b just from a coding perspective. Put even more bluntly, currently we don't even know whether your modeling is correct or not. Anyway, your coding can be improved, here are some suggestions:

  1. NDSolve is able to handle high-order derivative without any difficulty, there's no need to rewrite y''[t] to dy'[t], and as far as I can tell, rewriting y'[t] as dy[t] won't make NDSolve work any better.

  2. On the other hand, eliminating H[t] from the equation set will improve the performance of NDSolve, because there's no derivative for H[t] in the equation set, NDSolve will treat the equations as a DAE system. (Indeed, in some cases NDSolve will transform the DAE system to an ODE system, but not always, and the warning about IDA you encountered is exactly an evidence for the triggering of DAE solver, just search IDA in the document. ) The DAE solver of NDSolve is weaker compared to the ODE solver, so it's good to take the expression of H[t] out of NDSolve.

  3. You can make use of ParametricNDSolveValue since a and b is under-determined.

  4. Since the parameters are extreme, we need to raise the WorkingPrecision of NDSolve. Also, we need to use Rationalize to raise the precision of the parameters because you've used MachinePrecision numbers. (Actually in this case we just need to raise AccuracyGoal to e.g. AccuracyGoal -> 16, see this post for more information. Given AccuracyGoal is a subtle option, I'll stick to WorkingPrecision in the rest part of this post. )

  5. The graph you posted is a logarithmic graph, so we can use LogPlot.

  6. The position of Evaluate is wrong, it should be outside of the first argument of Plot, or it won't "break through" the HoldAll attribute of Plot.

  7. With all due respect, you should put a bit more effort into learning the core language of Mathematica.

The following is an improved version of the code:

Mp = 2.4353 10^18;
m = 1.8 10^13;
Rm = (1.8 10^13)/(2.4353 10^18);
σ = 2.24 10^19;
σm = (2.24 10^19)/(2.4353 10^18);
λ = ((1.8 10^13)/(2.4353 10^18))^4;
tf = 10^9;
H[t_] = Sqrt[1/3 (0.5 y'[t]^2 + λ (y[t]^2 - σm^2)^2 + ρ[t])]
sol = ParametricNDSolveValue[
  Rationalize[{y''[t] + 303 H[t] y'[t] + 4 λ (y[t]^2 - σm^2) y[t] == 0, 
    ρ'[t] + 4 H[t] ρ[t] == 300 H[t] y'[t]^2, y[0] == -a, 
    y'[0] == b, ρ[0] == Rm^4}, 0], {y, ρ}, {t, 0, tf}, {a, b}, 
  WorkingPrecision -> 32]

LogPlot[{H[t], (3/10 ρ[t] π^2)^(1/4)} /. 
   Thread[{y, ρ} -> sol[200, 10^-1]] // Evaluate, {t, 0, tf/100}, PlotRange -> All]

Mathematica graphics

As we can see, under the parameters chosen above, though not crossing curve of $H$, the curve of $T$ bends. Though the goal isn't achieved, it's a small step forward, I think.

  • $\begingroup$ If the Evaluate is outside then what is it evaluating? Why is there a // beside Evaluate? Also, what does Thread do here? I've read the documentation but I don't really understand what it does. Where did you specify the initial conditions, is it sol[200,10^-1]? $\endgroup$
    – mathemania
    Commented Apr 28, 2017 at 5:53
  • $\begingroup$ @mathemania 1. As said above, Evaluate here evaluates the whole first argument i.e. {H[t], (3/10 ρ[t] π^2)^(1/4)} /. Thread[{y, ρ} -> sol[200, 10^-1]]. 2. For example, Thread[{a, b} -> {c, d}] evaluates to {a -> b, c -> d}, that's exactly what I've done here. Actually if you feel confused about what Thread has done here, you can simply copy the Thread[{y, ρ} -> sol[200, 10^-1]] out and execute it and observe. 3. Yes, don't forget the sol here is a ParametricFunction. 4. I believe you'll find this post interesting. $\endgroup$
    – xzczd
    Commented Apr 28, 2017 at 6:01
  • $\begingroup$ @mathemania 5. As to the meaning of //, just move your cursor to the right of it (or to the middle of it, or simply select it by double click) and press F1 and read its document. Also, check this post. $\endgroup$
    – xzczd
    Commented Apr 28, 2017 at 6:17
  • $\begingroup$ 1. So // is somehow the same as [ ], right? {H[t], (3/10 ρ[t] π^2)^(1/4)} /. Thread[{y, ρ} -> sol[200, 10^-1]] // Evaluate is equivalent to Evaluate[{H[t], (3/10 ρ[t] π^2)^(1/4)} /. Thread[{y, ρ} -> sol[200, 10^-1]]] 2. So Thread is like Evaluate but in a more specific way, such that if there are many things that needs to be evaluated, like a matching type of thing. $\endgroup$
    – mathemania
    Commented Apr 28, 2017 at 6:49
  • $\begingroup$ @mathemania 1. Yes. 2. No, Thread is a function for list manipulation, as shown in the document, it can change f[{a, b, c}, {x, y, z}] to {f[a, x], f[b, y], f[c, z]}. (Don't forget Rule(->) is a function i.e. it just plays the role of f here. BTW, you should also read those examples when reading the document. ) Evaluate is for adjusting evaluation order. Evaluate and Thread are totally different. $\endgroup$
    – xzczd
    Commented Apr 28, 2017 at 7:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.