Problem with DifferentialEquations & StaticArrays

Specific Domains Numerics
1

Hello,

I would like to solve the Lorenz-96 problem with static arrays. The state size is 40. I have followed the tutorial Optimizing DiffEq Code but I have the following error message:

using StaticArrays
using DifferentialEquations

function lorenz96ode!(du,u,p,t)
    F = 8.0
    n = size(u,1)
    du[1] = (u[2]-u[end-1])*u[end] - u[1] + F
    du[2] = (u[3]-u[end])*u[1] - u[2] + F
    du[end] = (u[1] - u[end-2])*u[end-1] - u[end] + F

    @inbounds for i=3:n-1
        du[i] = (u[i+1] - u[i-2])*u[i-1] - u[i] + F
    end
    SVector{40,Float64}(du)
end
u0 = SVector{40, Float64}(randn(40))
tspan = (0.0,100.0)
prob = ODEProblem(lorenz96ode!,u0,tspan)
 sol = solve(prob, Tsit5())

MethodError: no method matching OrdinaryDiffEq.Tsit5Cache(::SArray{Tuple{40},Float64,1,40}, ::SArray{Tuple{40},Float64,1,40}, ::SArray{Tuple{40},Float64,1,40}, ::SArray{Tuple{40},Float64,1,40}, ::SArray{Tuple{40},Float64,1,40}, ::SArray{Tuple{40},Float64,1,40}, ::SArray{Tuple{40},Float64,1,40}, ::SArray{Tuple{40},Float64,1,40}, ::SArray{Tuple{40},Float64,1,40}, ::MArray{Tuple{40},Float64,1,40}, ::MArray{Tuple{40},Float64,1,40}, ::MArray{Tuple{40},Float64,1,40}, ::OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64})
Closest candidates are:
OrdinaryDiffEq.Tsit5Cache(::#217#uType, ::#217#uType, ::#218#rateType, ::#218#rateType, ::#218#rateType, ::#218#rateType, ::#218#rateType, ::#218#rateType, ::#218#rateType, !Matched::#217#uType, !Matched::#217#uType, ::#219#uNoUnitsType, ::#220#TabType) where {#217#uType, #218#rateType, #219#uNoUnitsType, #220#TabType} at /home/mat/.julia/packages/OrdinaryDiffEq/5igHh/src/caches/low_order_rk_caches.jl:272

2

You can’t mutate static arrays. This should just be done with standard arrays.

1 Like
3

If I use a standard array for the initial condition, the code runs. Is there any advantage to use StaticArrays in the rhs then?

using StaticArrays
using DifferentialEquations

function lorenz96ode!(du,u,p,t)
    F = 8.0
    n = size(u,1)
    du[1] = (u[2]-u[end-1])*u[end] - u[1] + F
    du[2] = (u[3]-u[end])*u[1] - u[2] + F
    du[end] = (u[1] - u[end-2])*u[end-1] - u[end] + F

    @inbounds for i=3:n-1
        du[i] = (u[i+1] - u[i-2])*u[i-1] - u[i] + F
    end
    SVector{40,Float64}(du)
end
u0 = randn(40)
tspan = (0.0,100.0)
prob = ODEProblem(lorenz96ode!,u0,tspan)
 sol = solve(prob, Tsit5())
4

Not at that size. I have found that static arrays will begin to loose their performance benefit at around 20ish ODEs, so lower than it used to be, and I am not sure what happened in their dispatches that caused that. Definitely not 100 anymore and we should update the docs.

4 Likes
5

Thank you for your explanation.