Numerical artifact in solution of system of ODEs

Specific Domains Numerics
1

I am currently solving the following system of ODEs:

\begin{cases} \ddot\phi=-3H\dot{\phi}-\phi\\ \dot{a}=Ha\\ \ddot\sigma=-H\dot{\sigma}-\frac{m_{Eff}^2}{a^2}\sigma\\ \ddot{\rho}_0=\frac{4\rho_2}{a^2}-H\dot{\rho}_{0}-\frac{4\omega^2\rho_{0}}{a^2} \end{cases}

where the last equation is solved for different k\in[1, 10^4] and

\begin{cases} H=\sqrt{\frac{\dot{\phi }^2+\phi ^2}{6}}\\ R=\dot{\phi }^2-2\phi^2\\ m_{Eff}^2=-a^2R\left(\xi-\frac{1}{6} \right)\\ \omega ^2=k^2+m_{Eff}^2\\ \rho _2=\frac{1}{2\rho _{0}}\left(0.25 + \omega^2\rho _0^2+\frac{\dot{\rho }_{0}^2}{4}\right).\end{cases}

Here’s an MWE of it (with \xi = \frac{1}{5}):

using DifferentialEquations
using Plots

function solve_equations()
    kavlength = 599
    k = Vector(logrange(1, 10^4, kavlength + 1))
    kav = (k[1:end - 1] + k[2:end]) / 2
    iconds = initialConditions(kav, kavlength)
    tspan = (0.0, 1.5)
    teval = range(tspan[1], tspan[2], 501)
    ddrho0 = zeros(kavlength)
    omega2 = zeros(kavlength)

    prob = ODEProblem(diffEquations!, iconds, tspan, (kav, kavlength, ddrho0, omega2), saveat=teval)
    sol = solve(prob, Vern9(), abstol=1e-15, reltol=1e-15)
    # savefig(heatmap(sol.t, kav, log10.(abs.(sol[6:5 + kavlength, :])), yscale=:log10, dpi=600, title="\\delta\\rho_0", xlabel="t", ylabel="k"), "heatmap.svg")
end

function initialConditions(kav, kavlength)
    rho0 = ones(kavlength) ./ (2 .* sqrt.(16 .+ kav.^2))
    drho0 = zeros(size(rho0))
    arr = [15, -0.8, 1, 0.001, 0]
    initial = vcat(arr, rho0, drho0)
    return initial
end

function ddrho0!(a, H, omega2, rho0, drho0, ddrho0)
    for i in eachindex(omega2)
        rho2 = (0.25 + omega2[i] * rho0[i]^2 + drho0[i]^2 / 4) / (2 * rho0[i])
        ddrho0[i] = 4 * rho2 / a^2 - H * drho0[i] - 4 * omega2[i] * rho0[i] / a^2
    end
    nothing
end

function diffEquations!(du, u, args, t)
    kav, kavlength, ddrho0, omega2 = args

    infl, dinfl, a, sigma, dsigma = @view u[1:5]
    rho0 = @view u[6:6 + kavlength - 1]
    drho0 = @view u[6 + kavlength:6 + 2 * kavlength - 1]

    H = sqrt((dinfl^2 + infl^2) / 6)
    da = a * H
    R = dinfl^2 - 2 * infl^2
    
    m2Eff = -a^2 * R / 30
    @. omega2 = kav^2 + m2Eff
    ddrho0!(a, H, omega2, rho0, drho0, ddrho0)
    ddinfl = -3 * H * dinfl - infl
    ddsigma = - H * dsigma - m2Eff * sigma / a^2

    du[1] = dinfl
    du[2] = ddinfl
    du[3] = da
    du[4] = dsigma
    du[5] = ddsigma
    du[6:5 + kavlength] .= drho0
    du[6 + kavlength:5 + 2 * kavlength] .= ddrho0
    nothing
end

solve_equations()

The solver actually errors at t\approx1,35 with the error message

Warning: At t=1.3482231005865255, dt was forced below floating point epsilon 2.220446049250313e-16, and step error estimate = 14.835297241231595. 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).

From the solution one can create a heatmap of log_{10}(\delta\rho_0) which looks like this:

Heatmap of log10(deltarho0)
Heatmap of log10(deltarho0)3600Ă—2400 192 KB

There a numerical artifact can be seen starting from around t=0,6. I believe another one is appearing from the point the solver fails as well. The solver already had issues getting over the first line, slowing down considerably. Changing the density of k values doesn’t change the line and changing k_{min} and k_{max} doesn’t change its position either, just the part that can be seen. The solver tolerances are already at 1e-15, changing them one way or another doesn’t make a difference either. Digging into the values of the solution the only notable thing is a seemingly large gap in the derivative of \delta\rho_0 at where the line is. Here’s an example of the derivative of \delta\rho_0 evaluated for k\approx1:

t derivative of deltarho0
0.564 -0.140
0.567 -0.137
0.57 -0.134
0.573 -0.130
0.576 -0.127
0.579 -0.123
0.582 -0.116
0.585 0.114
0.588 0.116
0.591 0.116
0.594 0.115
0.597 0.114
0.6 0.112
0.603 0.110

How would a jump (and line) like this appear?

2

What are k and \xi in relation to the differential equation? They do not appear in the ODE.

3

They are effectively constants. \xi is a parameter of the system, in the MWE I’ve set it to 0,2. k is also a constant but the fourth DE is evaluated with different values of k within the range k\in [1, 10^4].

4

In the code, \xi = 1/5 = 0.2 (using the English decimal point…), thus, \xi - 1/6 = 1/30.

5

Are you sure it’s not a singularity or instability in your equation? e.g. \rho_2 blows up whenever \rho_0 passes through zero, which could account for the apparent discontinuity in \dot{\rho}_0.

2 Likes
6

I’m quite confident in that (albeit less so after seeing it staying at the same place with different k_{min} and k_{max}). \rho_0 shouldn’t, and indeed doesn’t, go negative, at least on the points that were saved. It is possible, however, that I’ve misstyped something in the equations that’s creating this effect.

7

I’m not 100% sure if this is a viable idea because I don’t know much about solving systems of ODEs, but you could try using my new package PrecisionCarriers.jl to see if your functions maybe become numerically unstable around that line. There are a few subtractions in your functions that could lead to catastrophic cancellations for example.

I already tried using your code myself and solving the entire thing using the package, but it seems that DifferentialEquations.jl casts the values to Float64 somewhere (in the Vern9 cache?), which unfortunately gets rid of the carried precision. I’m not sure if this is something I could fix for PrecisionCarriers.

1 Like
8

Could you break the code down to a smaller example? If I understand you correctly, k is a fixed parameter in the set of ODEs you solve, but in the code it looks like you solve the ODEs for many different values of k simultaneously in one big equation.

If these equations all decouple, I’d wager it’s faster and more stable to solve them individually with one problem per value of k (perhaps with an EnsembleProblem). Otherwise, you are constraining the solver for all values of k to the equations with the “most unfavorable” k-values.

In any case, a simpler MWE would make it easier to dig further into the issue. Unless you are actually solving a discretized PDE, in which case other things might be going on.

9

What’s going on with a? It looks unstable, because \dot{a} = Ha and H is a positive square root.

All sorts of weird things happen, for example m^2_{Eff}=-a^2... should get big (not sure if positive or negative), and so \omega^2 should get big, and so \rho_2 should get big even if \rho_0 doesn’t pass through zero?

Several differential right-hand-sides multiply and then divide by a^2 in code, which doesn’t seem like good policy as a blows up. All this could be wrong depending on the constants, but it’s hard to see what’s going on, hence questions.

The dynamics look suspicious to me. I wonder if there is a simplification of the equations that eliminates need to multiply and divide by an unstable quantity.

As with @Sevi, I recommend ditching the sweep of k, since you get the artifact for all small values seen in the plot (and probably larger ones as well). Just choose one k, and plot all the states over time to see what they do. And straighten out the situation with a.

4 Likes
10

Thank you for your replies. Unfortunately, as so it often happens, the problem was between the keyboard and the chair. While writing the MWE for a single k I went over the equations once more and noticed that I had forgotten a single a^2 term from the equation for \rho_2 (which should’ve been \rho_2 = \frac{1}{2\rho_0}\left(0.25+\omega^2\rho_0^2+\frac{a^2\dot\rho_0^2}{4}\right)). With this correction the “artifact” disappeared and @stevengj was indeed correct, it was an instability of the system.

Notwithstanding, I edited the calculations slightly based on @apo383 's suggestions, mainly changing m_{Eff}^2 for \frac{m_{Eff}^2}{a^2} = -R\left(\xi - \frac{1}{6}\right) and \omega^2 for \frac{\omega^2}{a^2} = \frac{k^2}{a^2} + \frac{m_{Eff}^2}{a^2} and merging the equations for \ddot\rho_0 and \rho_2 to \ddot\rho_0 = \frac{1}{2\rho_0}\left(\dot\rho_0^2 + \frac{1}{a^2}\right) - H\dot\rho_0 -2\rho_0\frac{\omega^2}{a^2}. Hopefully they’ll help with the stability with larges values of t.

rho0 without a huge rift in the middle
rho0 without a huge rift in the middle3600Ă—2400 183 KB

5 Likes
11

It looks like a good idea to get rid of \rho_2. An additional suggestion is to replace dividing by a^2 with multiplying by b^2, where \dot{b}=-Hb. The problem with a is that it’s unstable and heading to an overflow error. Whereas b is decaying exponentially, and will gracefully underflow to zero.

Perhaps \rho_0 is unstable so you’re headed for trouble no matter what, but still seems best to isolate a single instability, otherwise it’s a race to see which one overflows first.

In general I also find it helpful to set expectations with a toy problem. The fact that you complained about the line artifact suggests you didn’t have an expectation on the overall behavior. A toy problem is one with trivial parameter values, where you set some to zero or infinity or something easy to reason with. Right now it’s hard to reason about your system because it’s high order and coupled. Perhaps \phi and \sigma are supposed to be stable, so can you set them to zero or constant? If you thus decouple \phi_0, it may be evident whether it should be unstable.

1 Like
12

Did you make u0 and t be the type? Or just u0? Because it is really well tested that DiffEq only casts to types that the user gives, and I have seen people post an issue about this where they only change the ODE state but keep a time type just Float64. Those two types are kept separate because for example u0 can be complex while t cannot. So if you want to do something to trace position, you should put the tracer to both u0 and t and that should work just fine (and is really well tested with things like AD)

13

Hm, as I said I don’t have experience with DifferentialEquations.jl, but as far as I can tell and debug, every variable is a PrecisionCarrier.

I just took the OP’s code and put precify on everything:

using DifferentialEquations
using Plots
using PrecisionCarriers

function solve_equations()
    kavlength = 599
    k = precify(Vector(logrange(1, 10^4, kavlength + 1)))
    kav = precify((k[1:end - 1] + k[2:end]) / 2)
    iconds = precify(initialConditions(kav, kavlength))
    tspan = precify((0.0, 1.5))
    teval = range(tspan[1], tspan[2], 501) # already a precified range from tspan
    ddrho0 = precify(zeros(kavlength))
    omega2 = precify(zeros(kavlength))

    prob = ODEProblem(diffEquations!, iconds, tspan, (kav, kavlength, ddrho0, omega2), saveat=teval)
    sol = solve(prob, Vern9(), abstol=precify(1e-15), reltol=precify(1e-15))
    # savefig(heatmap(sol.t, kav, log10.(abs.(sol[6:5 + kavlength, :])), yscale=:log10, dpi=600, title="\\delta\\rho_0", xlabel="t", ylabel="k"), "heatmap.svg")
end

function initialConditions(kav, kavlength)
    rho0 = ones(kavlength) ./ (2 .* sqrt.(16 .+ kav.^2))
    drho0 = zeros(size(rho0))
    arr = [15, -0.8, 1, 0.001, 0]
    initial = vcat(arr, rho0, drho0)
    return initial
end

function ddrho0!(a, H, omega2, rho0, drho0, ddrho0)
    for i in eachindex(omega2)
        rho2 = (0.25 + omega2[i] * rho0[i]^2 + drho0[i]^2 / 4) / (2 * rho0[i])
        ddrho0[i] = 4 * rho2 / a^2 - H * drho0[i] - 4 * omega2[i] * rho0[i] / a^2
    end
    nothing
end

function diffEquations!(du, u, args, t)
    kav, kavlength, ddrho0, omega2 = args

    infl, dinfl, a, sigma, dsigma = @view u[1:5]
    rho0 = @view u[6:6 + kavlength - 1]
    drho0 = @view u[6 + kavlength:6 + 2 * kavlength - 1]

    H = sqrt((dinfl^2 + infl^2) / 6)
    da = a * H
    R = dinfl^2 - 2 * infl^2
    
    m2Eff = -a^2 * R / 30
    @. omega2 = kav^2 + m2Eff
    ddrho0!(a, H, omega2, rho0, drho0, ddrho0)
    ddinfl = -3 * H * dinfl - infl
    ddsigma = - H * dsigma - m2Eff * sigma / a^2

    du[1] = dinfl
    du[2] = ddinfl
    du[3] = da
    du[4] = dsigma
    du[5] = ddsigma
    du[6:5 + kavlength] .= drho0
    du[6 + kavlength:5 + 2 * kavlength] .= ddrho0
    nothing
end

solve_equations()

If I do the evil thing and force a convert(Float64, PrecisionCarrier) to return a PrecisionCarrier again, I get the following error, because somewhere it seems to cache values into a Vector{Float64}:

TypeError: in typeassert, expected Float64, got a value of type PrecisionCarrier{Float64}

Stacktrace:
  [1] setindex!(A::Vector{Float64}, x::PrecisionCarrier{Float64}, i::Int64)
    @ Base ./array.jl:987
  [2] macro expansion
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/calculate_residuals.jl:107 [inlined]
  [3] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [4] calculate_residuals!
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/calculate_residuals.jl:106 [inlined]
  [5] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, true, Vector{PrecisionCarrier{Float64}}, Nothing, PrecisionCarrier{Float64}, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, Vector{Vector{PrecisionCarrier{Float64}}}, ODESolution{PrecisionCarrier{Float64}, 2, Vector{Vector{PrecisionCarrier{Float64}}}, Nothing, Nothing, Vector{PrecisionCarrier{Float64}}, Vector{Vector{Vector{PrecisionCarrier{Float64}}}}, Nothing, ODEProblem{Vector{PrecisionCarrier{Float64}}, Tuple{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}}, true, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{saveat::LinRange{PrecisionCarrier{Float64}, Int64}}, SciMLBase.StandardODEProblem}, Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Vector{PrecisionCarrier{Float64}}}, Vector{PrecisionCarrier{Float64}}, Vector{Vector{Vector{PrecisionCarrier{Float64}}}}, Nothing, OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.DEOptions{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{PrecisionCarrier{Float64}, DataStructures.FasterForward}, DataStructures.BinaryHeap{PrecisionCarrier{Float64}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, LinRange{PrecisionCarrier{Float64}, Int64}, Tuple{}}, Nothing, Float64, Nothing, OrdinaryDiffEqCore.DefaultInit, Nothing}, cache::OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, repeat_step::Bool)
    @ OrdinaryDiffEqVerner ~/.julia/packages/OrdinaryDiffEqVerner/rIAEX/src/verner_rk_perform_step.jl:1152
  [6] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEqVerner/rIAEX/src/verner_rk_perform_step.jl:1043 [inlined]
  [7] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, true, Vector{PrecisionCarrier{Float64}}, Nothing, PrecisionCarrier{Float64}, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, Vector{Vector{PrecisionCarrier{Float64}}}, ODESolution{PrecisionCarrier{Float64}, 2, Vector{Vector{PrecisionCarrier{Float64}}}, Nothing, Nothing, Vector{PrecisionCarrier{Float64}}, Vector{Vector{Vector{PrecisionCarrier{Float64}}}}, Nothing, ODEProblem{Vector{PrecisionCarrier{Float64}}, Tuple{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}}, true, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{saveat::LinRange{PrecisionCarrier{Float64}, Int64}}, SciMLBase.StandardODEProblem}, Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Vector{PrecisionCarrier{Float64}}}, Vector{PrecisionCarrier{Float64}}, Vector{Vector{Vector{PrecisionCarrier{Float64}}}}, Nothing, OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqVerner.Vern9Cache{Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.DEOptions{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PrecisionCarrier{Float64}, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{PrecisionCarrier{Float64}, DataStructures.FasterForward}, DataStructures.BinaryHeap{PrecisionCarrier{Float64}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, LinRange{PrecisionCarrier{Float64}, Int64}, Tuple{}}, Nothing, Float64, Nothing, OrdinaryDiffEqCore.DefaultInit, Nothing})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/zs1s7/src/solve.jl:620
  [8] #__solve#62
    @ ~/.julia/packages/OrdinaryDiffEqCore/zs1s7/src/solve.jl:7 [inlined]
  [9] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/zs1s7/src/solve.jl:1 [inlined]
 [10] solve_call(_prob::ODEProblem{Vector{PrecisionCarrier{Float64}}, Tuple{PrecisionCarrier{Float64}, PrecisionCarrier{Float64}}, true, Tuple{Vector{PrecisionCarrier{Float64}}, Int64, Vector{PrecisionCarrier{Float64}}, Vector{PrecisionCarrier{Float64}}}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(diffEquations!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{saveat::LinRange{PrecisionCarrier{Float64}, Int64}}, SciMLBase.StandardODEProblem}, args::Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{abstol::PrecisionCarrier{Float64}, reltol::PrecisionCarrier{Float64}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:667
 [11] solve_call
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:624 [inlined]
 [12] #solve_up#45
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:1199 [inlined]
 [13] solve_up
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:1177 [inlined]
 [14] #solve#43
    @ ~/.julia/packages/DiffEqBase/yhgdI/src/solve.jl:1089 [inlined]
 [15] solve_equations()
    @ Main ~/repos/PrecisionCarriers.jl/notebooks/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W0sZmlsZQ==.jl:16
 [16] top-level scope
    @ ~/repos/PrecisionCarriers.jl/notebooks/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W0sZmlsZQ==.jl:64

I realize that convert should always return the requested type, so I patched that behavior out of my package now. With the “normal” convert behavior, it now runs but reports basically zero epsilons everywhere, which I find hard to believe. Also, it still apparently saves a Vector{Float64} somewhere.

You can get the crashing behavior again by changing the first convert overload in src/functionality/conversion.jl to

function Base.convert(::Type{T}, p::P) where {T <: AbstractFloat}
    return P{T}(convert(T, p.x), p.big)
end