The Unique Features and Performance of DifferentialEquations.jl¶
### JuliaCon 2017
### Chris Rackauckas ### University of California, Irvine
Quick About Me¶
- PhD Candidate in Mathematics at UC Irvine
- Research in:
- Numerical methods for stochastic ODEs/PDEs
- Multiscale and spatial biological simulations
- Effects of stochasticity in gene regulatory networks
- I created DiffEq because it was necessary:
- No software was using the latest methods
- I wanted the latest and fastest methods to benchmark against
- I am too lazy to change my code around: everything needs one interface
- I wanted to be able to use these methods to build large biological models
- I wanted parallel computing to be part of the design
- I wanted tools to be able to match models to data
DifferentialEquations.jl is about modeling¶
Differential Equations vs Machine Learning¶
- The models are different
- Machine learning uses general non-mechanistic models
- Differential equations utilize knowledge about the structure of the problem
- The reliance on data is different
- Machine learning needs to learn from data
- Machine learning models stack simple models together
- DE models are complex and nonlinear
- DE can model without data, or can calibrate to data
- DE have fewer parameters, and can be learned with less data
- Both are used to create predictions
- Ordinary and partial differential equations are the standard tools of physics
- Stochastic differential equations are a standard tool in stock market analysis
There are many types of differential equations¶
- Ordinary differential equations evolve over one variable (let's call this time)
- Systems of differential equations evolve multiple quantities at the same time
- Stochastic differential equations add randomness (referred to as "noise")
- Delay differential equations allow for delayed interactions
- Gillespie simulations evolve discrete quantities via a stochastic model
- Partial differential equations evolve over multiple variables (space and time)
DifferentialEquations.jl incorporates all of these under a single interface
Problem 1: Ordinary Differential Equations¶
http://docs.juliadiffeq.org/latest/tutorials/ode_example.html
du(t)dt=f(u(t)p,t)where u(t0)=u0 is the known initial value. We want to know u(t) for t∈[t0,tf]
Exponential Growth Model¶
- My money has been temporarily stolen by the banking cartel (it's in the bank)
- My money grows at a rate of 98% per year due to interest (really good bank)
- I start with $1: How much will I have in 5 years?
Solve for u(t) for t∈[0,5]
In [1]:
using DifferentialEquations
f(u,p,t) = p*u
u0=1.0
p = 0.98
tspan = (0.0,5.0)
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob);
In [2]:
using Plots; gr()
# Use the plot recipe, along with some plotting attributes
plot(sol,linewidth=5,title="My Bank Account",
xaxis="Time (t)",yaxis="Money!",legend=false)
Out[2]:
Controlling the Solvers¶
See http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html
In [3]:
sol = solve(prob)
println(length(sol.t))
sol = solve(prob,abstol=1e-6,reltol=1e-6)
println(length(sol.t))
sol = solve(prob,saveat=0.5,reltol=1e-6)
println(sol.t)
sol = solve(prob,tstops=[0.5],reltol=1e-4)
println(sol.t)
sol = solve(prob,reltol=1e-6,save_everystep=false)
println(sol.t)
11 15 [0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0] [0.0,0.0634774,0.221501,0.43294,0.5,0.696266,0.923443,1.20696,1.53056,1.89898,2.30449,2.74501,3.21551,3.71248,4.23208,4.77109,5.0] [0.0,5.0]
Analyzing the Solution¶
In [4]:
sol = solve(prob,Tsit5())
println(sol.t[5]) # 5th timepoint
println(sol[5]) # 5th timepoint value
println([t+2u for (t,u) in zip(sol.t,sol.u)])
1.106052911995544 2.956279956560924 [2.0,2.30727,3.17657,4.63952,7.01861,11.2071,18.9435,34.2493,66.3366,137.407,273.572]
DifferentialEquations.jl comes with interpolations on the solution type
In [5]:
sol(0.45) # The value of the solution at t=0.45
Out[5]:
1.5542610480525965
Choosing Solvers¶
http://docs.juliadiffeq.org/latest/solvers/ode_solve.html
DifferentialEquations.jl is composed of many (30+?) packages which together provide its functionality under one interface. For Ordinary Differential Equations, we have:
- OrdinaryDiffEq.jl: native Julia library. Fully featured, and in many cases faster than the classic C/Fortran libraries.
- Sundials.jl: provides Adams and BDF methods. A standard library
- LSODA.jl: provides the stability-detecting LSODA algorithm
- ODEInterface.jl: provides classic Fortran algorithms like
dopriandradau - ODE.jl: backwards compatibility
DifferentialEquations.jl also provides automatic detection of the appropriate algorithm
In [6]:
# Automatically uses CVODE_BDF from Sundials
sol = solve(prob,alg_hints=[:stiff],reltol=1e-8,abstol=1e-8)
# Use Tsit5 from OrdinaryDiffEq
sol = solve(prob,Tsit5());
# Use CVODE_BDF from Sundials
sol = solve(prob,CVODE_BDF());
using ODEInterfaceDiffEq
# Use radau from ODEInterfaceDiffEq
sol = solve(prob,radau());
The *DiffEq Libraries¶
The *DiffEq libraries (OrdinaryDiffEq.jl, StochasticDiffEq.jl, DelayDiffEq.jl) are the native Julia solvers. They have a lot of unique features:
- They routinely benchmark as the fastest implementations
- They have cheap/free high order interpolations (will show in a second!)
- They have expressive event handling
- They can handle generic
NumberandAbstractArraytypes
The wrapped libraries, while classics (Sundials.jl, ODEInterfaceDiffEq.jl) are restricted to Vector{Float64}.
http://docs.juliadiffeq.org/latest/basics/compatibility_chart.html
Method Recommendations¶
BS3()for fast low accuracy non-stiffTsit5()for standard non-stiff. This is the first algorithm to try in most cases.Vern7()for high accuracy non-stiffRosenbrock23()for stiff equations with Julia-defined types, events, etc.radau()for stiff equations at higher accuracy onVector{Float64}CVODE_BDF()for large stiff equations (discretizations of PDEs) onVector{Float64}
Translations from MATLAB/Python/R¶
ode23-->BS3()ode45/dopri5-->DP5(), though in most casesTsit5()is more efficientode23s-->Rosenbrock23()ode113-->CVODE_Adams(), though in many casesVern7()orVern8()are more efficientdop853-->DP8(), though in most casesVern7()orVern8()are more efficientode15s/vode-->CVODE_BDF(), though in many casesradau()is more efficientode23t-->Trapezoid()lsoda-->lsoda()(requiresPkg.add("LSODA"); using LSODA)ode15i-->IDA()
Example 2: Systems of Differential Equations¶
dxdt=σ(y−x)dydt=x(ρ−z)−ydzdt=xy−βzIn [ ]:
function lorenz(du,u,p,t)
du[1] = 10.0(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,25.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,Tsit5());
The Plot Recipe Interpolates and Does Phase Plots¶
In [4]:
plot(sol,vars=(1,2,3))
Out[4]:
In [ ]:
plot(sol,vars=(1,2,3),denseplot=false)
ParameterizedFunctions¶
dxdt=σ(y−x)dydt=x(ρ−z)−ydzdt=xy−βzIn [10]:
g = @ode_def LorenzExample begin
dx = σ*(y-x)
dy = x*(ρ-z) - y
dz = x*y - β*z
end σ ρ β
p = [10.0,28.0,8/3]
Out[10]:
(::LorenzExample) (generic function with 15 methods)
DifferentialEquations.jl respects the user's input type¶
In [11]:
using StaticArrays, DifferentialEquations
A = big.(@SMatrix [ 1.0 0.0 0.0 -5.0
4.0 -2.0 4.0 -3.0
-4.0 0.0 0.0 1.0
5.0 -2.0 2.0 3.0])
u0 = big.(@SMatrix rand(4,2))
tspan = (0.0,1.0); f_SA(u,p,t) = A*u
prob = ODEProblem(f_SA,u0,tspan);
sol = solve(prob,Tsit5()); plot(sol)
Out[11]:
In [12]:
using Unitful
f_units(y,p,t) = 0.5*y / 3.0u"s"
u0 = 1.0u"N"
prob = ODEProblem(f_units,u0,(0.0u"s",1.0u"s"))
sol =solve(prob,Tsit5());
println(sol.t[end])
println(sol[end])
sol(0.5u"s")
1.0 s 1.1813604129914028 N
Out[12]:
1.0869040537266441 N
In [13]:
using DifferentialEquations, SymEngine
y0 = symbols(:y0)
u0 = y0
f_sym(u,p,t) = 2u
prob = ODEProblem(f_sym,u0,(0.0,1.0))
sol = solve(prob,RK4(),dt=1/10);
sol = solve(prob,RK4(),dt=1/1000)
end_solution = lambdify(sol[end])
println(end_solution(2)) # 14.778112197857341
println(2exp(2)) # 14.7781121978613
println(end_solution(3)) # 22.167168296786013
println(3exp(2)) # 22.16716829679195
14.778112197857341 14.7781121978613 22.167168296786013 22.16716829679195
DifferentialEquations.jl is unique when speed matters¶
DifferentialEquations.jl is a research lab which includes many new methods
In [14]:
using DifferentialEquations, Plots, ODEInterfaceDiffEq, ODE
# 2D Linear ODE
function test_f(du,u,p,t)
for i in 1:length(u)
du[i] = 1.01*u[i]
end
end
function test_f(::Type{Val{:analytic}},u₀,p,t)
u₀*exp(1.01*t)
end
tspan = (0.0,10.0)
prob = ODEProblem(test_f,rand(100,100),tspan)
abstols = 1./10.^(3:13)
reltols = 1./10.^(0:10);
setups = [Dict(:alg=>DP5())
Dict(:alg=>ode45())
Dict(:alg=>dopri5())
Dict(:alg=>Tsit5())]
names = ["DifferentialEquations";"ODE";"ODEInterface";"DifferentialEquations Tsit5"]
wp = ode_workprecision_set(prob,abstols,reltols,setups;names=names,save_everystep=false);
In [15]:
plot(wp)
Out[15]:
Callbacks and Event Handling¶
You can directly apply dynamic behavior to an ODE through its callbacks and event handling
In [16]:
function b_condition(u,t,integrator) # Event when event_f(t,u) == 0
u[1]
end
function b_affect!(integrator)
integrator.u[2] = -integrator.u[2]
end
cb = ContinuousCallback(b_condition,b_affect!)
bb = @ode_def_bare BallBounce begin
dy = v
dv = -g
end g=9.81
u0 = [50.0,0.0]
tspan = (0.0,15.0)
prob = ODEProblem(bb,u0,tspan)
sol = solve(prob,Tsit5(),callback=cb);
In [17]:
plot(sol,plotdensity=1000)
Out[17]:
Pretty much anything can be changed with the integrator¶
Let's model a growing cell population
- Each cell has an internal protein
X - When the concentration of
Xhits 1, the cell divides - When there are more cells, there is less food and so
Xis created less rapidly
In [18]:
using DifferentialEquations, Plots
function growth(du,u,p,t)
for i in 1:length(u)
du[i] = 0.3u[i]/length(u)
end
end
function condition(u,t,integrator) # Event when event_f(t,u) == 0
1-maximum(u)
end
function affect!(integrator)
u = integrator.u
resize!(integrator,length(u)+1)
maxidx = findmax(u)[2]
Θ = rand()
u[maxidx] = Θ
u[end] = 1-Θ
length(u) == 10 && terminate!(integrator)
nothing
end
callback = ContinuousCallback(condition,affect!)
u0 = [0.2]
tspan = (0.0,100.0)
prob = ODEProblem(growth,u0,tspan)
sol = solve(prob,callback=callback);
In [19]:
plot(sol.t,map((x)->length(x),sol[:]),lw=3,
ylabel="Number of Cells",xlabel="Time")
Out[19]:
In [20]:
ts = linspace(0,sol.t[end],1000)
plot(ts,map((x)->x[1],sol.(ts)),lw=3,
ylabel="Amount of X in Cell 1",xlabel="Time")
Out[20]:
This flexibility takes DiffEq from equation solver to simulation engine¶
- Abstract types can be used to encode complex models
- Integrators and event handling gives dynamic interactions
- Speed, control, adaptivity, plotting, etc. then all comes for free!
This leads to more sustainable scientific software development:
- Many scientific projects develop their own internal solvers for their models
- This requires a lot of extra time!
- How well tested is it? How well optimized is it? How maintainable is it?
- Is it using all of the latest PI-stepsize adaptivity? Verner's coefficients? etc.?
- DifferentialEquations.jl is thus ideal as a core engine in complex models
- Developer-friendly dependency management, MIT licensed
- Active developer channels (Gitter, or anywhere...)
Example: Multiscale Biological Models¶
In [21]:
using MultiScaleArrays
using OrdinaryDiffEq, DiffEqBase, StochasticDiffEq
# Define a hierarchy
immutable Cell{B} <: AbstractMultiScaleArrayLeaf{B}
x::Vector{B}
end
immutable Population{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
x::Vector{T}
y::Vector{B}
end_idxs::Vector{Int}
end
immutable Tissue{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
x::Vector{T}
y::Vector{B}
end_idxs::Vector{Int}
end
immutable Embryo{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
x::Vector{T}
y::Vector{B}
end_idxs::Vector{Int}
end
In [22]:
cell1 = Cell([1.0;2.0;3.0])
cell2 = Cell([4.0;5])
population = construct(Population,deepcopy([cell1,cell2])) # Make a Population from cells
cell3 = Cell([3.0;2.0;5.0])
cell4 = Cell([4.0;6])
population2 = construct(Population,deepcopy([cell3,cell4]))
tissue1 = construct(Tissue,deepcopy([population,population2])) # Make a Tissue from Populations
tissue2 = construct(Tissue,deepcopy([population2,population]))
embryo = construct(Embryo,deepcopy([tissue1,tissue2])); # Make an embryo from Tissues
In [23]:
function cell_ode(dcell,cell,p,t)
m = mean(cell)
for i in eachindex(cell)
dcell[i] = -m*cell[i]
end
end
function embryo_ode(dembryo,embryo,p,t)
for (cell,y,z) in LevelIterIdx(embryo,2)
cell_ode(@view dembryo[y:z],cell,p,t)
end
end
# Have a cell divide at t=0.5
tstop = [0.5]
function grow_condition(u,t,integrator)
t ∈ tstop
end
function grow_affect!(integrator)
add_daughter!(integrator,integrator.u[1,1,1],1,1)
end
growing_cb = DiscreteCallback(grow_condition,grow_affect!)
prob = ODEProblem(embryo_ode,embryo,(0.0,1.0));
sol = solve(prob,Tsit5(),callback=growing_cb,tstops=tstop)
sol[end]
Out[23]:
Embryo{Tissue{Population{Cell{Float64},Float64},Float64},Float64}(Tissue{Population{Cell{Float64},Float64},Float64}[Tissue{Population{Cell{Float64},Float64},Float64}(Population{Cell{Float64},Float64}[Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.290911,0.581821,0.872732]),Cell{Float64}([1.16364,1.45455]),Cell{Float64}([0.0,0.0,0.0])],Float64[],[3,5,8]),Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.600023,0.400015,1.00004]),Cell{Float64}([0.80003,1.20005])],Float64[],[3,5])],Float64[],[8,13]),Tissue{Population{Cell{Float64},Float64},Float64}(Population{Cell{Float64},Float64}[Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.600023,0.400015,1.00004]),Cell{Float64}([0.80003,1.20005])],Float64[],[3,5]),Population{Cell{Float64},Float64}(Cell{Float64}[Cell{Float64}([0.250004,0.500008,0.750013]),Cell{Float64}([1.00002,1.25002])],Float64[],[3,5])],Float64[],[5,10])],Float64[],[13,23])
Stochastic Differential Equations¶
du=f(t,u)dt+∑igi(t,u)dWi- The f term is the deterministic or "drift" term
- The g term is the noise or "diffusion term
Well-known model: the Black-Scholes Model
- My stock is growing at an average rate of 101%
- However, the value of my money has fluctuatations of 87% of its current value
In [19]:
using DifferentialEquations, Plots; gr()
α=1.01; β=0.87; u₀=1.0
bs_f(u,p,t) = α*u
bs_g(u,p,t) = β*u
bs_f(::Type{Val{:analytic}},u₀,p,t,W) = u₀*exp((α-(β^2)/2)*t+β*W)
prob = SDEProblem(bs_f,bs_g,u₀,(0.0,1.0))
sol = solve(prob,EM(),dt=1/2^(4))
plot(sol,plot_analytic=true)
Out[19]:
Fast adaptive methods are unique to DifferentialEquations.jl¶
In [25]:
sol =solve(prob,SRIW1())
plot(sol,plot_analytic=true)
Out[25]:
They tend to be orders of magnitude faster than standard methods¶
Many models cannot be solved without these new methods. Even those that can are faster with them.
Gillespie Simulations¶
Gillespie simulations are discrete models which change stochastically over time.
- A protein binds to DNA at a rate X
- A protein is created at a rate Y
These types of models are related to differential equations
- When the number of objects gets high, it converges for an SDE on the concentrations
- When the number gets really high, it converges to a deterministic ODE
In [10]:
prob = DiscreteProblem([100,1,0],(0.0,100.0))
rate = (u,p,t) -> (0.1/100.0)*u[1]*u[2]
jump_affect! = function (integrator)
integrator.u[1] -= 1
integrator.u[2] += 1
end
jump = ConstantRateJump(rate,jump_affect!)
rate = (u,p,t) -> 0.01u[2]
jump_affect! = function (integrator)
integrator.u[2] -= 1
integrator.u[3] += 1
end
jump2 = ConstantRateJump(rate,jump_affect!)
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
Out[10]:
DiffEqJump.JumpProblem{P,A,C,J<:Union{DiffEqJump.AbstractJumpAggregator,Void},J2} with problem DiffEqBase.DiscreteProblem{uType,tType,isinplace,F,C} and aggregator DiffEqJump.Direct
Number of constant rate jumps: 2
Number of variable rate jumps: 0
In [11]:
sol = solve(jump_prob,Discrete(apply_map=false));
plot(sol)
Out[11]:
Mix jumps and differential equations¶
In [ ]:
diff_f = function (du,u,p,t)
du[4] = u[2]*u[3]/1000 - u[1]*u[2]/100000
end
diff_g = function (du,u,p,t)
du[4] = 0.1u[4]
end
prob = SDEProblem(diff_f,diff_g,[100.0,1.0,0.0,1.0],(0.0,100.0))
jump_prob = JumpProblem(prob,Direct(),jump,jump2)
sol = solve(jump_prob,SRIW1(),abstol=1e-1,reltol=1e-1); plot(sol)
These features together allow for expressive scientific modeling¶
Let's go back to the multiscale embyro model:
- We can have cell proteins change according to a stochastic differential equation
- A rate of division can be dependent on a specific "growth factor"
- "Jumps" using this rate can cause cell division, growing our multiscale model
- We get high performance algorithms and optimized implementations the whole way down.
<3 Julia
Addons¶
Up until now, I have shown how to solve many different equations and encode a wide variety of models to simulate using the various solvers in DifferentialEquations.jl. But because all of this uses a single interface, an addon suite was able to be developed.
- Tools for automatic parallelization of Monte Carlo simulations
- Uncertainty quantification
- Parameter estimation and data fitting
Monte Carlo Simulations¶
Many times you need to solve the same equation multiple times. This is especially true for stochastic simulations.
In [29]:
pf = function (du,u,p,t)
du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
du[2] = -3 * u[2] + u[1]*u[2]
end
pf = ParameterizedFunction(pf_func,)
pg = function (du,u,p,t)
du[1] = p[3]*u[1]
du[2] = p[4]*u[2]
end
p = [1.5,1.0,0.1,0.1]
prob = SDEProblem(pf,pg,[1.0,1.0],(0.0,10.0),p)
prob_func = function (prob,i)
problem_new_parameters(prob,[1.5,1.0,0.3rand(),0.3rand()])
prob
end
monte_prob = MonteCarloProblem(prob,prob_func=prob_func)
sim = solve(monte_prob,SRIW1(),num_monte=10);
In [30]:
plot(sim,linealpha=0.6,color=:blue,vars=(0,1),title="The Timeseries'")
plot!(sim,linealpha=0.6,color=:red,vars=(0,2))
Out[30]:
In [31]:
summ = MonteCarloSummary(sim,0:0.1:10)
gr() # Note that plotly does not support ribbon plots
plot(summ,fillalpha=0.4)
Out[31]:
Uncertainty Quantification¶
Estimate the uncertainty in the solution due to numerical error
In [32]:
fitz = @ode_def FitzhughNagumo begin
dV = c*(V - V^3/3 + R)
dR = -(1/c)*(V - a - b*R)
end a b c
p = [0.2,0.2,3.0]
u0 = [-1.0;1.0]
tspan = (0.0,20.0)
prob = ODEProblem(fitz,u0,tspan,p)
cb = AdaptiveProbIntsUncertainty(5)
monte_prob = MonteCarloProblem(prob);
In [33]:
sim = solve(monte_prob,Tsit5(),num_monte=100,
callback=cb,abstol=1e-3,reltol=1e-1);
using Plots; gr(); plot(sim,vars=(0,1),linealpha=0.4)
Out[33]:
In [34]:
sim = solve(monte_prob,Tsit5(),num_monte=100,
callback=cb,abstol=1e-6,reltol=1e-3)
plot(sim,vars=(0,1),linealpha=0.4)
Out[34]:
Case Study: Parameter Estimation¶
Instead of simply solving differential equation models, one can ask the inverse problem: given that I see data like y, what needs to be true about the model such that we can match the data?
The act of estimating the parameters for a mechanistic differential equation model is called parameter estimation.
- Parameter estimation required repeated solving of differential equations
- Solving equations must be fast
- It requires access to optimization routines (global optimization is necessary)
- Machine learning tools (regularization, cross validation, etc.) are a must
Thus DifferentialEquations.jl is well-suited for handling these kinds of problems.
Problem Setup¶
In [1]:
using DifferentialEquations, Plots; gr()
lotka = @ode_def_nohes LotkaVolterraTest begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a b c d
p = [1.5,1.0,3.0,1.0]
u0 = [1.0;1.0]; tspan = (0.0,10.0);
prob = ODEProblem(lotka,u0,tspan,p)
sol = solve(prob,Vern8(),abstol=1e-14,reltol=1e-14);
t = collect(linspace(0,10,200))
randomized = [(sol(t[i]) + .01randn(2)) for i in 1:length(t)]
using RecursiveArrayTools;
data = vecvec_to_mat(randomized);
In [2]:
plot(sol,plotdensity=1000); scatter!(t,data)
Out[2]:
We can tell it to build an objective function that internally uses the DiffEq solver¶
This objective function has dispatches to work with all of the common optimization packages in Julia (JuMP/MathProgBase, Optim, BlackBoxOptim, ...)
In [3]:
obj = build_loss_objective(prob,Tsit5(),CostVData(t,data),maxiters=10000,verbose=true)
Out[3]:
(::DiffEqObjective) (generic function with 2 methods)
Let's recover the parameters using the NLopt.jl package
In [4]:
import NLopt
opt = NLopt.Opt(:LN_BOBYQA, 4)
NLopt.lower_bounds!(opt,zeros(4))
NLopt.upper_bounds!(opt,5ones(4))
NLopt.min_objective!(opt, obj)
NLopt.maxeval!(opt, Int(1e5))
(minf,minx,ret) = NLopt.optimize(opt,ones(4))
Out[4]:
(0.004039644894588035,[1.5,0.999476,2.99947,0.99869],:ROUNDOFF_LIMITED)
What's Next?¶
The next focus for DifferentialEquations.jl development is PDEs and stiff problems
- Native Julia stiff solvers are currently a lower order of accuracy then the top methods (@miguelraz's GSoC project will help here)
- Parallel-in-time Neural Network methods for large systems of stiff ODEs/SDEs (@akaysh's GSoC project)
- Fast matrix-free operators for FDM derivative discretizations (@shivin9's GSoC project)
- A wrapper for FEniCS's FEM toolbox (@ysimillides's GSoC project)
- High order Rosenbrock methods
- Stiffness detection and switching (already implemented, but no default methods)
Other areas which are under development are:
- New high order (+ adaptive) methods for SDEs will be published soon (orders of magnitude improvement)
- Integration of regularization, maximum likelihood, and Bayesian (MCMC) techniques for parameter estimation (@ayush-iitkgp's GSoC project)
- Boundary value problem solvers (@YingboMa's GSoC project)
- More geometric and symplectic integrators
- Integration of DiffEq with applications (large scale biological models, Pk/Pd estimation, tools for physical and financial modeling)
Acknowledgments¶
I am deeply indebted to every JuliaDiffEq contributor. I would like to especially thank:
- Gabriel Gellner
- Ethan Levien
- Scott P. Jones
- Lyndon White
- @finmod
- @dextorious
- Viral Shah
- Vijay Ivaturi
- Andreas Rossler
- Ernst Hairer
- Lawrance F. Shampine
The Unique Features and Performance of DifferentialEquations.jl ### JuliaCon 2017 ### Chris Rackauckas
### University of California, Irvine