DifferentialEquations.jl 2.0
This marks the release of ecosystem version 2.0. All of the issues got looked over. All (yes all!) of the API suggestions that were recorded in issues in JuliaDiffEq packages have been addressed! Below are the API changes that have occurred. This marks a really good moment for the JuliaDiffEq ecosystem because it means all of the long-standing planned API changes are complete. Of course new things may come up, but there are no more planned changes to core functionality. This means that we can simply work on new features in the future (and of course field bug reports as they come). A blog post detailing our full 2.0 achievements plus our 3.0 goals will come out at our one year anniversary. But for now I want to address what the API changes are, and the new features of this latest update.
The main changes in 2.0 is that this finishes a lot of the fundamental design of DiffEq. Many new problems were added, and the existing problems were modified to accommodate a much larger domain of problems. This leaves the next developments in the ecosystem to be fleshing out the algorithms for solving all of the problem types.
API Changes
This post is to detail the API changes that are happening in the DifferentialEquations.jl ecosystem. Most of them should not be very disruptive, but are things that you should take notice of.
Solution indexing change
The solution object now indexing has been changed to match Julian APIs better.
The linear index sol[i]
returning the value at the i
th timepoint is the same.
However, the “matrix indexing” has been flipped. Now, the j
th component at
timepoint i
is given as sol[j,i]
whereas before it was sol[i,j]
. This means that
sol[j,i]==sol[i][j]
. The reason for this is that in our setup “columns” are timepoints,
not across time. Given Julia’s column-major format, this gave the wrong impression
about contiguousness of the array.
Now, the DESolution
type is simply a subtype of the new AbstractVectorOfArray
type in RecursiveArrayTools.jl. That means that, in addition to this change,
many other features were added that make the solution type act more like an
array.
save_timeseries
changed to save_everystep
Many users were confused about what the save_timeseries
argument did, so it
was changed to save_everystep
. Now it’s clear that it means that every internal
solver step is counted for saving. This can still be controlled by timeseries_steps
to skip integer numbers of intervals.
saveat
=> no save_everystep
Before, save_everystep
always defaulted to true. However, this tripped up many
users who were interested in saving the output at specific timepoints: why would
the solver give me points I wasn’t asking for when I was already saying what I
wanted?
To fix this, the default for save_everystep
was changed to be whether isempty(saveat)
.
The result is that
sol = solve(prob,alg,saveat=save_timepoints)
is equivalent to what used to be
sol = solve(prob,alg,saveat=save_timepoints,save_everystep=false)
Note that you can recover the previous default behavior by explicitly setting
save_everystep=true
:
sol = solve(prob,alg,saveat=save_timepoints,save_everystep=true)
We believe that this default should be more sensible to most users.
saveat
using Numbers
A common usage for saveat
is to get out an evenly-spaced grid. For example,
0.0:0.1:1.0
if you want 10 evenly-spaced points. Actually, that has some floating
point issues, so we thought it would be better to help you out.
sol = solve(prob,alg,saveat=0.1)
will now make the solver save at tspan[1]:0.1:tspan[2]
timepoints. This should
be a nice shorthand for a very standard operation.
save_idxs
A new argument was added which allows you to choose which indices to save. So far
this is only implemented in the *DiffEq
packages (OrdinaryDiffEq.jl, StochasticDiffEq.jl,
DelayDiffEq.jl), but can be made compatible with other solvers as well (please open
and issue if you would like to have this in Sundials.jl or other packages. Or a
PR! This isn’t a difficult change!). For example:
sol = solve(prob,alg,saveat=0.1,save_idxs=1:2:5)
will save an array of three values (dependent variables 1, 3, and 5) at each
0.1
points in time. Note this can be used without saveat
.
save_start
Another new output control. save_start
notes whether the initial condition should
be appended to the start of the solution type. This defaults as true, and most of
the solver packages implement this option.
Callbacks in problems
The problem types now can hold callbacks. This is to associate some behavior like
events with the problem instead of the integrator. They do not act differently,
and are just a callback
keyword argument in the problem types.
ParameterizedFunctions
Although there was discussions about changing the @ode_def
macro syntax, I ultimately
decided against it. Instead, I polished up the edges of ParameterizedFunctions
and added to the documentation how you can do all of the things which users
had issues about. Some of the functionality in there is new. There was a large
change and an upgrade in the version of SymEngine that is used, so now more
functions can be differentiated. Importantly, no @eval
s are used anymore, and
this means that in the future all of DifferentialEquations.jl can be statically
compiled. But now this means that DifferentialEquations.jl can be precompiled.
As for where this is going next, ParameterizedFunctions is going to maintenance and
bugfix mode. It seems everything that can be done in that architecture is done.
All of the new cool ideas will become an @diffeq
DSL in DiffEqDSL.jl. You
can read the following issues for the plan:
https://github.com/JuliaDiffEq/DiffEqDSL.jl/issues/1 https://github.com/JuliaDiffEq/ParameterizedFunctions.jl/issues/17
It will be exciting, but will be separate from the current macro and thus you have no reason to expect breakage.
New Features
EulerHeun and RKMil interpretation
The EulerHeun
method is our first SDE method for Stratonovich differential equations.
We note that Ito vs Stratonovich is a difference in the interpretation of the integral,
not in the SDE itself. Thus there is no switch for “choosing” Ito vs Stratonovich
(or whatever else), instead it’s just noted in the docs for which interpretation
the integration technique converges.
For some methods, there are ways to make a different version compatible with different
interpreations. So for example, we have in development an interpretation
syntax
for RKMil
, where RKMil(interpretation=:Stratonovich)
converges with first-order
in the Stratonovich sense (defaults to :Ito
). We hope to extend our support as
time progresses.
Non-Diagonal Noise
We finally accept a form of non-diagonal noise. For SDEs, you can give a
noise_rate_prototype
which will be the array for the du
slot in the SDE’s
noise function g
. This defaults to nothing
, which will make du
and array
whose geometry matches u
and does diagonal noise, that is du.*dW
. However,
when you specify du
, it changes it do the standard form for noise: du*dW
.
This allows for commutative and non-commutative noise by specifying the array
to be a matrix like noise_rate_prototype=rand(n,m)
where m
is the Ito dimension
(number of independent Brownian motions), and it will automatically generate dW
as a vector of size m
.
Importantly, this works on any AbstractMatrix
. So for example, if your reactions
are sparse on the vector dW
, you can choose noise_rate_prototype=sparse(A)
,
some sparse matrix A
which has the correct set of non-zeros for your problem,
and then du
will be that sparse matrix. Even special forms for the matrix,
like Tridiagonal
, are accepted. This means you can specify this in a way that
is suitable to your problem yet retain efficiency by using specialized matrix
types.
Currently only EM
and EulerHeun
are compatible with non-diagonal noise,
but this area should be growing.
NoiseProcess Improvements
A NoiseProcess
is how you create colored noise in the SDE and RODE problems.
In fact, it’s general enough to create any noise process! For example,
WienerProcess(t0,W0,Z0=nothing)
WienerProcess!(t0,W0,Z0=nothing)
are used to build a WienerProcess
either in its inplace form or in its out-of-place
form. The resulting type W
is designed to use the RSwM adaptivity algorithms to
allow for very fast generation of the processes. In addition, the process is a
condition function W(t)
which will automatically interpolate at the necessary
points. This allowed for the creation of a wrapper
NoiseWrapper(W::NoiseProcess)
that lets you re-use the same noise process in a different stochastic simulation. In addition, you can also create spatial colored noise easily by passing in a constant covariance matrix:
CorrelatedWienerProcess(Γ,t0,W0,Z0=nothing)
CorrelatedWienerProcess!(Γ,t0,W0,Z0=nothing)
The noise processes now have their own package DiffEqNoiseProcess.jl and the interface is documented so that way the creation of new types of noise processes can be done as needed.
RODEs
Random Ordinary Differential Equations (RODEs) are a rapidly growing area where,
for some stochsatic process y(t)
, you have a differential equation:
u' = f(t,u,y)
It can be shown that, under certain assumptions, SDEs are a form of RODE. There is
a new problem type and our first integrator, RandomEM()
, which is the random
Euler-Maruyama method. If no NoiseProcess
is given, then y(t)
defaults to
white noise, and this converges to the Ito interpretation of the random integral.
More methods to come. One thing which will help will be the release of Kloeden’s new book on methods for RODEs. This has the same compatibility throughout the ecosystem as SDEs, so it works for things like parameter estimation. We hope this exciting and one of a kind feature helps your research take new and interesting directions!
Steady State Calcuations
In many cases you want an easy way to calculate the steady states of a differential
equation, that is the state u
such that
0 = f(u)
The new library DiffEqSteadyState.jl has solvers which make this simple. You can
generate a new steady state problem from an f
and an inital guess:
SteadyStateProblem(f,u0,mass_matrix=I)
or directly convert an ODEProblem
into a SteadyStateProblem
:
SteadyStateProblem(prob::ODEProblem)
This can then be solved using the new libraries solve function dispatch. Right now
only SSRootfind
is implemented, which will solve for the steady states using
a rootfinding algorithm. However, we plan to implement accelerated methods specific
for differential equations (and partial differential equations) that will speed
up these calculations. In fact, a new algorithm for stochastic steady states
may be coming as a new publication, with an implmentation here…
Complex Plotting
This new update adds the recipe library DimensionalPlotRecipes.jl to the fray. This allows for complex numbers to be directly plotted, and many other controls are provided. We plan to extend this to work with Quaternions and allow for automatic dimensional reduction for high dimensional plots. Documentation for this can be found at the library’s README and is now linked to in the DiffEq plotting section.
Monte Carlo Parallelism Types
The Monte Carlo functionality now has support for different parallelism types,
so you can tell it to parallelize using threading, pmap, @parallel
, or no
parallelism at all. Additionally, there is a split_threads
option which will
use threading separately on each process, allowing one to using threading on every
node of a cluster for maximum performance.
“Refined Problem Types”
This update introduces “refined problem types”. These are problem types which give
the solver additional information. For example, a SplitODEProblem
is of the form:
du = f1(t,u) + f2(t,u) + ... + fn(t,u)
and the solvers can use the information from this split problem to have more optimized algorithms. Additionally, a partitioned form
du1 = f1(t,u1,u2,...,uN)
du2 = f2(t,u1,u2,...,uN)
...
dun = fn(t,u1,u2,...,un)
which can be used for things like symplectic algorithms. The underlying machinery
is the new type ArrayPartition
in RecursiveArrayTools.jl and allows for heterogeneous
types to be used with type-stable broadcasting. Notably, this means that second
order problems:
u'' = f(t,u,u')
can be solved with the proper units without loss of efficiency (once broadcast
is used internally). Thus there are problem types for these new problems, including
a helper type for SecondOrderODEProblem
.
While there are not many algorithms implemented which use all of this yet, the machinery is all in place in OrdinaryDiffEq.jl, DelayDiffEq.jl, and StochasticDiffEq.jl, so implementing new algorithms which solve these types of problems is now easy. This will be the basis of some PDE algorithms as well.
Mass Matrices
Mass matrices were added to the appropriate problem types. Most of the solvers don’t support using them yet, but now all of the solvers have a way of recieving a user-defined mass matrix. Updates will come which will then make use of this functionality.
linsolve/nlsolve Choice
Now one can explicitly choose the linear solver and nonlinear solver functions
for each of the *DiffEq
solvers. This documentation page
explains how to do this. This means you can tell the linear solving to occur on
the GPU, or using PETSc, etc., and replace the nonlinear solver code with one of
your own choosing.
Full v0.6 Compatibility
DifferentialEquations.jl is now compatible with Julia’s v0.6. Additionally, the deprecation warnings have all been cleaned up (except for the deprecation warnings from the dependencies NLsolve.jl and IterativeSolvers.jl… go bug them to merge my PRs and tag a new version if you want those depwarns gone :)). This means that you should feel free to start using DifferentialEquations on v0.6 without a hitch.
But there is one caveat to mention. Plots.jl is not v0.6 compatible right now, and so you should stay on v0.5 until plotting is ready if that’s a functionality that you need.
Near Future Changes
These are changes which didn’t quite make the 2.0 release, but will be coming soon after.
Full Precompilation
Before, DifferentialEquations.jl could not be precompiled because anything related to ParameterizedFunctions.jl could not precompiled. This was due to the dynamic use of SymEngine.jl. However, SymEngine will release an update which makes this no longer necessary, and thus all of this can now be precompiled. The result is that the entirety of DiffEq will be precompile friendly. Not only that, it will be statically compliable.
Experimental: DiffEqIO
IO functionality has been experimentally added through IterableTables.jl. This allows one to easily save solution types to DataFrames, CSVs, and more. This will be added to the documentation when it’s out of the experimental phase.
Extensive PDE Tools
This will be detailed in the follow-up post on the state of the ecosystem.