ChrisRackauckas

Brusselator Stiff Partial Differential Equation Benchmark: Julia DifferentialEquations.jl vs Python SciPy

Tested is DifferentialEquations.jl vs Python's SciPy ODE solvers. Notes:

• Stiff ODE solvers are used since they are required to solve this problem effectively.
• The Python code is vectorized with for maximum performance
• All of the performance features are tried: automatic sparsity detection, preconditioners, etc.

Results

ChrisRackauckas

retcode: Success
Interpolation: 3rd order Hermite
t: [0.0, 0.5, 1.0]
u: Vector{Num}[[x0, y0], [x0 + 0.08333333333333333((1.5x0) + (1.5(x0 + (0.5((1.5(x0 + (0.25((1.5(x0 + (0.25((1.5x0) - (x0*y0))))) - ((x0 + (0.25((1.5x0) - (x0*y0))))*(y0 + (0.25((x0*y0) - (3y0))))))))) - ((x0 + (0.25((1.5(x0 + (0.25((1.5x0) - (x0*y0))))) - ((x0 + (0.25((1.5x0) - (x0*y0))))*(y0 + (0.25((x0*y0) - (3y0))))))))*(y0 + (0.25(((x0 + (0.25((1.5x0) - (x0*y0))))*(y0 + (0.25((x0*y0) - (3y0))))) - (3(y0 + (0.25((x0*y0) - (3y0))))))))))))) + (2((1.5(x0 + (0.25((1.5x0) - (x0*y0))))) + (1.5(x0 + (0.25((1.5(x0 + (0.25((1.5x0) - (x0*y0))))) - ((x0 + (0.25((1.5x0) - (x0*y0))))*(y0 + (0.25((x0*y0) - (3y0))))))))) - ((x0 + (0.25((1.5x0) - (x0*y0))))*(y0 + (0.25((x0*y0) - (3y0))))) - ((x0 + (0.25((1.5(x0 + (0.25((1.5x0) - (x0*y0))))) - ((x0 + (0.25((1.5x0) - (x0*y0))))*(y0 + (0.25((x0*y0) - (3y0))))))))*(y0 + (0.25(((x0 + (0.25((1.5x0) - (x0*y0))))*(y0 + (0.25((x0*y0) - (3y0))))) - (3(y0 + (0.25((x0*y0) - (3y0))))))))))) - (x0*y0) - ((x0 + (0.5(

ChrisRackauckas

using OrdinaryDiffEq
using Plots
using Flux, DiffEqFlux, Optim

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end

ChrisRackauckas

# --------------------------------------------------------------------------------
#   Ewing model
#   translation by: slwu89@berkeleu.edu (July 2020)
# --------------------------------------------------------------------------------

using DifferentialEquations
using Plots
using LabelledArrays
using StaticArrays

ChrisRackauckas

Conservative Package Compiler Results

1. Using a ton of standard pervasive packages only helps 10%
2. Getting DiffEqBase helps 50%
3. Getting OrdinaryDiffEq without extra precompiles gets to 1.3 seconds, vs 0.6 seconds on the second run. That's pretty good!

ChrisRackauckas

DifferentialEquations.jl Runge-Kutta Integrators vs Taylor Integration

This is against a Taylor integration method that utilizes taylor-mode automatic differentiation in a specifically optimized way for ODE integration via its @taylorize macro (for more details, see the documentation). This demonstration is done on the classic Pleiades n-body problem, a non-stiff ODE since the Taylor integration methods are only for non-stiff equations.

ChrisRackauckas

fjac = (var"##MTIIPVar#585", var"##MTKArg#583")->begin
        @inbounds begin
                let (u₁ˏ₁ˏ₁, u₂ˏ₁ˏ₁, u₃ˏ₁ˏ₁, u₄ˏ₁ˏ₁, u₅ˏ₁ˏ₁, u₆ˏ₁ˏ₁, u₇ˏ₁ˏ₁, u₈ˏ₁ˏ₁, u₉ˏ₁ˏ₁, u₁₀ˏ₁ˏ₁, u₁₁ˏ₁ˏ₁, u₁₂ˏ₁ˏ₁, u₁₃ˏ₁ˏ₁, u₁₄ˏ₁ˏ₁, u₁₅ˏ₁ˏ₁, u₁₆ˏ₁ˏ₁, u₁ˏ₂ˏ₁, u₂ˏ₂ˏ₁, u₃ˏ₂ˏ₁, u₄ˏ₂ˏ₁, u₅ˏ₂ˏ₁, u₆ˏ₂ˏ₁, u₇ˏ₂ˏ₁, u₈ˏ₂ˏ₁, u₉ˏ₂ˏ₁, u₁₀ˏ₂ˏ₁, u₁₁ˏ₂ˏ₁, u₁₂ˏ₂ˏ₁, u₁₃ˏ₂ˏ₁, u₁₄ˏ₂ˏ₁, u₁₅ˏ₂ˏ₁, u₁₆ˏ₂ˏ₁, u₁ˏ₃ˏ₁, u₂ˏ₃ˏ₁, u₃ˏ₃ˏ₁, u₄ˏ₃ˏ₁, u₅ˏ₃ˏ₁, u₆ˏ₃ˏ₁, u₇ˏ₃ˏ₁, u₈ˏ₃ˏ₁, u₉ˏ₃ˏ₁, u₁₀ˏ₃ˏ₁, u₁₁ˏ₃ˏ₁, u₁₂ˏ₃ˏ₁, u₁₃ˏ₃ˏ₁, u₁₄ˏ₃ˏ₁, u₁₅ˏ₃ˏ₁, u₁₆ˏ₃ˏ₁, u₁ˏ₄ˏ₁, u₂ˏ₄ˏ₁, u₃ˏ₄ˏ₁, u₄ˏ₄ˏ₁, u₅ˏ₄ˏ₁, u₆ˏ₄ˏ₁, u₇ˏ₄ˏ₁, u₈ˏ₄ˏ₁, u₉ˏ₄ˏ₁, u₁₀ˏ₄ˏ₁, u₁₁ˏ₄ˏ₁, u₁₂ˏ₄ˏ₁, u₁₃ˏ₄ˏ₁, u₁₄ˏ₄ˏ₁, u₁₅ˏ₄ˏ₁, u₁₆ˏ₄ˏ₁, u₁ˏ₅ˏ₁, u₂ˏ₅ˏ₁, u₃ˏ₅ˏ₁, u₄ˏ₅ˏ₁, u₅ˏ₅ˏ₁, u₆ˏ₅ˏ₁, u₇ˏ₅ˏ₁, u₈ˏ₅ˏ₁, u₉ˏ₅ˏ₁, u₁₀ˏ₅ˏ₁, u₁₁ˏ₅ˏ₁, u₁₂ˏ₅ˏ₁, u₁₃ˏ₅ˏ₁, u₁₄ˏ₅ˏ₁, u₁₅ˏ₅ˏ₁, u₁₆ˏ₅ˏ₁, u₁ˏ₆ˏ₁, u₂ˏ₆ˏ₁, u₃ˏ₆ˏ₁, u₄ˏ₆ˏ₁, u₅ˏ₆ˏ₁, u₆ˏ₆ˏ₁, u₇ˏ₆ˏ₁, u₈ˏ₆ˏ₁, u₉ˏ₆ˏ₁, u₁₀ˏ₆ˏ₁, u₁₁ˏ₆ˏ₁, u₁₂ˏ₆ˏ₁, u₁₃ˏ₆ˏ₁, u₁₄ˏ₆ˏ₁, u₁₅ˏ₆ˏ₁, u₁₆ˏ₆ˏ₁, u₁ˏ₇ˏ₁, u₂ˏ₇ˏ₁, u₃ˏ₇ˏ₁, u₄ˏ₇ˏ₁, u₅ˏ₇ˏ₁, u₆ˏ₇ˏ₁, u₇ˏ₇ˏ₁, u₈ˏ₇ˏ₁, u₉ˏ₇ˏ₁, u₁₀ˏ₇ˏ₁, u₁₁ˏ₇ˏ₁, u₁₂ˏ₇ˏ₁, u₁₃ˏ₇ˏ₁, u₁₄ˏ₇ˏ

ChrisRackauckas

function fast_gm(t,u,p,du)
  a,α,Abar,β,D1,D2 = p

  @inbounds for i in 2:N-1, j in 2:N-1
    du[i,j,1] = D1*(u[i-1,j,1] + u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
              a*u[i,j,1]^2/u[i,j,2] + Abar - α*u[i,j,1]
  end

  @inbounds for i in 2:N-1, j in 2:N-1
    du[i,j,2] = D2*(u[i-1,j,2] + u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +

ChrisRackauckas

:((var"##MTIIPVar#414", var"##MTKArg#412")->begin
          @inbounds begin
                  @sync begin
                          let (u₁, u₂, u₃, u₄, u₅, u₆, u₇, u₈, u₉, u₁₀) = (var"##MTKArg#412"[1], var"##MTKArg#412"[2], var"##MTKArg#412"[3], var"##MTKArg#412"[4], var"##MTKArg#412"[5], var"##MTKArg#412"[6], var"##MTKArg#412"[7], var"##MTKArg#412"[8], var"##MTKArg#412"[9], var"##MTKArg#412"[10])
                              begin
                                  Threads.@spawn begin
                                          (var"##MTIIPVar#414").nzval[1] = -2
                                          (var"##MTIIPVar#414").nzval[2] = 1
                                          (var"##MTIIPVar#414").nzval[3] = 1
                                          (var"##MTIIPVar#414").nzval[4] = -2

ChrisRackauckas

multithreadedfjac = (var"##MTIIPVar#581", var"##MTKArg#579")->begin
        @inbounds begin
                let (u₁ˏ₁ˏ₁, u₂ˏ₁ˏ₁, u₃ˏ₁ˏ₁, u₄ˏ₁ˏ₁, u₅ˏ₁ˏ₁, u₆ˏ₁ˏ₁, u₇ˏ₁ˏ₁, u₈ˏ₁ˏ₁, u₉ˏ₁ˏ₁, u₁₀ˏ₁ˏ₁, u₁₁ˏ₁ˏ₁, u₁₂ˏ₁ˏ₁, u₁₃ˏ₁ˏ₁, u₁₄ˏ₁ˏ₁, u₁₅ˏ₁ˏ₁, u₁₆ˏ₁ˏ₁, u₁ˏ₂ˏ₁, u₂ˏ₂ˏ₁, u₃ˏ₂ˏ₁, u₄ˏ₂ˏ₁, u₅ˏ₂ˏ₁, u₆ˏ₂ˏ₁, u₇ˏ₂ˏ₁, u₈ˏ₂ˏ₁, u₉ˏ₂ˏ₁, u₁₀ˏ₂ˏ₁, u₁₁ˏ₂ˏ₁, u₁₂ˏ₂ˏ₁, u₁₃ˏ₂ˏ₁, u₁₄ˏ₂ˏ₁, u₁₅ˏ₂ˏ₁, u₁₆ˏ₂ˏ₁, u₁ˏ₃ˏ₁, u₂ˏ₃ˏ₁, u₃ˏ₃ˏ₁, u₄ˏ₃ˏ₁, u₅ˏ₃ˏ₁, u₆ˏ₃ˏ₁, u₇ˏ₃ˏ₁, u₈ˏ₃ˏ₁, u₉ˏ₃ˏ₁, u₁₀ˏ₃ˏ₁, u₁₁ˏ₃ˏ₁, u₁₂ˏ₃ˏ₁, u₁₃ˏ₃ˏ₁, u₁₄ˏ₃ˏ₁, u₁₅ˏ₃ˏ₁, u₁₆ˏ₃ˏ₁, u₁ˏ₄ˏ₁, u₂ˏ₄ˏ₁, u₃ˏ₄ˏ₁, u₄ˏ₄ˏ₁, u₅ˏ₄ˏ₁, u₆ˏ₄ˏ₁, u₇ˏ₄ˏ₁, u₈ˏ₄ˏ₁, u₉ˏ₄ˏ₁, u₁₀ˏ₄ˏ₁, u₁₁ˏ₄ˏ₁, u₁₂ˏ₄ˏ₁, u₁₃ˏ₄ˏ₁, u₁₄ˏ₄ˏ₁, u₁₅ˏ₄ˏ₁, u₁₆ˏ₄ˏ₁, u₁ˏ₅ˏ₁, u₂ˏ₅ˏ₁, u₃ˏ₅ˏ₁, u₄ˏ₅ˏ₁, u₅ˏ₅ˏ₁, u₆ˏ₅ˏ₁, u₇ˏ₅ˏ₁, u₈ˏ₅ˏ₁, u₉ˏ₅ˏ₁, u₁₀ˏ₅ˏ₁, u₁₁ˏ₅ˏ₁, u₁₂ˏ₅ˏ₁, u₁₃ˏ₅ˏ₁, u₁₄ˏ₅ˏ₁, u₁₅ˏ₅ˏ₁, u₁₆ˏ₅ˏ₁, u₁ˏ₆ˏ₁, u₂ˏ₆ˏ₁, u₃ˏ₆ˏ₁, u₄ˏ₆ˏ₁, u₅ˏ₆ˏ₁, u₆ˏ₆ˏ₁, u₇ˏ₆ˏ₁, u₈ˏ₆ˏ₁, u₉ˏ₆ˏ₁, u₁₀ˏ₆ˏ₁, u₁₁ˏ₆ˏ₁, u₁₂ˏ₆ˏ₁, u₁₃ˏ₆ˏ₁, u₁₄ˏ₆ˏ₁, u₁₅ˏ₆ˏ₁, u₁₆ˏ₆ˏ₁, u₁ˏ₇ˏ₁, u₂ˏ₇ˏ₁, u₃ˏ₇ˏ₁, u₄ˏ₇ˏ₁, u₅ˏ₇ˏ₁, u₆ˏ₇ˏ₁, u₇ˏ₇ˏ₁, u₈ˏ₇ˏ₁, u₉ˏ₇ˏ₁, u₁₀ˏ₇ˏ₁, u₁₁ˏ₇ˏ₁, u₁₂ˏ₇ˏ₁, u₁