# If almost all ODE solvers fail, how can we trust the results from the other solvers?
**URL:** https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299
**Category:** General Usage
**Tags:** question
**Created:** [January 30, 2021, 9:48pm UTC](https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299 "2021-01-30T21:48:35Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Xiaodong\_Liang](https://discourse.julialang.org/u/Xiaodong_Liang)
#### Post date: [January 30, 2021, 9:48pm UTC](https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299/1 "2021-01-30T21:48:35Z")
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).
```julia
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.
---
### Author:  [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [January 30, 2021, 10:12pm UTC](https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299/2 "2021-01-30T22:12:26Z")
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…
```julia
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.]
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 31, 2021, 1:27am UTC](https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299/4 "2021-01-31T01:27:58Z")
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.
```julia
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.
> [@Xiaodong\_Liang](#):
>
> CVODE\_BDF()
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 😉 .
---
### Author:  [@Xiaodong\_Liang](https://discourse.julialang.org/u/Xiaodong_Liang)
#### Post date: [January 31, 2021, 7:45am UTC](https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299/5 "2021-01-31T07:45:12Z")
Thanks. Yes, it is a problem of parabolic equations with diffusion and a reaction term.
---
### Author:  [@Xiaodong\_Liang](https://discourse.julialang.org/u/Xiaodong_Liang)
#### Post date: [January 31, 2021, 7:49am UTC](https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299/6 "2021-01-31T07:49:00Z")
Hi, I did not withdraw this post, and I actually get very good answers/suggestions/advices.
---
### Author:  [@vavrines](https://discourse.julialang.org/u/vavrines)
#### Post date: [January 31, 2021, 8:02am UTC](https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299/7 "2021-01-31T08:02:57Z")
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.
---
### Author:  [@Xiaodong\_Liang](https://discourse.julialang.org/u/Xiaodong_Liang)
#### Post date: [January 31, 2021, 8:13am UTC](https://discourse.julialang.org/t/if-almost-all-ode-solvers-fail-how-can-we-trust-the-results-from-the-other-solvers/54299/8 "2021-01-31T08:13:32Z")
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.