
A Comparison of Differential Equation Solver Suites - one-more-minute
http://www.stochasticlifestyle.com/comparison-differential-equation-solver-suites-matlab-r-julia-python-c-fortran/
======
ViralBShah
Chris is an encyclopedia on ODEs. Take a look at this. No other system has
this thorough a study of so many methods in one place.

[https://github.com/JuliaDiffEq/DiffEqDevTools.jl/blob/master...](https://github.com/JuliaDiffEq/DiffEqDevTools.jl/blob/master/src/ode_tableaus.jl)

This is the kind of thing we imagined experts being able to do when we
designed Julia. Having implemented some such solvers in the past, I can safely
say this stuff is easily the best out there.

------
3JPLW
Chris — the author of this blog post and lead developer of the
DifferentialEquations.jl ecosystem — had a really nice workshop at this year's
JuliaCon. You can watch it online here:
[https://www.youtube.com/watch?v=75SCMIRlNXM](https://www.youtube.com/watch?v=75SCMIRlNXM).
His intro there is really approachable and even if you just watch the first 10
minutes you can come away with something tangible.

------
scotto3394
I see those Bombay Sapphires have been working towards a good cause!

But in all seriousness, the comparison was far more thorough than I expected
going in. I learned quite a bit, and it's making me consider giving Julia
another shot. I've been tinkering with mostly Rust lately, very happily, but
this post makes me miss something with the strong scientific bend of Julia. So
job well done and thanks for taking the time to write this.

------
throwaway613834
Intel has an ODE solver [1] that seems _unbelievably_ good at solving certain
kinds of differential equations. I'm surprised/disappointed it wasn't even
mentioned. I would've loved to see an evaluation of it against MATLAB (or
other solvers) more than any other pair in the list.

[1] [https://software.intel.com/en-us/articles/intel-ordinary-
dif...](https://software.intel.com/en-us/articles/intel-ordinary-differential-
equations-solver-library)

~~~
keldaris
Could you elaborate on what's "unbelievably good" about it? From a quick
glance at the site and documentation I don't see any particularly
distinguishing features.

~~~
throwaway613834
Sure. Digging through my old emails, here's the particular example (a very
stiff Van der Pol equation) I was looking at, taken from Intel's own examples,
which for MATLAB would be:

    
    
      tic;
      [t, x] = ode15s(
        @(t, x) [x(2); ((1 - x(1)^2) * x(2) - x(1)) * 10^6],
        [0, 160], [2; 0], odeset('RelTol', 1e-5));
      toc;
      x(size(x, 1), :)
    

MATLAB takes 7.2 seconds to output the result [1.872, -0.748].

Intel's dodesol_rkm9mkn() gave me [1.877, -0.743], _which is a more accurate
result_ , in less than the blink of an eye: in _0.035 seconds_.

That's _more_ than 200x faster, by mere virtue of using a better algorithm.
Apparently the algorithm is "Merson's method", which I'm completely unfamiliar
with and which I'd never heard of until then.

~~~
ChrisRackauckas
I went ahead and ran some stuff on my laptop.

    
    
        using DifferentialEquations
    
        van = @ode_def_noinvhes VanDerPol2 begin
          dy = μ*((1-x^2)*y - x)
          dx = 1*y
        end μ=>1.
    
        prob = ODEProblem(VanDerPol2(μ=1e6),[0;2.],(0.0,160.0))
      
      
      solve(prob,CVODE_BDF(),reltol=1/10^5,save_everystep=false)
       solve(prob,Rodas5(),reltol=1/10^5,save_everystep=false)
    
    
        @time sol1 = solve(prob,CVODE_BDF(),reltol=1/10^5,save_everystep=false)
        @time sol2 = solve(prob,Rodas5(),reltol=1/10^5,save_everystep=false)
        sol3 = solve(prob,CVODE_BDF(),abstol=1/10^10,reltol=1/10^10)
    

On my laptop I'm getting ~0.1 seconds for both, with Sundials being a little
less accurate and OrdinaryDiffEq.jl being a little more accurate against the
10^-10 reference solution than the reported Intel. But "spot checking" like
this makes it hard to just timings without issues of more/less error involved,
which is why more sophisticated benchmarks use work-precision diagrams. So the
timings are in the same ballpark as Intel and I'd assume my (not so fast)
laptop doesn't match what you benchmarked with, the point being that MATLAB is
just slow here instead of Intel being unbelievably exceptional (though
admirable). Still, maybe Intel does pull off some extra efficiency here, but
then of course no method is universally good so I wouldn't expect that
everywhere. In the end, it's probably a well optimized (probably in the
DiffEq.jl/Sundials/FATODE range) single method without event handling or other
extra features, but importantly it is abandoned. Probably would've been
interesting if development kept up, but without it I think I'm just not going
to try to get it installed.

Thanks for pointing me to it though. It's definitely an interesting nugget of
ODE history that I didn't know existed.

~~~
throwaway613834
> Probably would've been interesting if development kept up, but without it I
> think I'm just not going to try to get it installed.

Out of curiosity, why is this an issue? Wouldn't you have just downloaded the
library and run it regardless of whether someone is actively developing it or
not? I don't really understand this argument (and I see many people making it,
not just you).

~~~
ChrisRackauckas
Well vanilla ODE solvers are a dime a dozen. Sufficiently optimized stiff ODE
solvers tend to get around the same efficiency, with some methods doing better
on some problems, and others doing better on other problems. So it's all
problem dependent and this is just one implicit RK method, while things like
ARKODE already have quite a few, FATODE has quite a few, and
DifferentialEquations.jl has quite a few, etc, meaning I wouldn't expect it to
be too different.

What would make it exciting to play with if it added something. Free
parallelism or uses the Xeon Phi automatically, complex number support, event
handling (many examples where this is necessary. One: PDE solvers need this
for adding projection steps). But since development has halted, we're pretty
much guaranteed that none of this will ever come.

------
ivan_ah
For a more "systems" approach for ODEs, you should look at Modelica:
[https://en.wikipedia.org/wiki/Modelica](https://en.wikipedia.org/wiki/Modelica)
It's very powerful stuff.

~~~
ihnorton
Julia's meta-programming functionality is ideal for implementing such domain-
specific languages (the popular JuMP suite is a flagship example in a
different domain:
[https://jump.readthedocs.io/en/latest/](https://jump.readthedocs.io/en/latest/)).

Several groups are working on systems inspired by Modelica, including:

[http://tshort.github.io/Sims.jl/latest/](http://tshort.github.io/Sims.jl/latest/)

and

[https://www.modelica.org/events/modelica2017/proceedings/htm...](https://www.modelica.org/events/modelica2017/proceedings/html/submissions/ecp17132693_ElmqvistHenningssonOtter.pdf)
([https://github.com/ModiaSim/Modia.jl](https://github.com/ModiaSim/Modia.jl))

------
Myrmornis
Can anyone recommend any introductory material on real-world examples of
differential equations? I don't mean the first-year undergrad textbook stuff
with differential equations for population growth and chemical reactions and
mechanics that can be solved via separation of variables etc. I mean
descriptions of models used in science and engineering that result in needing
to use the fancy solvers that this thread is discussing.

~~~
snaky
Electronics is all about differential equations.

> Circuit simulation programs take a text netlist describing the circuit
> elements (transistors, resistors, capacitors, etc.) and their connections,
> and translate this description into equations to be solved. The general
> equations produced are nonlinear differential algebraic equations

[https://en.wikipedia.org/wiki/SPICE](https://en.wikipedia.org/wiki/SPICE)

All the way down to transistors and other basic elements

> These models are very computer intensive, involving detailed spatial and
> temporal solutions of coupled partial differential equations on three-
> dimensional grids inside the device

[https://en.wikipedia.org/wiki/Transistor_model](https://en.wikipedia.org/wiki/Transistor_model)

------
alimw
This seems to be all about ODEs. Has anyone got any recommendation for a good
PDE solver? One that might run a problem having more than three spatial
dimensions? (and ideally callable from Python)

~~~
ChrisRackauckas
Talking about PDEs is way more complex and I didn't get into it. What do you
mean by a good PDE solver? Generally a PDE solver is built by doing some
spatial discretization and then either solving that with some numerical linear
solver directly (if it's a time-independent problem, or if it's a fully-
implicit discretization), or using an ODE solver. There are many different
spatial discretization techniques, including finite element, finite
difference, finite volume, and spectral methods. There are libraries for
helping one discretize space using each of these methods.

But then there are "full-stop PDE solvers" which choose a specific PDE and
solve that. For example, MATLAB and Mathematica can numerically solve
diffusion-advection systems directly. And there are many user packages for
different PDEs.

So to cover this fully, I would need to go through the available tools for
numerical linear algebra and spatial discretization, but on the other hand
cover some full-stop PDE solvers as well. This will be a different blog post.
Without knowing your background or what you want to do, FEniCS is a good
Python finite element toolbox that can then be used to build the linear
equation or system of ODEs, from which you use other tools like those
mentioned here to complete the solving. As for full-stop PDE solvers in
Python, I am not sure and need to look into it myself...

~~~
alimw
Thanks! My problem is based on a straightforward constant-coefficient
diffusion-advection equation, but there are complicated boundaries and I'd
like it to work in six or seven dimensions (given time and a big computer!). I
have read a little about FEniCS but I think it will struggle just to generate
the mesh. I'm looking at Mathematica now. It seems like it will treat an
evolution equation specifically (method of lines)---perhaps that's what you
mean by a "full-stop solver"?

~~~
ChrisRackauckas
Yup, that's exactly what I mean. People have pre-written solvers for specific
equations. Sounds like you found one for you, so that's great! When you ask
about PDE tools though, people start talking about tools for building PDE
solvers like FEniCS becuase it's so rare to find something exactly for your
PDE, but there's no reason to use these tools unless you have to.

------
jxcole
It's a pity there are no JavaScript options listed.

~~~
tlarkworthy
Wouldn't it be a single entry, numeric.js?

~~~
ChrisRackauckas
Yeah, I'm not sure of any others. But looking at numeric.js, it has a single
ODE solver which is Dopri, and thus matches the tableau of ode45 and dopri5.
However, it's easy to see that their implementation is very different just by
this line:

    
    
        h = min(0.8*h*pow(tol/erinf,0.25),4*h);
    

That means that the timesteps are changes via Cechino's optimality theorem,
the "standard" timestepping method which isn't actually used in production
quality solvers. The reason is because this form of adaptivity causes solvers
to change their h too rapidly, which actually decreases the stability region
and increases the number of step rejections. So all of the other explicit RK
solvers mentioned (except when I single it out) all use a PI-controller. And
this makes a pretty big difference: I've found something like a 10x
performance change by using a PI-controller on many equations, and then
there's many other equations where at a "simple" tolerance like reltol=10^-3
the dopri-type method will diverge without a PI-controller on the timesteps.
But it does have event handling and dense output. However, since it's lacking
any solvers for stiff ODEs, I am not sure it really is a full suite.

------
amelius
I see no mention of multigrid methods.

~~~
ChrisRackauckas
Multigrid methods are linear solvers. Generally they are only not used in the
context of timestepping. PDEs will be a separate post.

~~~
remcob
I'm so looking forward to this post!. hp-FEM methods, in particular, caught my
attention, they look so promising. Do they work in practice?

Some time ago I was trying to make a 'hpq-FEM' using rational functions as a
basis, like a Padé approximation. Do you know if anyone has tried something
similar?

