DifferentialEquations.jl 4.5: ABC, Adaptive Multistep, Maximum A Posteriori
Once again we stayed true to form and didn’t solve the problems in the development list but adding a ton of new features anyways. Now that Google Summer of Code (GSoC) is in full force, a lot of these updates are due to our very awesome and productive students. Here’s what we got.
Approximate Bayesian Computation (ABC)
Marc Williams (@marcjwilliams1) contributed the abc_inference
function to
DiffEqBayes.jl which utilizes approximate Bayesian computation (also known as
ABC) from ApproxBayes.jl
to perform the estimation of parameter posteriors. Compared to standard
Bayesian methods, ABC is computationally cheaper and faster at the cost of
making some approximations (the name is quite appropriate). For models with lots
of parameters where you have a good prior guess for the posterior parameter
point estimates, ABC can be a great way to get posterior distributions.
Maximum A Posteriori Estimation
GSoC student Vaibhav Dixit (@Vaibhavdixit02) added maximum a posteriori estimation to our existing parameter estimation routines in DiffEqParamEstim.jl. Now these optimizationbased methods can take into account prior distributions which can help global optimizers stay in the right parameter ranges and improve fitting. All it takes is passing an optional prior distribution to any of the existing methods, so it’s very easy to add it your current work!
Multiple shooting objective
In other parameter estimation news, Vaibhav Dixit (@Vaibhavdixit02) added the
multiple_shooting_objective
which, like the build_loss_objective
method,
allows for fitting ODE results to any loss function. However, multiple shooting
methods are naturally more robust by solving simultaneously from many different
time points.
DynamicSS Solver
For SteadyStateProblem
types, a new solver DynamicSS
has been added that
utilizes the ODE solvers to find steady states. It builds in a callback that
will halt when the derivative is sufficiently small, allowing it to be a
robust automated steadystate finding machine.
Unified retcode handling
Takafumi Arakaki (@tkf) performed an awesome refactoring in the DiffEq core which allows all of the DiffEq solvers to utilize the same retcode handling code. This means that we have unified exit warning messages and the resulting retcodes. Each package should act the same for the common retcodes.
Fixed order variable time step Adams methods
Shubham Maddhashiya (@sipah00) contributed a variety of new adaptive AdamsBashforth and AdamsBashforthMoulton methods. These methods are fixed order 35 which minimize function evaluations and are designed for large nonstiff ODE discretizations. They utilize RungeKutta methods to hotstart the Adams methods, making them more efficient than variable order versions when only lower orders are required or when there are a lot of events.
ABDF2
Yingbo Ma (@YingboMa) contributed ABDF2
, an adaptive implementation of the
Backwards Differentiation Formula (BDF) order 2 method. BDF schemes are the
goto choice for PDE discretizations since they only have one step (one
nonlinear equation to solve) for each time step, minimizing the number of
function calculations. Unlike adaptive order BDF discretizations, this adaptive
BDF2 is ABLstable, meaning that it should be stable for any stiff ODE that is
thrown at it. This makes it a nice testing or anchor method: one that can always
be relied on to do wellenough.
In development
A lot of the next developments will come from our GSoC students. Here’s a list of things we are aiming for:

SABDF2
, which is a strong order 0.5 adaptive BDF2 implementation for stochastic differential equations which is order 2 for small noise SDEs. This will be the first implicit adaptive integrator for small noise SDEs and will be a great choice for SPDEs. 
Adaptive order AdamsBashforthMoulton. This will be JuliaDiffEq’s native Julia implementation of an AOAT Adams method. It will utilize variable coefficient formulas like Shampine’s
ddeabm
which has been shown in sources like Hairer’s Solving Ordinary Differential Equations to be more efficient than the more standardCVODE_Adams
since it’s more easily able to utilize higher order steps. 
Adaptive order Nordsieck methods. Nordsieck methods are the special implementation multistep methods which Sundials’ CVODE uses. This implementation utilizes a fixed leading coefficient which makes it more optimized for solving large PDE discretizations since it can reuse the Jacobian between steps while adapting time steps. Other problem solving environments such as MATLAB and SciPy utilize quasifixed timesteps to get a similar benefit, but this has the added cost of fixing time steps for many steps which reduces stepping efficiency. The reason this is done is simple: Nordsieck implementations are hard which is why the only existing ones (that I know of) are EPISODE and VODE/CVODE (and VODE is an adaption of EPISODE, so it’s really one code!). However, our GSoC student has been hammering away at this for months and we’ve made great progress, with the fixed order Adams Nordsieck method having already merged. The last step is to figure out order adaptivity and then the translation to BDF coefficients is trivial. While this methodology isn’t as optimized for nonstiff ODEs in its Adams form, this BDF form will likely be the new goto method in DiffEq for large stiff PDE discretizations.

Until exponential integrators! Xingjian Guo (@MSeeker1340) has been doing an extensive study of the numerical implementations of
expmv
andphimv
for building efficient exponential integrators. If you’re curious, you can read up on some of the discussions and here. He is quite close to having both efficient and numerically stable methods for lazy adaptive Krylovexpmv
andphimv
which will allow for fast c alculations of exponential integrators without computing full matrix exponentials. This will then be used in the current exponential integrator implementations, along with the adaptive high order EPIRK methods. This class of integrators will be one to keep aware of if you’re interested in timedependent PDEs. It will be great to compare and contrast between these methods and BDF integrators in these problems. 
Yiannis Simillides (@ysimillides) keeps making improvements to FEniCS.jl. At this part a large portion (a majority?) of the tutorial works from Julia.

Vaibhav Dixit (@Vaibhavdixit02) finished most of the parameter estimation methods on our wish list, so he’s onto implementing global sensitivity analysis methods, starting with the Morris method and then going on to the eFAST and Sobol method.

Mikhail Vaganov (@MikhailVaganov) is making good progress on his Nbody modeling language. This will make it easy to utilize DiffEq as a backend for molecular dynamics simulation. Follow the progress in DiffEqPhysics.jl
And here’s a quick view of the rest of our “in development” list:
 Preconditioner choices for Sundials methods
 Adaptivity in the MIRK BVP solvers
 More general Banded and sparse Jacobian support outside of Sundials
 IMEX methods
 Function input for initial conditions and time span (
u0(p,t0)
)  LSODA integrator interface
Projects
Are you a student who is interested in working on differential equations software and modeling? If so, please get in touch with us since we may have some funding after August for some student developers to contribute towards some related goals. It’s not guaranteed yet, but getting in touch never hurts!