8 ms·

A Comparison of Differential Equation Solver Suites

3JPLW 9y ago
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.
throwaway613834 9y ago
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-differential-equations-solver-library https://software.intel.com/en-us/articles/intel-ordinary-dif...
- keldaris 9y ago
  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 9y ago
    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.
    - RBerenguel 9y ago
      Do the initials stand for Runge-Kutta-Merson, order 9? I wonder about the other 3 (mkn)
      - throwaway613834 9y ago
        I believe those are the initials, yes. The "9" refers to the maximum number of "stages" which is not the same thing as the order of the solver. The (m,k) refer to the corresponding parameters of the method (I think). The "n" at the end is for "numerical" Jacobian (as opposed to "a" for user-provided/"analytical"). There should be a PDF in the download that explains these better than I possibly can (though it's still somewhat terse).
        RBerenguel 9y ago
        I had a look before asking but could not find anything before the download. Thanks, it piqued my curiosity (used to write & theoretically analyse several numerical integrators back in my teaching days)
    - ChrisRackauckas 9y ago
      I don't think this really backs up your statement though. I am not at all surprised by this result. In the post I mentioned that MATLAB wasn't really designed for efficiency, and in this repo I show how even with explicit Runge-Kutta methods (i.e. the exact same methods, and these methods are simple enough that there isn't too much variation in implementation details sans getting it all to SIMD/FMA/etc.) you can get orders of magnitude better than MATLAB: https://github.com/JuliaDiffEq/MATLABDiffEq.jl https://github.com/JuliaDiffEq/MATLABDiffEq.jl So since stiff solvers have a lot more moving parts and thus tend to have even larger differences, 200x is easily explained by just language differences. One way to test this is to instead test the Intel solver against Sundials. I wouldn't be surprised if Sundials actually was faster in this case. But there's another factor here. I explained in the blog post that multistep methods like BDF tend to do well on asymtopically large problems since they minimize the amount of user-function evaluations, but they don't tend to do well on smaller systems of ODEs. This is exactly that case: you're comparing an implicit RK (Intel) to MATLAB's BDF method, and see that the RK is more accurate and is fast. But this effect of the method choice on small (user-scale) problems is exactly why I emphasized it as a downside when solver suites only include multistep BDF methods as their stiff solvers. That isn't to say there isn't interesting here. Intel's nonstiff method is a 4th order method that can add extra 1st order stages to boost the stability region, similar to a Runge-Kutta Chevyshev method. I honestly have no idea what this stiff solver is so I bookmarked it to read when I get to a campus where the paper is free: http://www.sciencedirect.com/science/article/pii/0041555389900578 http://www.sciencedirect.com/science/article/pii/00415553899... But in general, RKC methods weren't widely adopted because usually using a full non-stiff or stiff solver ends of being better, and picking wildly different methods like this which are not widely known usually isn't a good sign in a very "crowded" field like ODE solvers where generally things are pretty figured out and it's all about the details.
- SloopJon 9y ago
  Intel's ODE solver library appears to be pre-release software, retired in 2011, that was only licensed for "internal testing and development."
  - throwaway613834 9y ago
    That's the ironic thing about it. At least for the examples I've tried (probably a biased subset, but still), it beats the crap out of the couple of commercial solvers I've tried.
- snaky 9y ago
  Might it be somewhat related to the 'Neuromorphic Chip' recently announced by Intel - which is, as discussed here on HN[1], being especially good for approximating dynamic systems and to use symbolic computation? [1] https://news.ycombinator.com/item?id=15334533 https://news.ycombinator.com/item?id=15334533
ViralBShah 9y ago
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/src/ode_tableaus.jl https://github.com/JuliaDiffEq/DiffEqDevTools.jl/blob/master... 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.
- deleted 9y ago
  [deleted]
jxcole 9y ago
It's a pity there are no JavaScript options listed.
- tlarkworthy 9y ago
  Wouldn't it be a single entry, numeric.js?
  - ChrisRackauckas 9y ago
    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.
alimw 9y ago
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 9y ago
  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 9y ago
    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"?
    - 9y ago
scotto3394 9y ago
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.
amelius 9y ago
I see no mention of multigrid methods.
- ChrisRackauckas 9y ago
  Multigrid methods are linear solvers. Generally they are only not used in the context of timestepping. PDEs will be a separate post.
  - remcob 9y ago
    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?
Myrmornis 9y ago
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 9y ago
  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
ivan_ah 9y ago
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 9y ago
  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/html/submissions/ecp17132693_ElmqvistHenningssonOtter.pdf https://www.modelica.org/events/modelica2017/proceedings/htm... (https://github.com/ModiaSim/Modia.jl https://github.com/ModiaSim/Modia.jl)

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.

throwaway613834 9y ago

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-differential-equations-solver-library https://software.intel.com/en-us/articles/intel-ordinary-dif...

keldaris 9y ago
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 9y ago
  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.
  - RBerenguel 9y ago
    Do the initials stand for Runge-Kutta-Merson, order 9? I wonder about the other 3 (mkn)
    - throwaway613834 9y ago
      I believe those are the initials, yes. The "9" refers to the maximum number of "stages" which is not the same thing as the order of the solver. The (m,k) refer to the corresponding parameters of the method (I think). The "n" at the end is for "numerical" Jacobian (as opposed to "a" for user-provided/"analytical"). There should be a PDF in the download that explains these better than I possibly can (though it's still somewhat terse).
      - RBerenguel 9y ago
        I had a look before asking but could not find anything before the download. Thanks, it piqued my curiosity (used to write & theoretically analyse several numerical integrators back in my teaching days)
  - ChrisRackauckas 9y ago
    I don't think this really backs up your statement though. I am not at all surprised by this result. In the post I mentioned that MATLAB wasn't really designed for efficiency, and in this repo I show how even with explicit Runge-Kutta methods (i.e. the exact same methods, and these methods are simple enough that there isn't too much variation in implementation details sans getting it all to SIMD/FMA/etc.) you can get orders of magnitude better than MATLAB: https://github.com/JuliaDiffEq/MATLABDiffEq.jl https://github.com/JuliaDiffEq/MATLABDiffEq.jl So since stiff solvers have a lot more moving parts and thus tend to have even larger differences, 200x is easily explained by just language differences. One way to test this is to instead test the Intel solver against Sundials. I wouldn't be surprised if Sundials actually was faster in this case. But there's another factor here. I explained in the blog post that multistep methods like BDF tend to do well on asymtopically large problems since they minimize the amount of user-function evaluations, but they don't tend to do well on smaller systems of ODEs. This is exactly that case: you're comparing an implicit RK (Intel) to MATLAB's BDF method, and see that the RK is more accurate and is fast. But this effect of the method choice on small (user-scale) problems is exactly why I emphasized it as a downside when solver suites only include multistep BDF methods as their stiff solvers. That isn't to say there isn't interesting here. Intel's nonstiff method is a 4th order method that can add extra 1st order stages to boost the stability region, similar to a Runge-Kutta Chevyshev method. I honestly have no idea what this stiff solver is so I bookmarked it to read when I get to a campus where the paper is free: http://www.sciencedirect.com/science/article/pii/0041555389900578 http://www.sciencedirect.com/science/article/pii/00415553899... But in general, RKC methods weren't widely adopted because usually using a full non-stiff or stiff solver ends of being better, and picking wildly different methods like this which are not widely known usually isn't a good sign in a very "crowded" field like ODE solvers where generally things are pretty figured out and it's all about the details.
SloopJon 9y ago
Intel's ODE solver library appears to be pre-release software, retired in 2011, that was only licensed for "internal testing and development."
- throwaway613834 9y ago
  That's the ironic thing about it. At least for the examples I've tried (probably a biased subset, but still), it beats the crap out of the couple of commercial solvers I've tried.
snaky 9y ago
Might it be somewhat related to the 'Neuromorphic Chip' recently announced by Intel - which is, as discussed here on HN[1], being especially good for approximating dynamic systems and to use symbolic computation? [1] https://news.ycombinator.com/item?id=15334533 https://news.ycombinator.com/item?id=15334533

ViralBShah 9y ago

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/src/ode_tableaus.jl https://github.com/JuliaDiffEq/DiffEqDevTools.jl/blob/master... 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.

deleted 9y ago
[deleted]

jxcole 9y ago

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

tlarkworthy 9y ago
Wouldn't it be a single entry, numeric.js?
- ChrisRackauckas 9y ago
  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.

alimw 9y ago

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 9y ago
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 9y ago
  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"?
  - 9y ago

scotto3394 9y ago

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.

amelius 9y ago

I see no mention of multigrid methods.

ChrisRackauckas 9y ago
Multigrid methods are linear solvers. Generally they are only not used in the context of timestepping. PDEs will be a separate post.
- remcob 9y ago
  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?

Myrmornis 9y ago

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 9y ago
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

ivan_ah 9y ago

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 9y ago
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/html/submissions/ecp17132693_ElmqvistHenningssonOtter.pdf https://www.modelica.org/events/modelica2017/proceedings/htm... (https://github.com/ModiaSim/Modia.jl https://github.com/ModiaSim/Modia.jl)