6 Months of DifferentialEquations.jl: Where We Are and Where We Are Going



Hello everyone,

I invite you to read the latest summary on what’s going on in DifferentialEquations.jl and the JuliaDiffEq universe. In this blog post I discuss the overarching design, what it has allowed us to achieve, and where we are going next. I hope this both inspires other developers, being more contributors to the ecosystem, and helps grow interest in Julia itself.

Lastly, if you find these developments interesting and would like to support our development efforts, please star the DifferentialEquations.jl repository. I hope to use these types of metrics to show adoption of the program in order to one day be funded and make this a full time thing (I love Julia!). Thank you for your support (especially those who are contributors to JuliaDiffEq!).


Wow, that looks really exciting, thanks for the writeup!


Thanks for the nice writeup. I especially like the idea of using callable model types instead of closures, very elegant and will be useful for me in other contexts.


This is starting to look like a fantastic effort, I’ll start exploring it next time I have the opportunity. And can I also say how grateful I am that you made DifferentialEquations.jl a super-package, given how generic the name is.


Hey Chris, just wanted to say that I had been reading the tutorials and was really impressed by the high quality of the work.


Hi Chris,

Sorry for the late reply, very impressive writeup!

Regarding the solution object that is returned: how would that work in case of e.g. a finite element solution, where you may have a solution over some 2D or 3D unstructured grid at every timestep? Is the design able to account for that?

@BenLauwens this was a very interesting read, do you think it’s possible to fit SimJulia.jl into the same framework, allowing it to be used as a method in solve?


Hi Chris, Bart

Indeed very inspiring!
The quantized state system solver that I am implementing in SimJulia looks at the same equations but from a very different perspective, e.g. the ParametrizedFunctions package is brilliant but I need other inputs and outputs than the ones that are provided.
What I intent to do is to mimic the philosophie of the packages in Julia DiffEq as a high level interface to SimJulia in a separate package, i.e. provide a solve function. This means that the quantized state system solver can be similarly used as the other ODE solvers. The coding and performance issues with event handling are however completely but differently solved within a quantized state system framework and porting the callback macros is not trivial.
I hope to finish the implementation for both non-stiff and stiff ODEs in SimJulia in January and to have the high level interface ready by February.


Hey everyone,
Thank you for all of the kind words!

You can check out what my first draft of the end game (finite elements for unstructured grids to solve stochastic partial differential equations) looks like at this tutorial page:


It’s completely usable and you can already do some cool stuff, but let me say up front that this part is only a first draft and not well-optimized for finite element PDE solvers (I mentioned that in the blog post as well). In fact, I wrote this before most of the ODE solvers, and so it served as a jumping off point for how to “do this better for a 2.0”, namely:

  1. It needs to be changed to inplace updates instead of the vectorization I used.
  2. It needs to utilize ODE solvers. I actually wrote each method directly in here. The reason was because not only could no ODE solver do splitting, but no ODE solver could use the fact that the Laplacian is a linear operator, and thus when split we need to solve a linear problem. The naive splitting for Crank-Nicholson (Trapezoid rule) in ODEs actually requires a nonlinear solve.
  3. It needs to have a way for users to choose how to solve the nonlinear and linear problems. I had to implement a bunch of options to make this pretty extensive.
  4. Choosing algorithms by symbols results in dynamic dispatch. It’s not a big deal since most of the time is spent in an inner function and so the dispatch cost is fairly minimal, but it’s still there.

If you notice, all four of these problems were addressed by the design of the ODE/SDE/DAE solvers, and so when the setup with all of the special problem types and promotion is complete, the finite element solvers will switch over to just discretizing the problem, creating the matrices for the operators, and calling the ODE solvers (giving users the option for which algorithm). Since they will be able to recognize linearity, then for example Trapezoid() would directly give the Crank-Nicholson method with linear solves you know and love for the Heat Equation and all of that jazz. In fact, the event framework for ODEs/SDEs/DAEs is general enough that you can use it to implement adaptive meshing as well: use whatever criteria you want to trigger a meshing change, and the events can change the size of the ODE and replace the linear operator. This case will need a slightly different solution object to save a timeseries of meshes as well, but that’s not difficult.

But this is why the API for finite element methods is not quite up to speed with the other solvers: it’s old and there’s no reason to update it until this big change occurs.

All of that aside, a finite element problem discretizes into an ODE on matrices, where each column vector is the spatial discretization of a variable. So sol[i] is the matrix of values at the ith timestep, sol[i,j,:] are all of the variables at the jth location in space, and you have sol.mesh.elem[j] which gives [n1 n2 n3] as the nodes for the triangles, and sol.mesh.node[n] gives the (x,y) (and in the future z) coordinates for the node n.

So all of the information is there, but there’s probably a little bit of work to be done to make it a little bit more accessible. One thing which goes a long way here is plot recipes: the plot recipes take this and automatically make trisurf plots and animations. Note that there is a problem with Plots.jl that its trisurf doesn’t let you choose what node triples create triangles (like PyPlot does with its triangles argument), so the plots don’t look so great right now, but that issue was already noted a long time ago and fixes will come.

But together, this package shows a way to do finite element discretizations, and I’ll use it to provide some linear finite element discretizations. That should also serve as a good guideline for how other finite element discretizations could be done, which would have their own problem/solution types acting in a similar manner.


I see. I took a look at SimJulia and read up on quantized state systems and indeed it sounds like there’s a big gap in how they work. However, if your algorithms could take in an ODEProblem and spit out a solution type which has the same actions, it could integrate into the common interface quite well. It would be very interesting to benchmark and see where QSSs outperform classical methods, and have a quick switch to swap over to them.

Note that this is one reason among many why there is a new version of events/callbacks being implemented soon. Events are going to be split from callbacks since that higher-level API is much easier to map to different solvers, while callbacks will still exist but should be much more niche. Take a look at the proposal for the new event API:

Let me know if you would need any change to the proposal to make it work for your case.


Wow that’s ALOT of effort you put in there! I am totally shocked that this module had only started 6 months ago back then.

I just have one question for you, is there any plan to fully support finite element analysis? From a google search, I found the following modules which seem like silo efforts made by different people, but I think they are all still far from mature.

  1. JuAFEM.jl
  2. JuliaFEM.jl
  3. EllipticFEM.jl
  4. FinElt.jl

I have also watched @kristoffer.carlsson’s talk in 2016 about JuAFEM.jl https://www.youtube.com/watch?v=30TUEhbGmuc but I wonder if there have been/will be any further development in this package to turn it into another super-packge perhaps.

Please let me know any thoughts you have on this (everyone), because I am thinking of contributing seriously some time soon. I will soon be starting my PhD and I will definitely need an FEA package with a lot of control, customizability and performance, and Julia seems like the perfect language for that. Seeing how you demonstrate Julia’s awesomeness into DifferentialEquations.jl, I am even more convinced and motivated to start something along the same lines but for more specific (or perhaps even generic) classes of PDEs. I especially like the ability to setup parametric functions that can be symbollically differentiated and used seemlessly with the rest of the code at almost no cost of elegance!

If you are wondering what my research interest is/will be, it is a blend of optimization theory and mechanical design, so at some point I would like to mix up the PDEs with other constraints and an objective in a so-called variational program, and hopefully develop some scalable algorithms to optimize those models, hence the customizability I need, and possible integration with JuMP. That is of course unless my supervisor thinks I have to work on something slightly different! But these are the thick lines anyways.

ETA: It seems that I missed the date this post was published and a few other details on FEA in the replies to previous comments, so please excuse my sleepiness :slight_smile:

ETA2: I notice there is documentation for solving the Heat equation and Poisson’s equation using FEM which certainly shows the intention, but I would also like to ask you if you think it is fundamentally possible to interface JuMP.jl and DifferentialEquations.jl for solving variational programs, basically discretizing the differential equations and integrals into a set of JuMP expressions rather than just directly solving them using the machinery of DifferentialEquations.jl.


Diggin’ up some old posts :slight_smile:.

To an extent, yes. Let me be clear about one thing. I specialize in timestepping methods. ODEs/SDEs/DAEs/DDEs, adaptivity, dense output, etc. What I know well is dealing with solving discretizations in time. This means that, while I do solve (S)PDEs in my own applications work, I am not an expert in spatial discretizations. Nor is it something I can devote all of my time to.

Building an FEM toolbox is something that can take an entire career. I have a very good niche for research and application in timestepping methods, and am making enough progress that I do not plan on shifting any time soon. That means I cannot devote the time to fully develop and FEM toolbox.

However, what I can do is provide a nice interface over FEM, spectral decomposition, and FDM (well, we’re making these since no spatial expert finds these interesting enough :slight_smile:) toolboxes that make solving common PDEs simple. So that’s what we’re doing. Right now, there are not sufficient tools in Julia, so as part of a GSoC project we are writing a wrapper to FEniCS. We are building automatic-generation of FDM operators (GSoC project), and interfacing with spectral decomposition tools like ApproxFun.jl. With all of that integrated, we’re putting the DiffEq interface on common PDEs so users just define a some functions, initial conditions, (a mesh), and bam get a result without needing to know the underlying tools. That’s what we’ll offer for PDEs: the timestepping tools required to solve PDEs with time, and an easy interface over more complex spatial discretization tools.

Hopefully some group will get a major FEM toolbox in Julia, and we’ll then use our tools as a front-end to that, and make sure our timestepping tools are compatible with it in every way possible.

Variational inference is a very broad term. Are you talking about optimal control? JuMP is a bit hard to interface with here because of its macro interface. Going directly to MPB is a bit easier. This is something that’s still being worked on. It should fundamentally be possible since all of the tools can be exposed, but it’ll take some work.

If you’re talking about having DiffEqs in the objective function, we already have that down with parameter estimation.



Thanks for your answer.

I see.

That will be quite helpful.

Cool this sounds awesome, but I guess for my own purposes FEniCS may not be what I am looking for, but for most FEM applications out there it is, so I am glad you are including wrapping it on your agenda.

I will be glad to contribute if such a group shows up, or perhaps JuAFEM.jl may be the center of such a group, seems like the most actively developed one of the lot. I beleive I should give the src a thorough read.

Variational inference and variational programming are 2 different things, the latter is a form of mathematical programming where some constraints are functionals over global time or space domains and other constraints are differential equations to hold for all time and/or space points. The techniques used to solve them revolve around formulating the KKT conditions of the discretized functionals and DEs and attempting to solve those. The discretization/derivation is done using some calculus of variations, and it is a form of mathematical programming hence the term variational programming. I am not a deep reader of optimal control, but I beleive they are heavily overlapping, perhaps optimal control leans more towards dynamic versions of some classes of variational programs.

Close but no, it’s not parameter estimation I am trying to do, I should probably read the src though, thanks for sharing. Ideally I am looking for something that spits out a mathematical program so I can apply some optimization algorithms on it, the constraints and/or objective could be functionals. Making the discretization by hand and writing the mathematical program using JuMP is another alternative, but I prefer a more generic approach that makes it more usable in the future. I will definitely be writing the optimization algorithms in Julia, and hopefully interfacing them with MPB and/or JuMP, so I was thinking the modelling should also converge to MBP or JuMP to hook all my work to Julia’s optimization ecosystem.

Well I suppose there is a lot of work I need to do then.

I will probably keep this forum updated to seek some advice and avoid reinventing the wheel, so I would appreciate your advice now and then :slight_smile:


Please add my new package to the list, https://github.com/PetrKryslUCSD/FinEtools.jl.




According to this comparison, everything else bites the dust compared to DifferentialEquations.jl!


Yep, very impressive! Nice work Chris!

(Although in that graph, as mentioned in the blog-text, maybe there should also be a line “Implicitly defined ODE/DAE” with a yellow box.)


Yes indeed this is a spot where we are weak and I’ve actually put this at the top of the list for “things for DiffEq 4.0” (the DiffEq 3.0 blogpost is coming in about a week). I’ll add a row for this and update the chart right now. Thanks for pointing it out. Anything else you noticed?

The purpose of this exercise was to look through every single other popular offering and find out what we don’t have in order to add it to the list of next developments. But it’s of course showing my implicit biases, so I put this out there so that people would correct my eye and lead me to cool new features we don’t have. So please point out everything you see that fell below the radar.

(And if anyone wants to work on custom allocators / the integrator interface for event handling with Sundials.jl, that would be the quickest way to bring this to “Good”. Just throwing that out there :slight_smile:. )


You could also include the PETSc solvers. Other than that all good, nice research. And I’m enjoying reading the interview with Hindmarsh you linked to (link).