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

General Usage
1

I want to introduce Julia to one course. We have been using MATLAB in this course for many years, which are very robust. For these special problems, I tried Python and it seems not as robust as MATLAB (I am also a beginner of Python), but I managed to solve them all. I heard from friends about Julia, and I tried, but apparently I can hardly solve the same problems (the warning/error is “Instability detected. Aborting”). Here below is the source code, and it seems only the solver Rosenbrock23 gives the expected results. This situation makes it difficult for me to explain to others (in MATLAB almost all solvers work, just the efficiency matters).

 using DifferentialEquations, Plots

function dyth_dtau(dyth, yth, par, t)    
    # nt = int(length(yth)/2)
    nt = par[1] # this can be calculated from the line given above
    dz = par[2] # this can be calculated from n (nt)
    gamma = par[3]
    beta = par[4]
    phi = par[5]
    Le = par[6]
    theta = yth[1:nt]
    y = yth[nt+1:2*nt]
    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] = 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)

for i in 1:length(Le_arr)
    Le = Le_arr[i]
    println(Le)
    par = (n, dz, gamma, beta, phi, Le)
    prob = ODEProblem(dyth_dtau, yth0, t_span, par)
    sol = solve(prob)
    if (i == 1)
        global plt=plot(sol.t, sol[1, :])
    else
        global plt=plot!(sol.t, sol[1, :])
    end
end

display(plt)

Here is the status of the package: " [0c46a032] DifferentialEquations v6.16.0 ".

It is very strange that not all ODE solvers are available. For example: CVODE_BDF() and RadauIIA() are missing.

Could you please give me some advices? Thank you.

2 Likes
2

I’ve had very good experience with the DifferentialEquations.jl solvers – when I learned to use them, so I’m a little bit surprised by your experience.

Your code would be much easier on the eye if you insert the code within a start line consisting of triple backticks followed by “julia”, and and a closing line consisting of triple backticks. Like…

using DifferentialEquations, Plots

function dyth_dtau(dyth, yth, par, t)
# nt = int(length(yth)/2)
nt = par[1] # this can be calculated from the line given above
dz = par[2] # this can be calculated from n (nt)
gamma = par[3]
beta = par[4]
phi = par[5]
Le = par[6]
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] = 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

[You model looks like a discretized set of 2 PDEs, e.g., a simultaneous temperature distribution (theta) and possibly concentration (y), with diffusion and some reaction term.]

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
5

Thanks. Yes, it is a problem of parabolic equations with diffusion and a reaction term.

6

Hi, I did not withdraw this post, and I actually get very good answers/suggestions/advices.

7

Never mind, the info appears since I withdraw my post.
The example code you post doesn’t run with gamma(1.0 - 1.0/theta[i]), I thought you need a Gamma function, but it seems to be par[3]. You just need to add a * between multiplicators, so I recalled that post.
For the algorithm usage, please refer Chris’ answer.

8

Hi Chris, thank you so much. I am completely new to Julia, and several colleagues of mine are very much excited about Julia, and highly recommend to take it as the first language for research and teaching (we have been using FORTRAN + MATLAB for quite some years). The first step is to introduce it in one course for which we have been using MATLAB. Thank you. Everything is working now.

3 Likes