Stream: helpdesk (published)

Topic: Tricks for solving an ODE with a large derivative?

Robbie Rosati (Mar 05 2021 at 01:55):

I'm trying to solve an ODE with DifferentialEquations.jl

$$\mathbf{\Gamma^\prime} = \mathbf{u}(t) \mathbf{\Gamma}$$

where \Gamma and u are matrices. u is a block matrix, and has an element or two that are large enough that Tsit5() and Rodas5() overflow using Float64s. Using BigFloats slows down the solvers to a crawl, and I haven't actually gotten the equation of motion to solve.
Are there any tricks for speeding this up and/or avoiding overflows? If these were scalars I would just solve the log of this equation, but matrix logs need commutativity to apply the usual identity.

Here's some working and probably badly optimized code:

using DifferentialEquations
using LinearAlgebra
using LambertW

const d = 3
const Tend = 75.0
const Tstart = Tend - 65.0
const k1 = 10.0
const L = 2/k1
const λ = 0.7
const ϕ₀ = 20.0
const πv = -3L
const H₀ = 1/(3-Eh(0.0))

function Eh(t)
    return 3L*λ/(2+L*λ)
end
const πw = sqrt(L*λ*(3-Eh(0.0)) - 9L^2)

function g(t)
    f = ϕ(t)
    return Diagonal([1,exp(k1*f),exp(k1*f)])
end

function V(t)
    return exp(λ*ϕ(t))
end

function H(t)
    return sqrt(V(t)/(3-Eh(t)))
end

function ϕ(t)
    return ϕ₀-2/λ*lambertw(3λ*L*t/(2*sqrt(3-Eh(t)))*exp(λ*ϕ₀/2))
end

function Mab(t)
    e = Eh(t)
    Vab = Diagonal((3-e).*[λ^2,k1*λ/2,k1*λ/2])
    piu = [πv, πw, 0]
    pil = g(t)*piu
    gvl = [λ*(3-e),0.0,0.0]
    gvu = inv(g(t))*gvl
    Rππ = (piu'.*pil - I*2*e)/L^2
    Pot = (3-e)*piu'.*pil + piu'.*gvl + gvu'.*pil
    return Vab + Rππ + Pot
end

function hp(t)
    return [ 0 exp(k1*ϕ(t))*k1*πw/2 0;
            k1*πw/2  k1*πv/2 0;
            0 0 k1*πv/2]
end

const loga = -127.0
function Gprime(dG,G,logk,t)
    @inbounds begin
        T = eltype(G)
        u = Array{T}(undef,2*d,2*d)
        u[1:d,1:d] .= zeros(T,(d,d))
        u[1:d,d+1:end] .= Matrix{T}(I,d,d)
        u[d+1:end,1:d] .= -Matrix{T}(I,d,d)*exp(2*(logk-t-loga-log(H(t)))) .- Mab(t)
        u[d+1:end,d+1:end] .= Matrix{T}(I,d,d)*(Eh(t)-3)
        Gπ = hp(t)
        G11 = @view G[1:d,:]
        G12 = @view G[d+1:end,:]
        dG .= u*G .- vcat(Gπ*G11, Gπ*G12)
    end
    nothing
end


prob = ODEProblem(Gprime,Matrix{Float64}(I,2d,2d),(Float64(Tstart),Float64(Tend)),-136.0)
solve(prob,Tsit5(),reltol=1e-10,abstol=1e-10,progress=true)

Typically I get retcodes of Unstable or DtLessThanMin.

Mason Protter (Mar 05 2021 at 02:02):

@Christopher Rackauckas if he's around is probably the ideal person to comment, but I can try to give some very generic advice.

Mason Protter (Mar 05 2021 at 02:02):

First off, are you sure that your differential equation is correct? E.g. are there some consistency checks you can do?

Christopher Rackauckas (Mar 05 2021 at 02:03):

Yes, first make sure it's actually the equation you wanted to solve. That's 99.9% of issues.

Christopher Rackauckas (Mar 05 2021 at 02:04):

If it does solve with BigFloat but not Float64, then you might want to look into alternative number types. ArbFloats.jl, DoubleFloats.jl, etc.

Christopher Rackauckas (Mar 05 2021 at 02:04):

There's a large set of higher precision number types that are faster than the generic BigFloat.

Mason Protter (Mar 05 2021 at 02:05):

It may also be that there's a clever way you can rescale your problem so that you can just use Float64

Mason Protter (Mar 05 2021 at 02:06):

E.g. a very simple example would be that if the solution is of the form

f(x) = 1e100 + g(x)

then you can often just subtract off the big part and solve for g and get more stable numerics.

Christopher Rackauckas (Mar 05 2021 at 02:08):

Using a non-generic ODE solver could help too. Defining the problem with a DiffEqArrayOperator and using the Magnus integrators could help too https://diffeq.sciml.ai/stable/solvers/nonautonomous_linear_ode/#State-Independent-Solvers

Mason Protter (Mar 05 2021 at 02:08):

But seriously, there's been so many times where I was tryng to solve an equation and getting a solution that oscillated like crazy or was super stiff or blew up in other ways and like the fourth time I double checked the equation I realized I had made a mistake defining the diffeq.

Robbie Rosati (Mar 05 2021 at 02:09):

Mason Protter said:

It may also be that there's a clever way you can rescale your problem so that you can just use Float64

Yeah, if this equation is correct that's my hope. I'll triple check it and make sure.

Mason Protter (Mar 05 2021 at 02:19):

One trick you can use is try to do an asymptotic expansion on your function to see various behaviour and factor it out. E.g. you may be able to show that for large $x$, your solution goes as

$$f(x) \approx e^{-x}$$

in which case, you might be able to write

$$f(x) = e^{-x} g(x)$$

then you can stick this form of $f$ into your differential equation and solve for $g$ instead.

Andrew Dolgert (Mar 05 2021 at 04:34):

@Robbie Rosati Have you estimated the condition number of the u-matrix? If some numbers are very large, then you may have two problems: both overflow and stiffness. Mason suggests factorization can help. If it isn't stiff, then, yes, rescale the whole thing, but it's most likely stiff. There are some other ways to take Mason's advice: Picard iterations will approximate an integral over an interval. Operator splitting will let you combine a stable integration with an unstable one. Basically, find that condition number, and it may send you to a chapter on stiff equations.

Robbie Rosati (Mar 05 2021 at 17:44):

Thanks everyone for the replies. So, some updates:

• I rechecked the example I made and did find a mistake, but it was in a subleading term and doesn't seem to change the behavior much.
• I changed the example to a 4x4 version (as small as I can go), also without much change qualitative change.
• @Andrew Dolgert The $u$ matrix has a huge condition number $||u||*||u^{-1}|| \sim 10^{19}$ (and shrinking as time evolves)
• I can get it to solve with with BigFloats or DoubleFloats.jl if I increase Tend. This also decreases the condition number by a lot.
• I tried changing the solver to a stiff oneRosenbrock23(autodiff=false) (with the original Tend) and I get DtLessThanMin with Float64s. It seems to be working with Double64s but is very slow. Is there a better solver to try, or is there a better way to massage the equation?

I think there's a paper working on a similar model that estimated the growth of the solution, I might try @Mason Protter 's suggestion and see if I can estimate the growth and factor it out.

Mason Protter (Mar 05 2021 at 17:52):

You may want to give https://github.com/dzhang314/MultiFloats.jl a try. Should be about 2x faster than Double64, but with worse NaN handling.

Mason Protter (Mar 05 2021 at 17:53):

Though they don't support transcendental functions yet, so that might be a dealbreaker.

Last updated: Oct 02 2023 at 04:34 UTC