If almost all ODE solvers fail, how can we trust the results from the other solvers?

General Usage
4

Looks like it’s just a tolerance thing. In MATLAB as well (and in general), tolerance can effect stability. Here, the exp definitely amplifies that. But with the right tolerances it looks like every algorithm is fine.

This thing takes about 0.3 seconds to solve. Simplest seems to be CVODE, swap it to GMRES, 0.4 seconds. And then optimizing the RHS a bit brings it to 0.3 seconds.

using DifferentialEquations, Plots, Sundials

function dyth_dtau(dyth, yth, par, t)
    nt = Int(length(yth)/2)
    dz = par[1] # this can be calculated from n (nt)
    gamma = par[2]
    beta = par[3]
    phi = par[4]
    Le = par[5]
    theta = yth[1:nt]
    y = yth[nt+1:2nt]
    rhs = [phi^2 * y[i] * exp(gamma*(1.0 - 1.0/theta[i])) for i in 1:nt]
    # dtheta/dtau
    dyth[1] = 2 * (theta[2] - theta[1])/dz^2 + beta * rhs[1]
    for i in 2:nt-1
        dyth[i] = (theta[i+1] - 2 * theta[i] + theta[i-1])/dz^2 + beta * rhs[i]
    end
    dyth[nt] = (1 - 2 * theta[nt] + theta[nt-1])/dz^2 + beta * rhs[nt]
    dyth[1:nt] ./= Le
    # dy/dtau
    dyth[nt+1] = 2 * (y[2] - y[1])/dz^2 - rhs[1]
    for i in 2:nt-1
        dyth[nt+i] = (y[i+1] - 2 * y[i] + y[i-1])/dz^2 - rhs[i]
    end
    dyth[2*nt] = (1 - 2 * y[nt] + y[nt-1])/dz^2 - rhs[nt]
end

n = 99
dz = 1/(n+1)
gamma = 30.0
beta = 0.15
phi = 1.1
Le_arr = [1.0, 0.5, 0.25, 0.15, 0.1]
t_span = (0.0, 1.0)
yth0 = ones(2*n)
alg = Rosenbrock23()

#alg = TRBDF2()
#alg = Rados5()
#alg = KenCarp4()
#alg = KenCarp47()
#alg = CVODE_BDF()
#alg = ROCK2()
#alg = RadauIIA5()

for i in 1:length(Le_arr)
    Le = Le_arr[i]
    println(Le)
    par = (dz, gamma, beta, phi, Le)
    prob = ODEProblem(dyth_dtau, yth0, t_span, par)
    sol = solve(prob,alg,abstol=1e-8,reltol=1e-6)

    if i == 1
        global plt=plot(sol.t, sol[1, :])
    else
        global plt=plot!(sol.t, sol[1, :])
    end
end

# Accelerate the user's dydt function because it's very unoptimized!

using ModelingToolkit
par = (dz, gamma, beta, phi, Le_arr[end])
prob = ODEProblem(dyth_dtau, yth0, t_span, par)
sys = modelingtoolkitize(prob)
fastprob = ODEProblem(sys,yth0,t_span,par,jac=true,sparse=true)

@time sol = solve(fastprob,CVODE_BDF(linear_solver=:GMRES))
@time sol = solve(fastprob,Rosenbrock23())
@time sol = solve(prob,RadauIIA5())
@time sol = solve(fastprob,KenCarp47())
@time sol = solve(fastprob,Rodas5(),abstol=1e-8,reltol=1e-6)

The runner up seems to be Rosenbrock23() with sparsity and the RHS rewritten (the RHS function is very clearly unoptimal).

I might want to take a deeper look later (or @YingboMa might) to see why SDIRK methods in particular seem to require lower tolerances here. My guess is it has to do with tolerance of the nonlinear solve. I’d be curious to see what MATLAB’s ode23tb acts like on this.

Quoting the documentation:

Note that this setup is not automatically included with DifferentialEquations.jl. To use the following algorithms, you must install and use Sundials.jl:

]add Sundials
using Sundials

And it’s RadauIIA5(). I found the spot in the docs you must’ve been looking at and fixed the typo. All other instances than 1 spot had it right, so I can see what list you were looking at :wink: .

8 Likes