# I can't get to the "end of inflation" with DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697>\
**Category:** Performance\
**Tags:** question, diffeq, first-steps, physics\
**Created:** [November 24, 2020, 4:23pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697 "2020-11-24T16:23:57Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 24, 2020, 4:23pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/1 "2020-11-24T16:23:57Z")

</div>

Hello all, I’ve been struggling to get the full solution to a set of ODE’s that describe an inflationary cosmology, where the equations, the value of the constants and the initial values are given below on the code. In these equations, u[1] is proportional the volume of the universe (i.e. proportional to the cubic scale factor).  
I know from previous literature that this system undergoes very fast inflation around t=2929 but when u[1] is around 1e80 inflation stops and the system stops growing as rapidly. I’ve tried all sorts of different compilers but I usually get the error that the MaxIters are too low. When I manually change the MaxIters to be higher, the solution doesn’t get any deeper into inflation but the computation time does increase a lot. Does anybody else has any experience with solving rapidly-growing systems with Julia and can give me some advice? Do you recommend a specific algorithm (I’ve basically tried all the ones in DifferentialEquations.jl)? Any help would me most appretiatied. Below is my code so you can take a look

```julia
using Plots
using DifferentialEquations
using Roots

m=1.2e-6;
γ=0.2375;
Δ=3/(8*π*(γ^2)*0.41);
Δs=sqrt(Δ);

ϕ0=0.97;
v0=1.0;
β0=π/(2*Δs);

Ω(v,β)=(v*sin(Δs*β))/Δs

Heff(v,β,ϕ,p)=(8π*((p^2/v)+v*m^2*ϕ^2)-(6*Ω(v,β)^2/(v^2*γ^2)))

Heff0(p)=Heff(v0,β0,ϕ0,p)

p0=find_zero(Heff0, (0,1))

function geometria(du,u,p,t)
    du[1]=(3(u[1]^(1/3))^4*sin(2*Δs*(u[2])))/(2*Δs*γ)
    du[2]=2*π*γ*(u[1]^(1/3))*(m^2*(u[3])^2-((u[4])^2/u[1]^2))-((3*(u[1]^(1/3))*sin(Δs*(u[2]))^2)/(2*Δ*γ))
    du[3]=(u[4])/(u[1]^(1/3))^2
    du[4]=-(u[1]^(1/3))^4*m^2*(u[3])
end

uih = [v0;β0;ϕ0;p0]
tspan = (0.0,3000.0)
prob = ODEProblem(geometria,uih,tspan)
sol = DifferentialEquations.solve(prob,alghints=[:stiff],reltol=1e-6);

```

The following graph is as far as the solution gets and may help you visualize the problem

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/5/d50a0c320f85d47eee9e5c44cd159349ba5eb370.png)

---

<div class="post-metadata">

**Author:** ![YingboMa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yingboma/32/2181_2.png) [@YingboMa](https://discourse.julialang.org/u/YingboMa)\
**Post date:** [November 24, 2020, 10:42pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/2 "2020-11-24T22:42:33Z")

</div>

How is this model solved in the literature? 1e80 is a crazy large number. I tried to convert the initial condition to `BigFloat`s, but it doesn’t seem to help.

Setting `abstol` and `reltol` helps a bit (letting it give up the error control for `u[1]` and `u[4]`).

```julia
julia> prob = ODEProblem(geometria,big.(uih),tspan);

julia> sol = @time solve(prob, Rodas4(), reltol=[Inf, 1e-6, 1e-6, Inf], abstol=[Inf, 1e-6, 1e-6, Inf], maxiters=1000);

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase ~/.julia/dev/DiffEqBase/src/integrator_interface.jl:329
  0.560747 seconds (2.83 M allocations: 141.568 MiB, 4.10% gc time)

julia> sol.t[end]
2937.761057284057

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 25, 2020, 11:37am UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/4 "2020-11-25T11:37:33Z")

</div>

I’d double check your ODE system. At very low tolerances it seems to give you this result, which indicates that it’s truly the result of the ODE that you wrote down.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 25, 2020, 11:50am UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/6 "2020-11-25T11:50:12Z")

</div>

Any chance you could solve for log(u[1]) rather than u[1]? That would improve the scaling. Tracking dynamics in the inflation part of the solution looks really hard. I can see how the solver in an implicit method would get into real trouble there.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 25, 2020, 11:52am UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/7 "2020-11-25T11:52:19Z")

</div>

The thing is the solution acts exactly as I’d expect it to, I just can’t get any part of the solution past this point. I fully expect the solution to stop growing as rapidly when it gets to 1e80. In fact, if I choose the ImplicitEuler method with the default error tolerances I get the following graph

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/4/0484c8e2d56baf5686b4a1fed99c10ae46e2f578.png) The “only” issue here is that the system enters inflation way too early.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 25, 2020, 11:56am UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/8 "2020-11-25T11:56:13Z")

</div>

I remember trying it like that a while back (more specifically I tried solving for log(u[1])/3) and I run into the same issue. I imagine this happens because the value of the derivatives doesn’t really change so the Jacobian is just as large and the timesteps also tend to zero. Still, I will try once again for good measure.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 25, 2020, 12:00pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/9 "2020-11-25T12:00:14Z")

</div>

As far as I know, the first place where these equations were solved was [https://arxiv.org/pdf/1601.01716.pdf](https://arxiv.org/pdf/1601.01716.pdf) but there’s no mention of the numerical method. A similar set of equations in [https://arxiv.org/pdf/1401.5256.pdf](https://arxiv.org/pdf/1401.5256.pdf) was solved with the adaptive Runge-Kutta method 89 (which I think is Vern9(), for which I’ve tried both the lazy and non-lazy version and no luck, but I’ll look in the reference of the article to see if it’s something different). I tried what you did with relaxing the error and using `BigFloat` but the only thing that happened was that inflation started a little bit later (you can check that the maximum values of u[1] are of the same order) but got stuck at the same place.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 25, 2020, 12:00pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/10 "2020-11-25T12:00:48Z")

</div>

I’ve done graphs with what I have so far (that is, before the end of inflation) and they look identical to those in the literature so I think the ODE system is correct.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 25, 2020, 12:03pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/11 "2020-11-25T12:03:45Z")

</div>

For anyone wondering “how” inflation ends, or more especifically, how u[1] stops growing, u[4] becomes negative somewhere around t=2100 and then also grows very large (negatively, of course).

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 25, 2020, 12:23pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/12 "2020-11-25T12:23:04Z")

</div>

` done graphs with what I have so far (that is, before the end of inflation) and they look identical to those in the literature so I think the ODE system is correct.`

Can you convince the integrator to evaluate the Jacobian at every time step? If your solution is rapidly varying and the Jacobian is ill-conditioned, the logic in some integrators can get confused and not update the Jacobian when it should. We have run into this with DASSL years ago. We had to hack the code, but I don’t think that is necessary anymore. I think, but my memory on this may be wrong, that pvode and cvode let you do this.

@ChrisRackauckas, is there a flag @MrSantiagoPrado can set in DifferentialEquation.jl to force Jacobian updates at every time step?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 25, 2020, 12:24pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/13 "2020-11-25T12:24:13Z")

</div>

Try Rosenbrock23 or Rodas5. It’s much more natural of a method for that case.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 25, 2020, 12:33pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/15 "2020-11-25T12:33:20Z")

</div>

I’ve tried both with the same results. I will try relaxing the errors in u[1] and u[4] to see if that does it.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 25, 2020, 12:34pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/16 "2020-11-25T12:34:57Z")

</div>

I’m a bit new to using numeriacal algoriths, and DifferentialEquations.jl, but, is there any advantage in me calculating the exact Jacobian by hand and feeding it to the solver?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [November 25, 2020, 12:54pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/17 "2020-11-25T12:54:54Z")

</div>

I’ve noticed the r.h.s. of your equations mostly uses the diameter of the universe (i.e. cbrt(volume) ). If that’s the case maybe you could make that your dynamic variable and get a bit smaller numbers to work with. You can always take a power of three in your analysis

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 25, 2020, 1:30pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/18 "2020-11-25T13:30:27Z")

</div>

Most integrators default to a finite diffrerence Jacobian. That is generally ok and you’d need very little effects from changing to an analytic Jacobian. However, your Jacobian is ill-conditioned and your solution is rapidly varying during inflation. So, the finite difference error might make a difference. If your integrator uses automatic differentiation, then you already have a (computer-generated) analytic Jacobian, so you won’t gain much. For small systems, the performance difference should be pretty minimial. A hand computed Jacobian is faster than anything else if the human time to compute it is reasonable.

The good news is that your system is small and the Jacobian is easy to compute by hand. I would give it a try in this case to see what happens. My guess is that it won’t matter much, but I would be interested in how it works out. Please reply to this message after you do it.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 25, 2020, 1:43pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/19 "2020-11-25T13:43:00Z")

</div>

Do you know if you have finite time explosion (not numerically)? I mean your flow is globally defined?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [November 25, 2020, 2:55pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/21 "2020-11-25T14:55:52Z")

</div>

Ok, I think I got it: The mechanism by which V stops growing is that \beta \to 0. However, numerically there is a small residual for \beta that doesn’t vanish (it stays around 10^{-16}). So you need to fix that somehow if you want stop V from exploding.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 25, 2020, 5:41pm UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/22 "2020-11-25T17:41:47Z")

</div>

> [@ctkelley](#):
>
> Most integrators default to a finite diffrerence Jacobian.

Most don’t. Only the ones that call C/Fortran solvers do.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 26, 2020, 11:35am UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/23 "2020-11-26T11:35:22Z")

</div>

@ChrisRackauckas @ctkelley Problem solved. I just used BigFloats to define the tspan. Sorry for all the trouble.

---

<div class="post-metadata">

**Author:** ![MrSantiagoPrado](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrsantiagoprado/32/19702_2.png) [@MrSantiagoPrado](https://discourse.julialang.org/u/MrSantiagoPrado)\
**Post date:** [November 26, 2020, 11:42am UTC](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697/24 "2020-11-26T11:42:51Z")

</div>

This here was the solution. The thing I had to put in big() were the tspan though so the timesteps could be as low as needed without loss of accuracy.

[Next page](https://discourse.julialang.org/t/i-cant-get-to-the-end-of-inflation-with-differentialequations-jl/50697.md?page=2)