Tag Archives: matlab

A Comparison Between Differential Equation Solver Suites In MATLAB, R, Julia, Python, C, Mathematica, Maple, and Fortran

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/comparison-differential-equation-solver-suites-matlab-r-julia-python-c-fortran/

Many times a scientist is choosing a programming language or a software for a specific purpose. For the field of scientific computing, the methods for solving differential equations are what’s important. What I would like to do is take the time to compare and contrast between the most popular offerings.
This is a good way to reflect upon what’s available and find out where there is room for improvement. I hope that by giving you the details for how each suite was put together (and the “why”, as gathered from software publications) you can come to your own conclusion as to which suites are right for you.

(Full disclosure, I am the lead developer of DifferentialEquations.jl. You will see at the end that DifferentialEquations.jl does offer pretty much everything from the other suite combined, but that’s no accident: our software organization came last and we used these suites as a guiding hand for how to design ours.)

Quick Summary Table

If you just want a quick summary, I created a table which has all of this information. You can find it here (click for PDF):

Comparison Of Differential Equation Solver Software

MATLAB’s Built-In Methods

Due to its popularity, let’s start with MATLAB’s built in differential equation solvers. MATLAB’s differential equation solver suite was described in a research paper by its creator Lawerance Shampine, and this paper is one of the most highly cited SIAM Scientific Computing publications. Shampine also had a few other papers at this time developing the idea of a “methods for a problem solving environment” or a PSE. The idea is pretty simple: users of a problem solving environment (the examples from his papers are MATLAB and Maple) do not have the same requirements as more general users of scientific computing. Instead of focusing on efficiency, they key for this group is to have a clear and neatly defined (universal) interface which has a lot of flexibility.

The MATLAB ODE Suite does extremely well at hitting these goals. MATLAB documents its ODE solvers very well, there’s a similar interface for using each of the different methods, and it tells you in a table in which cases you should use the different methods.

But the modifications to the methods goes even further. Let’s take for example the classic ode45. That method just works and creates good plots, right? Well, Shampine added a little trick to it. When you solve an equation using ode45, the Runge-Kutta method uses a “free” interpolation to fill in some extra points. So between any two steps that the solver takes, it automatically adds in 4 extra points using a 4th order interpolation. This is because high order ODE solvers are good enough at achieving “standard user error tolerances” that they actually achieve quite large timesteps, and in doing so step too infrequently to make a good plot. Shampine’s scheme is a good quick fix to this problem which most people probably never knew was occurring under the hood!

There’s quite a bit of flexibility here. The methods allow you to use complex numbers. You’re given access to the “dense output function” (this is the function which computes the interpolations). There’s a few options you can tweak. Every one of these methods is setup with event handling, and there are methods which can handle differential-algebraic equations. There are also dde23 and ddesd for delay differential equations, and in the financial toolbox there’s an Euler-Maruyama method for SDEs.

While MATLAB does an excellent job at giving a large amount of easily available functionality, where it lacks is performance. There’s a few reasons for this. For one thing, these modifications like adding extra points to the solution array can really increase the amount of memory consumed if the ODE system is large! This actually has an impact in other ways. There’s a very good example of this in ode45. ode45 is based on the Dormand-Prince 5(4) pair. However, in 1999, the same year the MATLAB ODE Suite was published, Shampine released a paper with a new 5(4) pair which was more efficient than the Dormand-Prince method. So that begs the question, why wasn’t this used in the MATLAB ODE Suite (it’s clear Shampine new about it!)? (I actually asked him in an email) The reason is because its interpolation requires calculating some extra steps, so it’s less efficient if you are ALWAYS interpolating. But since ode45 is always interpolating in order to make the plots look nicer, this would get in the way. Essentially, it can be more efficient, but MATLAB sets things up for nice plotting and not pure efficiency.

But there are other areas where more efficient methods were passed up during the development phase of the ODE suite. For example, Hairer’s benchmarks in his book Solving Ordinary Differential Equations I and II (the second is for stiff problems), along with the benchmarks from the Julia DifferentialEquations.jl suite, consistently show that high order Runge-Kutta methods are usually the most efficient methods for high accuracy solving of nonstiff ODEs. These benchmarks both consistently show that, for the same error, high order Runge-Kutta methods (like order >6) can solve the equation much faster than methods like Adams methods. But MATLAB does not offer high order Runge-Kutta methods and only offers ode113 (an Adams method) for high-accuracy solving.

Some of this is due to a limitation within MATLAB itself. MATLAB’s ODE solver requires taking in a user-defined function, and since this function is defined in MATLAB its function calls are very inefficient and expensive. Thus MATLAB’s ODE solver suite can become more efficient by using methods which reduce the number of function calls (which multistep methods do). But this isn’t the only case where efficient methods are missing. Both Hairer and the JuliaDiffEq benchmarks show that high-order Rosenbrock methods are the most efficient for low-medium accuracy stiff ODEs, but MATLAB doesn’t offer these methods. It does offer ode23s, a low-order Rosenbrock method, and ode15s which is a multistep method. ode15s is supposed to be the “meat and potatoes of the MATLAB Suite”, but it cannot handle all equations of since its higher order methods (it’s adaptive order) are not L-stable (and not even A-stable). For this reason there’s a few other low order SDIRK methods (ode23tb, an ESDIRK method for highly stiff problems) which are recommended to fill in the gaps, but none of the higher order variants which are known to be more efficient for many equations.

This pattern goes beyond the ODE solvers. The DDE solvers are all low order, and in the case of ddesd, it’s a low accuracy method which is fast for getting plots correct but not something which converges to many decimal places all too well since it doesn’t explicitly track discontinuities. This is even seen in the paper on the method which shows the convergence only matches dde23 to graphical accuracy on a constant-delay problem. Again, this fits in with the mantra of the suite, but may not hit all demographics. Shampine specifically made a separate version of ddesd for Fortran for people who are interested in efficiency, which is another way of noting that the key of ddesd is features and automatic usage, and not hardcore scientific computing efficiency. The mentioned SDE solver from the financial toolbox is only order 0.5 and thus requires quite a small dt to be accurate.

And I can keep going, but I think you get the moral of the story. This suite was created with one purpose in mind: to make it very easy to solve a wide array of differential equations and get a nice plot out. It does a very good job at doing so. But it wasn’t made with efficiency in mind, and so it’s missing a lot of methods that may be useful if you want high efficiency or high accuracy. Adding these wouldn’t make sense to the MATLAB suite anyways since it would clutter the offering. MATLAB is about ease of use, and not efficiency, and it does extremely well at what it set out to do. For software that was first made before the Y2K crisis with just a few methods added later, it still holds up very well.

Hairer’s Solvers (Fortran)

Next I want to bring up some Fortran solvers because they will come up later. Hairer’s Fortran solvers are a set methods with similar interfaces that were designed with efficiency in mind. Many of these methods are classics: dopri5, dop853, radau, and rodas will specifically show up in many of the suites which are discussed later. These methods are not too flexible: they don’t allow event handling (though with enough gusto you can use the dense output to write your own), or numbers that aren’t double-precision floating point numbers (it’s Fortran). They have a good set of options for tweaking parameters to make the adaptive timestepping more efficient, though you may have to read a few textbooks to know exactly what they do. And that’s the key to them: they will solve an ODE, stiff or non-stiff, and they will do so pretty efficiently, but nothing more. But even then, they show some age which don’t make them “perfectly efficient”. These solvers include their own linear algebra routines which don’t multithread like standard BLAS and LINPACK implementations, meaning that they won’t make full use of modern CPU architectures. The computations don’t necessarily SIMD or use FMA. But most of all, to use it directly you need to use Fortran which would be turn off for many people.

There is some flexibility in the offering. There are symplectic solvers for second order ODEs, the stiff solvers allow for solving DAEs in mass matrix form, there’s a constant-lag nonstiff delay differential equation solver (RETARD), there is a fantastic generalization of radau to stiff state-dependent delay differential equations (RADAR5), and there’s some solvers specifically for some “mechanical ODEs” commonly found in physical problems. Of course, to get this all working you’ll need to be pretty familiar with Fortran, but this is a good suite to look at if you are.

ODEPACK and Netlib ODE Solvers (Fortran)

ODEPACK is an old set of ODE solvers which accumulated over time. You can learn a bit about its history by reading this interview with Alan Hindenmarsh. I also bundle the Netlib suite together with it. There’s a wide variety of methods in there. There’s old Runge-Kutta methods due to Shampine, some of Shampine’s old multistep methods ddebdf and ddeabm, etc. The reason why this pack is really important though is because of the Lawarance Livermore set of methods, specifically LSODE and its related methods (LSODA, LSODR, VODE, etc.). It includes methods for implicit ODEs (DAEs) as well. These are a group of multistep methods which are descendent from GEAR, the original BDF code. They are pretty low level and thus allow you to control the solver step by step, and some of them have “rootfinding capabilities”. This means that you can use these to add event handling, though an event handling interface will take some coding. There’s lots of options and these are generally pretty performant for a large array of problems, but they do show their age in the same way that the Hairer codes do. Their linear algebra setups do not make use of modern BLAS and LINPACK, which in practical terms means they don’t make full use of modern computer architectures to speed things up. They don’t always make use of SIMD and other modern CPU acceleration features. They are the vanilla ODE solvers which have existed through time. In particular, LSODA is popular because this is the most widely distributed method which automatically detects stiffness and swtiches between integration methods, though it should be pointed out that there is a performance penalty from this feature.

Sundials and ARKCODE (C++ and Fortran)

Sundials‘ CVODE is a re-write of VODE to C which is a descendent of LSODE which is a descendent itself of the original GEAR multistep code. Yes, this has a long history. But you should think of it as “LSODE upgraded”: it makes use of modern BLAS/LINPACK, and also a bunch of other efficient C/Fortran linear solvers, to give a very efficient Adams and BDF method. Its solver IDA is like CVODE but handles implicit ODEs (DAEs). The interface for these is very similar to the ODEPACK interface, which means you can control it step by step and use the rootfinding capabilities to write your own event handling interface. Since the Adams methods handle nonstiff ODEs and the BDF methods handle stiff ODEs, this performance plus flexibility makes it the “one-stop-shop” for ODE solving. Many different scientific computing software internally are making use of Sundials because it can handle just about anything. Well, it does have many limitations. For one, this only works with standard C++ numbers and arrays. There’s no stochasticity or delays allowed. But more importantly, Hairer and DifferentialEquations.jl’s show that these multistep methods are usually not the most efficient methods since high order Rosenbrock, (E)SDIRK, and fully implicit RK methods are usually faster for the same accuracy for stiff problems, and high order Runge-Kutta are faster than the Adams methods for the same accuracy. Multistep methods are also not very stable at their higher orders, and so at higher tolerances (lower accuracy) these methods may fail to have their steps converge on standard test problems (see this note in the ROBER tests), meaning that you may have to increase the accuracy (and thus computational cost) due to stiffness issues. But, since multistep methods only require a single function evaluation per step (or are implicit in only single steps), they are the most efficient for asymptotically hard problems (i.e. when the derivative calculation is very expensive or the ODE is a large 10,000+ system). For this reason, these methods excel at solving large discretizations of PDEs. To top it off, there are parallel (MPI) versions for using CVODE/IDA for HPC applications.

I also want to note that recently Sundials added ARKCODE, a set of Runge-Kutta methods. These include explicit Runge-Kutta methods and SDIRK methods, including additive Runge-Kutta methods for IMEX methods (i.e. you can split out a portion to solve explicitly so that the implicit portion is cheaper when you know your problem is only partly or semi stiff). This means that it covers the methods which I mentioned earlier are more efficient in many of the cases (though it is a bit lacking on the explicit tableaus and thus could be more efficient, but that’s just details).

If you are using C++ or Fortran and want to write to only one interface, the Sundials suite is a great Swiss Army knife. And if you have an asymtopically large problem or very expensive function evaluations, this will be the most efficient as well. Plus the parallel versions are cool! You do have to live with the limitations of low-level software forcing specific number and array types, along with the fact that you need to write your own event handling, but if you’re “hardcore” and writing in a compiled language this suite is a good bet.

SciPy’s Solvers (Python)

Now we come to SciPy’s suite. If you look at what it offers, the names should now start to sound familiar: LSODA, VODE, dopri5, dop853. That is write: SciPy is simply a wrapper over Hairer’s and ODEPACK’s methods. Because it writes a generic interface though, it doesn’t have the granularity that is offered by ODEPACK, meaning that you don’t have step-by-step control and no event handling. The tweaking options are very basic as well. Basically, they let you define a function in Python, say what timeframe to solve it on, give a few tolerances, and it calls Fortran code and spits back a solution. It only has ODE solvers, no differential-algebraic, delay, or stochastic solvers. It only allows the basic number types and does no event handling. Super vanilla, but gets the job done? As with the methods analysis, it has the high order Runge-Kutta methods for efficient solving of nonstiff equations, but it’s missing Rosenbrock and SDIRK methods entirely, opting to only provide the multistep methods.

I should note here that it has the same limitation as MATLAB though, namely that the user’s function is Python code. Since the derivative function is where the ODE solver spends most of its time (for sufficiently difficult problems), this means that even though you are calling Fortran code, you will lose out on a lot of efficiency. Still, if efficiency isn’t a big deal and you don’t need bells and whistles, this suite will do the basics.

deSolve Package (R)

There’s not much to say other than that deSolve is very much like SciPy’s suite. It wraps almost the same solvers, has pretty much the same limitations, and has the same efficiency problem since in this case it’s calling the user-provided R function most of the time. One advantage is that it does have event handling. Less vanilla with a little bit more features, but generally the same as SciPy.

PyDSTool (Python)

PyDSTool is an odd little beast. Part of the software is for analytic continuation (i.e. bifurcation plotting). But another part of it is for ODE solvers. It contains one ODE solver which is written in Python itself and it recommends against actually using this for efficiency reasons. Instead, it wraps some of the Hairer methods, specifically dopri5 and radau, and recommends these. But it’s different than SciPy in that it takes in the specification of the ODE as a string, and compiles it to a C function, and uses this inside the ODE solver. By doing so, it’s much more efficient. We still note that its array of available methods is small and it offers radau which is great for high accuracy ODEs and DAEs, but is not the most efficient at lower accuracy so it would’ve been nice to see Rosenbrock and ESDIRK methods. It has some basic event handling and methods for DDEs (again as wrappers to a Hairer method). This is a pretty good suite if you can get it working, though I do caution that getting the extra (non-Python) methods setup and compiled is nontrivial. One big point to note is that I find the documentation spectacularly difficult to parse. Together, it’s pretty efficient and has a good set of methods which will do basic event handling and solve problems at double precision.

JiTCODE and JiTCSDE (Python)

JiTCODE is another Python library which makes things efficient by compiling the function that the user provides. It uses the SciPy integrators and does something similar to PyDSTool in order to get efficiency. I haven’t tried it out myself but I’ll assume this will get you as efficient as though you used it from Fortran. However, it is lacking in the features department, not offering advanced things like arbitrary number types, event handling, etc. But if you have a vanilla ODE to solve and you want to easily do it efficiently in Python, this is a good option to look at.

Additionally, JiTCDDE is a version for constant-lag DDEs similar to dde23. JiTCSDE is a version for stochastic differential equations. It uses the high order (strong order 1.5) adaptive Runge-Kutta method for diagonal noise SDEs developed by Rackauckas (that’s me) and Nie which has been demonstrated as much more efficient than the low order and fixed timestep methods found in the other offerings. It employs the same compilation setup as JitCODE so it should create efficient code as well. I haven’t used this myself but it would probably be a very efficient ODE/DDE/SDE solver if you want to use Python and don’t need events and other sugar.

Boost ODEINT Solver Library (C++)

The Boost ODEINT solver library has some efficient implementations of some basic explicit Runge-Kutta methods (including high order RK methods) and some basic Rosenbrock methods (including a high order one). Thus it can be pretty efficient for solving the most standard stiff and nonstiff ODEs. However, its implementations do not make use of advanced timestepping techniques (PI-controllers and Gustofsson accleration) which makes it require more steps to achieve the same accuracy as some of the more advanced software, making it not really any more efficient in practice. It doesn’t have event handling, but it is flexible with the number and array types you can put in there via C++ templates. It and DifferentialEquations.jl are the only two suites that are mentioned that allow for solving the differential equations on the GPU. Thus if you are familiar with templates and really want to make use of them, this might be the library to look at, otherwise you’re probably better off looking elsewhere like Sundials.

GSL ODE Solvers (C)

To say it lightly, the GSL ODE solvers are kind of a tragedy. Actually, they are just kind of weird. When comparing their choices against what is known to be efficient according to the ODE research and benchmarks, the methods that they choose to implement are pretty bizarre like extrapolation methods which have repeatedly been shown to not be very efficient, while not included other more important methods. But they do have some of the basics like multistep Adams and BDF methods. This, like Boost, doesn’t do all of the fancy PI-controlled adaptivity and everything though, so YMMV. This benchmark, while not measuring runtime and only uses function evaluations (which can be very very different according to more
sophisticated benchmarks like the Hairer and DifferentialEquations.jl ones!), clearly shows that the GSL solvers can take way too many function evaluations because of this and thus, since it’s using methods similar to LSODA/LSODE/dopri5, probably have much higher runtimes than they should.

Mathematica’s ODE Solvers

Mathematica’s ODE solvers are very sophisticated. It has a lot of explicit Runge-Kutta methods, including newer high order efficient methods due to Verner and Shampine’s efficient method mentioned in the MATLAB section. These methods can be almost 5x faster than the older high order explicit RK methods which themselves are the most efficient class of methods for many nonstiff ODEs, and thus these do quite well. It
includes interpolations to make their solutions continuous functions that plot nicely. Its native methods are able to make full use of Mathematica and its arbitrary precision, but sadly most of the methods it uses are wrappers for the classic solvers. Its stiff solvers mainly call into Sundials or LSODA. By using LSODA, it tends to be able to do automatic stiffness detection by default. It also wraps Sundials’ IDA for DAEs. It uses a method of steps implementation over its explicit Runge-Kutta methods for solving nonstiff DDEs efficiently, and includes high order Runge-Kutta methods for stochastic differential equations (though it doesn’t do adaptive timestepping in this case). One nice feature is that all solutions come with an interpolation to make them continuous. Additionally, it can use the symbolic form of the user-provided equation in order to create a function for the Jacobian, and use that (instead of a numerical differentiation fallback like all of the previous methods) in the stiff solvers for more accuracy and efficiency. It also has symplectic methods for solving second order ODEs, and its event handling is very expressive. It’s very impressive, but since it’s using a lot of wrapped methods you cannot always make use of Mathematica’s arbitrary precision inside of these numerical methods. Additionally, its interface doesn’t give you control over all of the fine details of the solvers that it wraps.

Maple’s ODE Solvers

Maple’s ODE solvers are pretty basic. It defaults to a 6th order explicit RK method due to Verner (not the newer more efficient ones though), and for stiff problems it uses a high order Rosenbrock method. It also wraps LSODE like everything else. It has some basic event handling, but not much more. As another symbolic programming language, it computes Jacobians analytically to pass to the stiff solvers like Mathematica, helping it out in that respect, but its offerings pale in comparison to Mathematica.

FATODE (Fortran)

FATODE
is a set of methods written in Fortran. It includes explicit Runge-Kutta methods, SDIRK methods, Rosenbrock methods and fully implicit RK methods. Thus it has something that’s pretty efficient for pretty much every case. What makes it special is that it includes the ability to do sensitivity analysis calcuations. It can’t do anything but double-precision numbers and doesn’t have event handling, but the sensitivity calculations makes it quite special. If you’re a FORTRAN programmer, this is worth a look, especially if you want to do sensitivity analysis.

DifferentialEquations.jl (Julia)

Okay, time for DifferentialEquations.jl. I left it for last because it is by far the most complex of the solver suites, and pulls ideas from many of them. While most of the other suite offer no more than about 15 methods on the high end (with most offering about 8 or less), DifferentialEquations.jl offers 200+ methods and is continually growing. Like the standard Python and R suites, it offers wrappers to Sundials, ODEPACK, and Hairer methods. However, since Julia code is always JIT compiled, its wrappers are more akin to PyDSTool or JiTCODE in terms of efficiency. Thus all of the standard methods mentioned before are available in this suite.

But then there are the native Julia methods. For ODEs, these include explicit Runge-Kutta methods, (E)SDIRK methods, and Rosenbrock methods. In each of these categories it has a huge amount of variety, offering pretty much every method from the other suites along with some unique methods. Some unique methods to point out are that it has the only 5th order Rosenbrock method, it has the efficient Verner methods discussed in the Mathematica part, it has newer 5th order methods (i.e. it includes the Bogacki-Shampine method discussed as an alternative to ode45’s tableau, along with an even newer tableau due to Tsitorious which is even more efficient). It has methods specialized to reduce interpolation error (the OwrenZen methods), and methods which are strong-stability presurving (for hyperbolic PDEs). It by default the solutions as continuous functions via a high order interpolation (though this can be turned off to make the saving more efficient). Each of these implementations employ a lot of extra tricks for efficiency. For example, the interpolation is “lazy”, meaning that if the method requires extra function evaluations for the continuous form, it will only do those extra calculations when the continuous function is used (so when you directly ask for it or when you plot). This is just a peak at the special things the library does to gain an edge.

The native Julia methods benchmark very well as well, and all of the benchmarks are openly available. Essentially, these methods make use of the native multithreading of modern BLAS/LINPACK, FMA, SIMD, and all of the extra little compiler goodies that allows code to be efficient, along with newer solver methods which theoretically reduce the amount of work that’s required to get the same error. They even allow you to tweak a lot of the internals and swap out the linear algebra routines to use parallel solvers like PETSc. The result is that these methods usually outperform the classic C/Fortran methods which are wrapped. Additionally, it has ways to symbolically calculate Jacobians like Mathematica/Maple, and instead of defaulting to numerical differentiation the stiff solvers fall back to automatic differentiation which is more efficient and has much increased accuracy.

Its event handling is the most advanced out of all of the offerings. You can change just about anything. You can make your ODE do things like change its size during the solving if you want, and you can make the event handling adjust and adapt internal solver parameters. It’s not a hyperbole to say that the user is given “full control” since the differential equation solvers themselves are written as essentially a method on the event handling interface, meaning that anything it can do internally you can do.

The variety of methods is also much more diverse than the other offerings. It has symplectic integrators like Harier’s suite, but has more high and low order methods. It has a range of Runge-Kutta Nystrom methods for efficiently solving second order ODEs. It has the same high order adaptive method for diagonal noise SDEs as JiTCSDE, but also includes high order adaptive methods specifically for additive noise SDEs. It also has methods for stiff SDEs in Ito and Stratanovich interpretations, and allows for event handling in the SDE case (with the full flexibility). It has DDE solvers for constant-lag and state-dependent delays, and it has stiff solvers for each of these cases. The stiff solvers also all allow for solving DAEs in mass matrix form (though fully implicit ODEs are possible, but can only be solved using a few methods like a wrapper to Sundials’ IDA and doesn’t include event handling here quite yet).

It allows arbitrary numbers and arrays like Boost. This means you can use arbitrary precision numbers in the native Julia methods, or you can use “array-like” objects to model multiscale biological organisms instead of always having to use simple contiguous arrays. It has addons for things like sensitivity analysis and parameter estimation. Also like Boost, it can solve equations on the
GPU by using a GPUArray.

It hits the strong points of each of the previously mentioned suites because it was designed from the get-go to do so. And it benchmarks extremely well. The only downside is that, because it is so feature and performance focused, its documentation is heavy. The beginning tutorial will give you the basics (making it hopefully as easy as MATLAB to pick up), but pretty much every control knob is available, making the extended portion of the documentation a long read.

Conclusion

Let’s take a step back and summarize this information. DifferentialEquations.jl objectively has the largest feature-set, swamping most of the others while wrapping all of the common solvers. Since it also features solvers for the non-ordinary differential equations and its unique Julia methods also benchmarks well, I think it’s clear that DifferentialEquations.jl is by far the best choice for “power-users” who are looking for a comprehensive suite.

As for the other choices from scripting languages, MATLAB wasn’t designed to have all of the most efficient methods, but it’ll handle basic equations with delays and events and output good plots. R’s deSolve is similar in most respects to MATLAB. SciPy’s offering is lacking in comparison to MATLAB and R’s due to the lack of event handling. But MATLAB/Python/R all have efficiency problems due to the fact that the user’s function is written in the scripting language. JiTCODE and PyDSTool are two Python offerings make the interface to the Fortran solvers more efficient than straight SciPy. Mathematica and Maple will do symbolic pre-calculations to speed things up and can JiT compile functions, along with offering pretty good event handling, and thus their wrappers are more like DifferentialEquations.jl in terms of flexibility and efficiency (and Mathematica had a few non-wrapper goodies mentioned as well). So in a pinch when not all of the bells and whistles are necessary, each of these scripting language suites will get you by. Behind DifferentialEquations.jl, I would definitely put Mathematica’s suite second for scripting languages with everything else much further behind.

If you’re already hardcore and writing in C++/Fortran, Sundials is a pretty good “one-stop-shop” to get everything you need, especially when you add in ARKCODE. Still, you’re going to have to write a lot of stuff yourself to make the rootfinding into an event handling interface, but if you put the work in
this will do you well. Hairer’s codes are a great set of classics which cover a wide variety of equations, and FATODE is the only non-DifferentialEquations.jl suite which offers a way to calculate sensitivity equations (and its sensitivity equations are more advanced). Any of these will do you well if you want to really get down-and-dirty in a compiled language and write a lot of the interfaces yourself, but they will be a sacrifice in productivity with no clear performance gain over the scripting language methods which also include some form of JIT compilation. With these in mind, I don’t really see a purpose for the GSL or Boost suites, and the ODEPACK methods are generally outdated.

I hope this review helps you make a choice which is right for you.

The post A Comparison Between Differential Equation Solver Suites In MATLAB, R, Julia, Python, C, Mathematica, Maple, and Fortran appeared first on Stochastic Lifestyle.

Julia iFEM3: Solving the Poisson Equation

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/julia-ifem3/

This is the third part in the series for building a finite element method solver in Julia. Last time we used our mesh generation tools to assemble the stiffness matrix. The details for what will be outlined here can be found in this document. I do not want to dwell too much on the actual code details since they are quite nicely spelled out there, so instead I will focus on the porting of the code. The full code will be available soon on a github repository, and so since most of it follows a straight translation from the linked document, I’ll leave it out of the post for you to find on the github.

The Ups, Downs, and Remedies to Math in Julia

At this point I have been coding in Julia for over a week and have been loving it. I come into each new function knowing that if I just change the array dereferencing from () to [] and wrap vec() calls around vectors being used as indexes (and maybe int()), I am about 95% done with porting a function. Then I usually play with cosmetic details. There are a few little details which make code in Julia a lot prettier. For example, for the FEM solver, we need to specify a function. In MATLAB, specifying such the function for which we want to solve -Delta u = f in-line would be done via anonymous functions like:

f = @(x)sin(2*pi.*x(:,1)).*cos(2*pi.*x(:,2));

I tend to have two problems with that code. For one, it’s not math, it’s programming, and so glancing at my equations and glancing at my code has a little bit of a translation step where errors can happen. Secondly, in many cases anonymous functions incur a huge performance decrease. This fact and the fact that MATLAB’s metaprogramming is restricted to string manipulation and eval destroyed a project I had tried a few years ago (general reaction-diffusion solver with a GUI, the GUI took in the reaction equations, but using anonymous functions killed the performance to where it was useless). However, in Julia functions can be defined inline and be first class, and have a nice appearance. For example, the same function in Julia is:

f(x) = sin(2π.*x[:,1]).*cos(2π.*x[:,2]);

Two major improvements here. For one, since the interpreter knows that variables cannot start with a number, it interprets 2π as the mathematical constant 2π. Secondly, yes, that’s a π! Julia uses allows unicode to be entered (In Juno you enter the latex pi and hit tab), and some of them are set to their appropriate mathematical constants. This is only cosmetic, but in the long-run I can see this as really beneficial for checking the implementation of equations.

However, not everything is rosy in Julia land. For one, many packages, including the FEM package I am working with, are in MATLAB. Luckily, the interfacing via MATLAB.jl tends to work really well. In the first post I showed how to do simple function calls. This did not work for what I needed to do since these function calls don’t know how to pass a function handle. However, digging into MATLAB.jl’s package I found out that I could do the following:

  put_variable(get_default_msession(),:node,node)
  put_variable(get_default_msession(),:elem,elem)
  put_variable(get_default_msession(),:u,u)
  eval_string(get_default_msession(),"sol = @(x)sin(2*pi.*x(:,1)).*cos(2*pi.*x(:,2))/(8*pi*pi);")
  eval_string(get_default_msession(),"Du = @(x)([cos(2*pi.*x(:,1)).*cos(2*pi.*x(:,2))./(4*pi) -sin(2*pi.*x(:,1)).*sin(2*pi.*x(:,2))./(4*pi)]);")
 
  eval_string(get_default_msession(),"h1 = getH1error(node,elem,Du,u);");
  eval_string(get_default_msession(),"l2 = getL2error(node,elem,sol,u);");
  h1 = jscalar(get_mvariable(get_default_msession(),:h1));
  l2 = jscalar(get_mvariable(get_default_msession(),:l2));

Here I send the variables node, elem, and u to MATLAB. Then I directly evaluate strings in MATLAB to make function handles. With all of the variables appropriately in MATLAB, I call the function getH1error and save its value (in MATLAB). I then use the get_mvariable to bring the result into MATLAB. That value is a value of MATLAB type, and so I use the MATLAB.jl’s conversion function jscalar to then get the scalar result h1. As you can see, using MATLAB.jl in this fashion is general enough to do any of your linking needs.

For very specialized packages, this is good. For testing the ported code for correctness, this is great. However, I hope not to do this in general. Sadly, every once in awhile I run into a missing function. In this example, I needed accumarray. It seems I am not the only MATLAB exile as once again Julia implementations are readily available. The lead Julia developer Stefan has a general answer:

function accumarray2(subs, val, fun=sum, fillval=0; sz=maximum(subs,1), issparse=false)
   counts = Dict()
   for i = 1:size(subs,1)
        counts[subs[i,:]]=[get(counts,subs[i,:],[]);val[i...]]
   end 
   A = fillval*ones(sz...) 
   for j = keys(counts)
        A[j...] = fun(counts[j])
   end
   issparse ? sparse(A) : A
end
  0.496260 seconds (2.94 M allocations: 123.156 MB, 8.01% gc time)
  0.536521 seconds (2.94 M allocations: 123.156 MB, 8.83% gc time)
  0.527007 seconds (2.94 M allocations: 123.156 MB, 9.41% gc time)
  0.544096 seconds (2.94 M allocations: 123.156 MB, 9.76% gc time)
  0.526110 seconds (2.94 M allocations: 123.156 MB, 12.22% gc time)

whereas Tim answer has less options but achieves better performance:

function accumarray(subs, val, sz=(maximum(subs),)) 
    A = zeros(eltype(val), sz...) 
    for i = 1:length(val) 
        @inbounds A[subs[i]] += val[i] 
    end 
    A 
end

Timings

  0.000355 seconds (10 allocations: 548.813 KB)
  0.000256 seconds (10 allocations: 548.813 KB)
  0.000556 seconds (10 allocations: 548.813 KB)
  0.000529 seconds (10 allocations: 548.813 KB)
  0.000536 seconds (10 allocations: 548.813 KB)
  0.000379 seconds (10 allocations: 548.813 KB)

Why is Julia “Missing” So Many Functions? And How Do We Fix It?

This is not the first MATLAB function I found myself missing. Off the top of my head I know I had to get/make versions of meshgrid, dot(var,dimension) just in the last week. To the dismay of many MATLAB exiles, many of the Julia developers are against “cluttering the base” with these types of functions. While it is easy to implement these routines yourself, many of these routines are simple and repeated by mathematical programmers around the world. By setting a standard name and implementation to the function, it helps code-reusability and interpretability.

However, the developers do make a good point that there is no reason for these functions to be in the Base. Julia’s Base is for functions that are required for general use and should be kept small in order to make it easier for the developers to focus on the core functionality and limit the resources required for a standard Julia install. This will increase the number of places where Julia could be used/adopted, and will help ensure the namespace isn’t too full (i.e. you’re not stepping on too many pre-made functions).

But mathematicians need these functions. That is why I will be starting a Julia Extended Mathematical Package. This package will be to hold the functions that are not essential language functions, but “essential math language” functions like you’d find in the base of MATLAB, R, numpy/scipy, etc., or even just really useful routines for mathematical programming. I plan on cleaning up the MATLAB implementations I have found/made ASAP and putting this package up on github for others to contribute to. My hope is to have a pretty strong package that contains the helper functions you’d expect to have in a mathematical programming language. Then just by typing using ExtendedMath, you will have access to all the special mathematical functionality you’re used to.

Conclusion

As of now I have a working FEM solver in Julia for Poisson’s equation with mixed Neumann and Dirichlet boundary conditions. This code has been tested for convergence and accuracy and is successful. However, this code has some calls to MATLAB. I hope to clean this up and after this portion is autonomous, I will open up the repository. The next stage will be to add more solvers: more equations, adaptive solvers, etc., as I go through the course. As, as mentioned before, I will be refactoring out the standard mathematical routines and putting that to a different library which I hope to get running ASAP for others to start contributing to. Stay tuned!

The post Julia iFEM3: Solving the Poisson Equation appeared first on Stochastic Lifestyle.

Julia iFEM 2: Optimizing Stiffness Matrix Assembly

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/julia-ifem2/

This is the second post looking at building a finite element method solver in Julia. The first post was about mesh generation and language bindings. In this post we are going to focus on performance. We start with the command from the previous post:

node,elem = squaremesh([0 1 0 1],.01)

which generates an array elem where each row holds the reference indices to the 3 points which form a triangle (element). The actual locations of these points are in the array node, and so node(1) gives the points in the (x,y)-plane for the ith point. What the call is saying is that these are generated for the unit square with mesh-size .01, meaning we have 10201 triangles.

The approach to building the stiffness matrix for the Poisson equation is described here. The general idea is that span our vector space by a basis of hat functions phi_{i}, and the so the stiffness matrix is found by the inner product (integral) between these basis functions. This translates to solving for the area of the triangles where two hat functions overlap, which we can do exactly since we chose the basis to be sufficiently nice. However, since the hat functions are zero except in a small range, most of these inner products are zero, meaning the resulting stiffness matrix is sparse. Our goal is to produce this sparse matrix as efficiently as possible. Lets get to it!

Building the Julia Version: Local Stiffness

The first function we need is titled localStiffness, which evaluates the inner product to give the “stiffness for one triangle”. In MATLAB the code is of the form:

function [At,area] = localstiffness(p)
    At = zeros(3,3);
    B = [p(1,:)-p(3,:); p(2,:)-p(3,:)];
    G = [[1,0]',[0,1]',[-1,-1]'];
    area = 0.5*abs(det(B));
    for i = 1:3
        for j = 1:3
            At(i,j) = area*((BG(:,i))'*(BG(:,j)));
        end
    end
end

It takes in a vector p of three points, solves for the area at the reference triangle, and transforms the area appropriately to give the stiffness for the triangle defined by the points of p. I want to make a few points about this code. First of all, it employs a “trick” for solving for the dot product. That is, it uses the transposed vector times another vector. I quote trick because to a mathematician, this is simply the definition of the dot product, and so it only seems natural to use it like this in MATLAB. However, two things to point out. First of all, in Julia, such an operation does not return a scalar, but a one dimensional vector. This in some cases will give unexpected errors, so it should be avoided. Not only that, but it is also inefficient. In both MATLAB and Julia, matrices are stored column-wise, that is they are stored as a array of pointers where the pointers go to an array of columns. Thus to access a row matrix, both MATLAB and Julia would have to access the pointer and then go to the array at which it points (a size 1 array), and take the value there. Notice this is an extra step. Therefore it’s more efficient to keep the vectors column-wise. (In reality, the vectors here are so small that it won’t make a difference, but this justifies that the change we will do for the Julia code in a performance-sense).

Two other small changes were made to this code. For one, Julia throws an error at the definition of G since it cannot read the transpose calls inside of an array declaration. This is easily fixed by simply transposing after the creation instead. Lastly, we change the pre-allocation of At from zeros to Julia’s general constructor. This is slightly more performant since it will allocate the space without doing an initial re-write step, saving the time it would take to loop through and set each value to zero. The result is the following:

function localstiffness(p)
  At = Array{Float64}(3,3);
  B = [p[1,:]-p[3,:]; p[2,:]-p[3,:]];
  G = [1 0 ; 0 1 ; -1 -1]';
  area = 0.5*abs(det(B));
  for i = 1:3
    for j = 1:3
      At[i,j] = area*dot(BG[:,i],BG[:,j]); 
    end
  end
  return(At,area)
end

which is ever so slightly more efficient, but in reality the same and just tweaked for the quirks.

Building the Julia Version: Matrix Assembly

Now we need to loop through each triangle and sum up the inner product between each pair of basis functions over each triangle. The intuitive code is:

function assemblingstandard(node,elem)
  N=size(node,1); NT=size(elem,1);
  A=zeros(N,N);
  for t=1:NT
    At,=localstiffness(node[vec(elem[t,:]),:]);
    for i=1:3
      @simd for j=1:3
        @inbounds A[elem[t,i],elem[t,j]]=A[elem[t,i],elem[t,j]]+At[i,j];
      end
    end
  end
  return(A)
end

I discussed previously the use of macros to speed up code without changing its style. The main problem that had to addressed in porting the code to Julia was that Julia will only take a vector for array referencing when the colon operator is used. Therefore since elem[t,:] returns a 1×3 Matrix of indices for the points associated with triangle t, once again a 1×3 Matrix is not an array in Julia so it throws an error. This is easy to fix by wrapping it in the vec() command, which others have tested to be the fastest method for conversion, and will actually do some fanciness in the background in order to not have to copy the array. This means that the cost of using vec is so small that I will use it liberally as a fix in these cases. Notice that within the loop vec is not required to reference A. This is because the issue only occurs when the colon operator is present.

However, since At is usually zero, we can improve this a lot by instead generating vectors to build a sparse matrix. What we will instead do is save the (i,j) pairs where the value should be stored, and the value, and use the sparse command to reduce. The sparse command will automatically add together the values from repeated (i,j) indices, effectively performing the update we had before. This gives us the code:

function assemblingsparse(node,elem)
  N = size(node,1); NT = size(elem,1);
  i = Array{Int64}(NT,3,3); j = Array{Int64}(NT,3,3); s = zeros(NT,3,3);
  index = 0;
  for t = 1:NT
    At, = localstiffness(node[vec(elem[t,:]),:]);
    for ti = 1:3, tj = 1:3
        i[t,ti,tj] = elem[t,ti];
        j[t,ti,tj] = elem[t,tj];
        s[t,ti,tj] = At[ti,tj];
    end
  end
  return(sparse(vec(i), vec(j), vec(s)));
end

Notice here that vec() needed to be used to build the sparse matrix because the easiest way to index within the loop was to use a 3-dimensional array and then flatten it via vec(). Notice a neat Julia trick where you can define multiple for loops in one line: “for ti = 1:3, tj = 1:3” is two loops and works as you’d expect. With some extra work we can get rid of the loop over the triangles by performing the calculations from localstiffness vector-wise, which gives us the vectorized form:

function assembling(node,elem)
  N = size(node,1); NT = size(elem,1);
  ii = Vector{Int64}(9*NT); jj = Vector{Int64}(9*NT); sA = Vector{Float64}(9*NT);
  ve = Array{Float64}(NT,2,3)
  ve[:,:,3] = node[vec(elem[:,2]),:]-node[vec(elem[:,1]),:];
  ve[:,:,1] = node[vec(elem[:,3]),:]-node[vec(elem[:,2]),:];
  ve[:,:,2] = node[vec(elem[:,1]),:]-node[vec(elem[:,3]),:];
  area = 0.5*abs(-ve[:,1,3].*ve[:,2,2]+ve[:,2,3].*ve[:,1,2]);
  index = 0;
  for i = 1:3, j = 1:3
     @inbounds begin
     ii[index+1:index+NT] = elem[:,i];
     jj[index+1:index+NT] = elem[:,j];
     sA[index+1:index+NT] = sum(ve[:,:,i].* ve[:,:,j],2) ./(4*area); # Replacing dot(ve[:,:,i],ve[:,:,j],2)
     index = index + NT;
     end
   end
  return(sparse(ii,jj,sA));
end

Here I note that Julia’s dot product will not act on matrices, only vectors. In order to do the row-wise dot product like we would do in MATLAB, we can simply use .* and sum the results in each row.

Dispelling a Myth: Vectorized Julia Rocks!

The most common complaint about Julia that people tend to have is that, in many cases, the code which gets the most performance is the de-vectorized code. “But the vectorized code can be so beautiful! Why would I want to change that?”. This myth seems to come from the alpha days or really bad tests, but it doesn’t seem to die. Instead, what I wish to show here is that vectorized code is also faster in Julia. First we run some basic timings:

@time assemblingstandard(node,elem);
@time assemblingsparse(node,elem);
@time assembling(node,elem);
 
0.538606 seconds (2.52 M allocations: 924.214 MB, 12.79% gc time)
0.322044 seconds (2.52 M allocations: 137.714 MB, 21.53% gc time)
0.015182 seconds (775 allocations: 28.650 MB, 16.97% gc time)

These timings pretty stability show this pattern. All of the methods were tried with parallelization and simd options with either no speedup or it being detrimental given the problem size. What this shows is that, out of the intuitive forms for solving these equations, the vectorized form was by far the fastest. The reason is that this is a highly vectorizable problem, whereas I discussed before the limitations that can cause vectorization to lose to devectorization.

But how does this fare against MATLAB? The “same” code was run in MATLAB (this is from iFEM, an optimized library Professor Long Chen) which gives the results:

tic; assemblingstandard(node,elem); toc;
tic; assemblingsparse(node,elem); toc;
tic; assembling(node,elem); toc;
 
Elapsed time is 0.874919 seconds.
Elapsed time is 0.698763 seconds.
Elapsed time is 0.020453 seconds.

To ensure the time difference between the vectorized versions, we had both problems solve it in a loop:

@time for i = 1:1000
  assembling(node,elem);
end
 
 9.312876 seconds (821.25 k allocations: 27.980 GB, 21.32% gc time)

vs MATLAB:

   tic; 
   for i=1:1000 assembling(node,elem); 
   end
   toc;
 
Elapsed time is 19.221982 seconds.

Thus, in line with what this coder found with R, the vectorized code ran more than twice as fast in Julia.

The Julia Optimization Mentality

This was a fun little exercise because I had no idea how it would turn out. Quite frankly, when I started I thought MATLAB would slightly edge out Julia when running such vectorized code. However, Julia continues to impress me. The only major problem that I had this time around was finding out to use the vec() function to deal with “1-dimensional matrices”, but once that was found I was able to get Julia to be faster than MATLAB, even though I know much more about MATLAB and this package itself is quite well optimized.

I think I should end on why I find the Julia philosophy compelling for scientific computing. The idea is not that “you have to devectorize to get the best code”, though there are situations where doing so can dramatically increase your speed. The idea is that you don’t have to contort yourself to vectorization to make everything work. In MATLAB, R, and Python, you have to vectorize in order to make your code to ever finish. That is not the case in Julia. Here, just write the code that seems natural and it will do really well. In this case, vectorized code was natural, and as you could see we got something that was even faster. To do better in MATLAB, we would at this stage have to start writing C/MEX code. In Julia, we could expand out the loops, play with adding SIMD calls, caching, etc. directly in the Julia language.

For scientific computing where we just want code that’s good enough to work, you can see it’s easier to get there with Julia. If you need to optimize it to be part of a library, you can optimize and devectorize it within Julia without having to go to C (many times by just adding macros throughout your code). Will it be as fast as C? No, but many tests show that you’ll at least get within a factor of two so, for almost every case, you might as well just code it in Julia and move on. Each blog post I do I am getting more and more converted!

The post Julia iFEM 2: Optimizing Stiffness Matrix Assembly appeared first on Stochastic Lifestyle.