faq.md

Frequently Asked Questions

This page is a compilation of frequently asked questions and answers.

Performance

Do you support GPUs? Multithreading? Distributed computation?

Yes. The *DiffEq.jl libraries (OrdinaryDiffEq.jl, StochasticDiffEq.jl, and DelayDiffEq.jl) are all written to be generic to the array and number types. This means they will adopt the implementation that is given by the array type. The in-place algorithms internally utilize Julia's broadcast (with some exceptions due to a Julia bug for now, see this issue) and Julia's A_mul_B! in-place matrix multiplication function. The out-of-place algorithms utilize standard arithmetical functions. Both additionally utilize the user's norm specified via the common interface options and, if a stiff solver, ForwardDiff/DiffEqDiffTools for the Jacobian calculation, and Base linear factorizations for the linear solve. For your type, you may likely need to give a better form of the norm, Jacobian, or linear solve calculations to fully utilize parallelism.

GPUArrays.jl, ArrayFire.jl, DistributedArrays.jl have been tested and work in various forms, where the last one is still not recommended for common use yet.

The next question is whether it matters. Generally, your system has to be large for parallelism to matter. Using a multithreaded array for broadcast we find helpful around N>1000, though the Sundials manual says N>100,000. For high order Runge-Kutta methods it's likely lower than the Sundials estimate because of more operations packed into each internal step, but as always that will need more benchmarks to be precise and will depend on the problem being solved. GPUs generally require some intensive parallel operation in the user's f function to be viable, for example a matrix multiplication for a stencil computation in a PDE. If you're simply solving some ODE element-wise on a big array it likely won't do much or it will slow things down just due to how GPUs work. DistributedArrays require parallel linear solves to really matter, and thus are only recommended when you have a problem that cannot fit into memory or are using a stiff solver with a Krylov method for the linear solves.

My ODE is solving really slow... what do I do?

First, check for bugs. These solvers go through a ton of convergence tests and so if there's a solver issue, it's either just something to do with how numerical methods work or it's a user-error (generally the latter, though check the later part of the FAQ on normal numerical errors). User-errors in the f function causing a divergence of the solution is the most common reason for reported slow codes.

If you have no bugs, great! The standard tricks for optimizing Julia code then apply. What you want to do first is make sure your function does not allocate. If your system is small (<=100 ODEs/SDEs/DDEs/DAEs?), then you should set your system up to use StaticArrays.jl. This is demonstrated in the ODE tutorial with static matrices. Static vectors/arrays are stack-allocated, and thus creating new arrays is free and the compiler doesn't have to heap-allocate any of the temporaries (that's the expensive part!). These have specialized super fast dispatches for arithmetic operations and extra things like QR-factorizations, and thus they are preferred when possible. However, they lose efficiency if they grow too large.

For anything larger, you should use the in-place syntax f(du,u,p,t) and make sure that your function doesn't allocate. Assuming you know of a u0, you should be able to do:

du = similar(u0)
@time f(du,u0,p,t)

and see close to zero allocations and close to zero memory allocated. If you see more, then you might have a type-instability or have temporary arrays. To find type-instabilities, you should do:

@code_warntype f(du,u,p,t)

and read the printout to see if there's any types that aren't inferred by the compiler, and fix them. If you have any global variables, you should make them const. As for allocations, some common things that allocate are:

• Array slicing, like u[1:5]. Instead, use @view u[1:5]
• Matrix multiplication with *. Instead of A*b, use A_mul_B!(c,A,b) for some pre-allocated cache vector c.
• Non-broadcasted expressions. Every expression on arrays should .= into another array, or it should be re-written to loop and do computations with scalar (or static array) values.

For an example of optimizing a function resulting from a PDE discretization, see this blog post.

The stiff solver takes forever to take steps for my PDE discretization! Why?

The solvers for stiff solvers require solving a nonlinear equation each step. In order to do so, they have to do a few Newton steps. By default, these methods assume that the Jacobian is dense, automatically calculate the Jacobian for you, and do a dense factorization. However, in many cases you may want to use alternatives that are more tuned for your problem.

First of all, when available, it's recommended that you pass a function for computing your Jacobian. This is discussed in the performance overloads section. Jacobians are especially helpful for Rosenbrock methods.

Secondly, if your Jacobian isn't dense, you shouldn't use a dense Jacobian! In the Sundials algorithm you can set linear_solver=:Band for banded Jacobians for example. More support is coming for this soon.

But lastly, you shouldn't use a dense factorization for large sparse matrices. Instead, if you're using a *DiffEq library you should specify a linear solver. For Sundials.jl, you should change the linear_solver option. See the ODE solve Sundials portion for details on that. Right now, Sundials.jl is the recommended method for stiff problems with large sparse Jacobians. linear_solver=:Band should be used if your Jacobian is banded and you can specify the band sizes. If you only know the Jacobian is sparse, linear_solver=:GMRES is a good option. Once again, a good reference for how to handle PDE discretizations can be found at this blog post.

Complicated Models

Can I switch my ODE function in the middle of integration?

There are a few ways to do this. The simplest way is to just have a parameter to switch between the two. For example:

function f(du,u,p,t)
  if p == 0
    du[1] = 2u[1]
  else
    du[1] = -2u[1]
  end
  du[2] = -u[2]
end

Then in a callback you can make the affect! function modify integrator.prob.p. For example, we can make it change when u[2]<0.5 via:

condition(t,u,integrator) = u[2] - 0.5
affect!(integrator) = integrator.prob.p = 1

Then it will change betweeen the two ODE choices for du1 at that moment. Another way to do this is to make the ODE functions all be the same type via FunctionWrappers.jl, but that is unnecessary. With the way that modern processors work, there exists branch prediction and thus execution of a conditional is free if it's predictable which branch will be taken. In this case, almost every call to f takes the p==0 route until the callback, at which point it is almost always the else route. Therefore the processor will effectively get rid of the computational cost associated with this, so you're likely over-optimizing if you're going further (unless this change happens every step, but even then this is probably the cheapest part of the computation...).

Numerical Error

The solver doesn't obey physical law X (e.g. conservation of energy)

Yes, this is because the numerical solution of the ODE is not the exact solution. There are a few ways that you can handle this problem. One way is to get a more exact solution. Thus instead of

sol = solve(prob,alg)

use

sol = solve(prob,alg,abstol=1e-10,reltol=1e-10)

Of course, there's always a tradeoff between accuracy and efficiency, so play around to find out what's right for your problem.

Another thing you can do is use a callback. There are some premade callbacks in the callback library which handle these sorts of things like projecting to manifolds and preserving positivity.

The symplectic integrator doesn't conserve energy?

Yes, symplectic integrators do not exactly conserve energy. It is a common misconception that they do. What symplectic integrators actually do is solve for a trajectory which rests on a symplectic manifold that is perturbed from the true solution's manifold by the truncation error. This means that symplectic integrators do not experience (very much) long time drift, but their orbit is not exactly the same as the true solution in phase space and thus you will see differences in energy that tend to look periodic. There is a small drift which grows linearly and is related to floating point error, but this drift is much less than standard methods. This is why symplectic methods are recommended for long time integration.

For conserving energy, there are a few things you can do. First of all, the energy error is related to the integration error, so simply solving with higher accuracy will reduce the error. The results in the DiffEqBenchmarks show that using a DPRKN method with low tolerance can be a great choice. Another thing you can do is use the ManifoldProjection callback from the callback library.

How do I get to zero error?

You can't. For floating point numbers, you shouldn't use below abstol=1e-14 and reltol=1e-14. If you need lower than that, use arbitrary precision numbers like BigFloats or ArbFloats.jl.

Autodifferentiation and Dual Numbers

Are the native Julia solvers compatible with autodifferentiation?

Yes! However, you should first look into sensitivity analysis as a (possibly in the future) more efficient method for calculating derivatives of the solution and functionals of the solution.

But if you still want to do it the naive way, here's how you do it. If the algorithm does not use adaptive time stepping, then you simply need to make the initial condition have elements of Dual numbers. If the algorithm uses Dual numbers, you need to make sure that time is also given by Dual numbers. A quick explanation of this is because changing the value of the initial condition will change the error in the steps, thus causing different steps to be taken changing the time values.

To show this in action, let's say we want to find the Jacobian of solution of the Lotka-Volterra equation at t=10 with respect to the parameters.

function func(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
function f(p)
  prob = ODEProblem(func,eltype(p).([1.0,1.0]),eltype(p).((0.0,10.0)),p)
  solve(prob,Tsit5(),save_everystep=false)[end]
end

This function takes in new parameters and spits out the solution at the end. We make the inital condition eltype(p).([1.0,1.0]) so that way it's typed to be Dual numbers whenever p is an array of Dual numbers, and we do the same for the timespan. Then we can take the Jacobian via ForwardDiff.jl:

using ForwardDiff
ForwardDiff.jacobian(f,[1.5,1.0])

2×2 Array{Float64,2}:
  2.20214   0.189782
 -6.2273   -0.700188

and compare it to Calculus.jl:

Calculus.jacobian(f,[1.5,1.0],:central)

2×2 Array{Float64,2}:
  2.20214   0.189782
 -6.2273   -0.700188

I get Dual number errors when I solve my ODE with Rosenbrock or SDIRK methods...?

This is because you're using a cache which is not compatible with autodifferentiaion via ForwardDiff.jl. For example, if we use the ODE function:

const tmp = zeros(4)
const A = rand(4,4)
function f(du,u,p,t)
  A_mul_B!(tmp,A,u)
  du .= tmp .+ u
end

Here we use a cached temporary array in order to avoid the allocations of matrix multiplication. When autodifferentiation occurs, the element type of u is Dual numbers, so A*u produces Dual numbers, so the error arises when it tries to write into tmp. There are two ways to avoid this. The first way, the easy way, is to just turn off autodifferentiation with the autodiff=false option in the solver. Every solver which uses autodifferentiation has this option. Thus we'd solve this with:

prob = ODEProblem(f,rand(4),(0.0,1.0))
sol = solve(prob,Rosenbrock23(autodiff=false))

and it will use a numerical differentiation fallback (DiffEqDiffTools.jl) to calculate Jacobians.

** Warning: Advanced **

The more difficult way is to create a Dual cache which dispatches for the cache choice based on the element type of u. This is done by the following:

using ForwardDiff
struct MyTag end
immutable DiffCache{T<:AbstractArray, S<:AbstractArray}
    du::T
    dual_du::S
end

function DiffCache{chunk_size}(T, size, ::Type{Val{chunk_size}})
    DiffCache(zeros(T, size...), zeros(ForwardDiff.Dual{nothing,T,chunk_size}, size...))
end

DiffCache(u::AbstractArray) = DiffCache(eltype(u),size(u),Val{ForwardDiff.pickchunksize(length(u))})

get_tmp{T<:ForwardDiff.Dual}(dc::DiffCache, ::Type{T}) = dc.dual_du
get_tmp(dc::DiffCache, T) = dc.du

Now we can get a cache that by dispatch either gives a cache array of Dual numbers or just floating point numbers:

const dual_cache = DiffCache(rand(4)) # Build the cache, this must match your IC
du = get_tmp(dual_cache,typeof(rand(4))) # Gives a Array{Float64}
dual_du = get_tmp(dual_cache,typeof(ForwardDiff.Dual(0.2,3.0))) # Gives Array{Dual}

Note that you have to make sure that your chunk size matches the choice in the ODE solver (by default, it uses ForwardDiff.pickchunksize(length(u)) as well, so you only need to change this if you explicitly set chunksize = ...). Now we can setup and solve our ODE using this cache:

function f(du,u,p,t)
  # Get du from cache
  tmp = get_tmp(dual_cache,eltype(u))
  # Fix tag
  _tmp = reinterpret(eltype(u),tmp)
  A_mul_B!(_tmp,A,u)
  du .= _tmp .+ u
end
prob = ODEProblem(f,rand(4),(0.0,1.0))
sol = solve(prob,Rosenbrock23())

Small explanation is in order. tmp = get_tmp(dual_cache,eltype(u)) makes tmp match u in terms of Dual or not, but it doesn't necessarily match the tag. So now we reinterpret our Dual array to put the right tag on there. Note that this simply changes type information and thus does not create any temporary arrays. Once we do that, our cached array is now typed correctly to hold the result of A_mul_B!.