# Tsit5() and DP5() stability

**URL:** <https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137>\
**Category:** New to Julia\
**Tags:** package, differentialequation\
**Created:** [May 15, 2026, 9:56pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137 "2026-05-15T21:56:37Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![ligeston](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ligeston/32/221913_2.png) [@ligeston](https://discourse.julialang.org/u/ligeston)\
**Post date:** [May 15, 2026, 9:56pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/1 "2026-05-15T21:56:37Z")

</div>

Hi guys,

I am wrapping up a script to solve the [TOV equation](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation), similar to the [TOVsolver\_Julia](https://github.com/jskMNMGCH/TOVsolver_Julia/tree/main), and I used ‘CommonSolve’ and ‘DifferentialEquations’ packages. In the solver, unlike the case in [TOVsolver\_Julia](https://github.com/jskMNMGCH/TOVsolver_Julia/tree/main), I used callbacks to terminate the integration when P\<0

```julia-auto
    condition(u, t, integrator) = u[1] < 0
    affect!(integrator) = terminate!(integrator)
    cb = DiscreteCallback(condition, affect!)

```

and applied Tsit5 and DP5 (the same as [TOVsolver\_Julia](https://github.com/jskMNMGCH/TOVsolver_Julia/tree/main)). However, I found some unsmooth results that did not occur in my hand-crafted RK4 method in the mass-radius curves for different initial conditions for the density at r=0,

**Using solve in CommonSolve**

 ![Screenshot 2026-05-15 at 10.54.24](https://global.discourse-cdn.com/julialang/original/3X/2/9/2915910f3b18137e90994bb38c76249dca050588.png)

In addition, I found out that in some cases, the integration becomes unstable, and a warning pops up for both Tsit5 and DP5 algorithms.

```julia-auto
> Warning: Verbosity toggle: dt_epsilon 
> │ At t= 9656.595735553243, dt was forced below floating point epsilon 1.8189894035458565e-12. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64.
> └ @ SciMLBase ~/.julia/packages/SciMLBase/TKP6E/src/integrator_interface.jl:736

```

Does anyone know if there is a better approach for this type of differential equation?

---

<div class="post-metadata">

**Author:** ![langestefan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/langestefan/32/207923_2.png) [@langestefan](https://discourse.julialang.org/u/langestefan)\
**Post date:** [May 15, 2026, 10:31pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/2 "2026-05-15T22:31:25Z")

</div>

Can you try `ContinuousCallback` instead? It can hit within the integration step, so it should be more stable

---

<div class="post-metadata">

**Author:** ![ligeston](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ligeston/32/221913_2.png) [@ligeston](https://discourse.julialang.org/u/ligeston)\
**Post date:** [May 15, 2026, 10:44pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/3 "2026-05-15T22:44:19Z")

</div>

> [@langestefan](#):
>
> ContinuousCallback

Thanks, just tried `ContinuousCallback` and also switched the callback off for another test. The same warning occurs.

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [May 16, 2026, 11:02am UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/4 "2026-05-16T11:02:40Z")

</div>

If RK4 is OK, why don’t you stick to it?

---

<div class="post-metadata">

**Author:** ![ligeston](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ligeston/32/221913_2.png) [@ligeston](https://discourse.julialang.org/u/ligeston)\
**Post date:** [May 16, 2026, 9:32pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/5 "2026-05-16T21:32:11Z")

</div>

There could be some discontinuities in the ODE, so I am not sure RK4 is the best approach.

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [May 17, 2026, 2:25am UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/6 "2026-05-17T02:25:28Z")

</div>

If there are discontinuities, lower order methods are generally expected to be better than higher order ones.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 17, 2026, 4:11am UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/7 "2026-05-17T04:11:33Z")

</div>

You might want to try out `BS3` to see how that does.

---

<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:** [May 17, 2026, 10:13am UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/8 "2026-05-17T10:13:17Z")

</div>

> [@ligeston](#):
>
> I found out that in some cases, the integration becomes unstable, and a warning pops up for both Tsit5 and DP5 algorithms.
> 
> ```julia-auto
> 
> ```

To be clear, this isn’t saying that it’s unstable, it means that it’s unable to achieve the accuracy to the requested tolerance.

This to me all looks and sounds like you did an invalid mutation. Is your `f` actually defining an ODE, i.e. do you get the same result from evaluating `f(du,u,p,t)` each time, or are you doing something like mutating `u` or trying to store some cache of past values? Without any code it’s hard to debug, but it seems like the one difference is the moment it starts to have a single step rejection it’s not able to recover, and that heavily points towards this.

This video could be helpful:

[![](https://global.discourse-cdn.com/julialang/original/3X/5/d/5d8be9be2d861eecc910a9bdda15b5b546dd839f.jpeg "How to Debug Julia Simulation Codes (ODEs, Optimization, etc.!) | Chris Rackauckas | JuliaCon 2022") ](https://www.youtube.com/watch?v=g-iOOhh2U6o)

---

<div class="post-metadata">

**Author:** ![ligeston](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ligeston/32/221913_2.png) [@ligeston](https://discourse.julialang.org/u/ligeston)\
**Post date:** [May 18, 2026, 10:55pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/9 "2026-05-18T22:55:50Z")

</div>

Thanks. I am going to have a look for the video. Below I attached my script. One of the purposes of this test script is that there will be a loop to change the initial condition. It seems to me that for a single run, the results look ok, but the loop for changing the parameter is also terminated.

```julia-auto
import OrdinaryDiffEq as DE
using OrdinaryDiffEqLowOrderRK: OrdinaryDiffEqLowOrderRK
using SciMLBase: ODEProblem, DiscreteCallback, terminate!
using Plots

    const planck_constant = 6.62607015e-34
    const ħ = planck_constant / (2*pi)
    const gravitational_constant = 6.67408e-11
    const speed_of_light_vacuum = 2.99792458e8
    const neutron_mass = 1.674927471e-27

EoS_Param = (
    K_crust = (3*π^2)^(2/3) * ħ^2 / (5*neutron_mass^(8/3)), # Polytropic constant for the crust
    γ_crust = 5.0/3.0, # Polytropic index for the crust
    γ_core = 3.0, # Polytropic index for the core
    ρ_b = 3e14*1e3, # Transition density between crust and core in kg/m^3
)

# Simulation Input Parameters
TOV_Sol_Input = (
    ρ0 = 1e16*1e3, # Initial central density in kg/m^3 (1e15 g/cm^3)
    dr = 0.01*1e3, # Radial step in meters
    Dr = 0.01*1e3, # Radial interval for recording values in meters
    r_end = 20e3, # Maximum radius to solve up to in meters
);

function EoS_two_component_polytrope(
    ParamIn::NamedTuple;
    report_transition_pressure::Bool = false,
)
    K_core = ParamIn.K_crust * ParamIn.ρ_b ^ (ParamIn.γ_crust - ParamIn.γ_core)

    function EoS_P_from_rho(ρ)
        P = if ρ < 0
            0.0
        elseif ρ <= ParamIn.ρ_b
            ParamIn.K_crust * ρ .^ ParamIn.γ_crust
        else
            K_core * ρ .^ ParamIn.γ_core
        end
        return P
    end

    # Pressure at the crust-core transition for continuity check
    P_b = EoS_P_from_rho(ParamIn.ρ_b)

    if report_transition_pressure
        println("Pressure at crust-core transition (Pb): ", P_b)
    end

    function EoS_rho_from_P(P)
        ρ = if P < 0
            0.0
        elseif P <= P_b
            (P / ParamIn.K_crust) .^ (1/ParamIn.γ_crust)
        else
            (P / K_core) .^ (1/ParamIn.γ_core)
        end
        return ρ
    end

    return EoS_P_from_rho, EoS_rho_from_P
end

function tov_eq!(EoS_rho_from_P::Function)

    function tov_inner!(du, u, paras, r)
        P = u[1]
        m = u[2]
        ρ = EoS_rho_from_P(P)
        du[1] = if r == 0.0
            0.0
        else
            -gravitational_constant / r^2 *
            (ρ + P/speed_of_light_vacuum^2) *
            (m + 4*pi*r^3*P/speed_of_light_vacuum^2) /
            (1 - 2*gravitational_constant*m/(r*speed_of_light_vacuum^2))
        end
        du[2] = 4*pi*r^2*ρ
    end
    return tov_inner!
end

EoS, EoS_inv= EoS_two_component_polytrope(EoS_Param);
u0 = [EoS(TOV_Sol_Input.ρ0); 0.0]; # Initial conditions: central pressure and enclosed mass
tov! = tov_eq!(EoS_inv);

condition(u, t, integrator) = u[1] < 0
affect!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition, affect!)

problem = ODEProblem(tov!, u0, (0.0, TOV_Sol_Input.r_end); callback = cb)
sol = DE.solve(problem, OrdinaryDiffEqLowOrderRK.DP5(); reltol = 1e-12)
#sol = DE.solve(problem, DE.Tsit5(); reltol = 1e-12)

ur = Array(sol)
r = sol.t
Pr = ur[1, :]
mr = ur[2, :]
ρr = EoS_inv.(Pr)
plot(r,ρr)

```

---

<div class="post-metadata">

**Author:** ![ligeston](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ligeston/32/221913_2.png) [@ligeston](https://discourse.julialang.org/u/ligeston)\
**Post date:** [May 18, 2026, 11:00pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/10 "2026-05-18T23:00:39Z")

</div>

Sounds a good idea, and I just tried. The same warning popped up. I had a look for the results by in the trimmed version (without a loop to vary parameters), a single run even with a warning looks actually fine. It may also terminate my loop I guess, so I could not get the results I want.

---

<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:** [May 19, 2026, 9:46am UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/11 "2026-05-19T09:46:56Z")

</div>

`sol = solve(prob, DP5(); reltol=1e-8, abstol=[1.0, 1e15])` seems to work fine: per component tolerances help a lot here since the values are wildly different.

```julia-auto
ρ0=1.0e18 retcode=Terminated R=13168.7m M=4.35e30 kg

```

The default absolute tolerances is not great for the value as unscaled. The other thing you could (and probably should) do is rescale the model towards 1.

Or you can put the continuous callback just before the surface.

---

<div class="post-metadata">

**Author:** ![ligeston](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ligeston/32/221913_2.png) [@ligeston](https://discourse.julialang.org/u/ligeston)\
**Post date:** [May 19, 2026, 10:19am UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/12 "2026-05-19T10:19:58Z")

</div>

Goodness. Thank you so much. And yes, I will seriously think about how to deal with those very large and very small constants in the equation.

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [May 19, 2026, 2:15pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/13 "2026-05-19T14:15:47Z")

</div>

> [@ligeston](#):
>
> And yes, I will seriously think about how to deal with those very large and very small constants in the equation.

You’ll probably want to convert your equations to [Natural Units](https://en.wikipedia.org/wiki/Natural_units)  
i.e. `\hbar = 1,...`

---

<div class="post-metadata">

**Author:** ![sablonl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sablonl/32/217805_2.png) [@sablonl](https://discourse.julialang.org/u/sablonl)\
**Post date:** [May 19, 2026, 2:32pm UTC](https://discourse.julialang.org/t/tsit5-and-dp5-stability/137137/14 "2026-05-19T14:32:02Z")

</div>

ModelingToolkit.jl is also quite handy for this kind of problem when it comes to rescaling. And since you have one system that you will solve many time for different parameters and/or initial conditions, you can clearly afford the mtk-compilation once and then use `EnsembleProblem` [Parallel Ensemble Simulations Interface · SciMLBase.jl](https://docs.sciml.ai/SciMLBase/dev/interfaces/Ensembles/)