Empty container reduce error: diff eq

Specific Domains Modelling & Simulations
1

Hi all,

I’m trying to use the DifferentialEquations.jl package to perform explicit integration of a mass-spring-damper system.

Here are the function I’m solving, and how I passed it into the SecondOrderODEProblem() function.

func = inv(M) * (f - (C *v + K* d))


SecondOrderODEProblem(func, v, d, 0:tsim)

This gives me the following error:

ERROR: LoadError: ArgumentError: reducing over an empty collection is not allowed

The stack trace didn’t give much insight into which collection was empty, but I’ve just checked, and all the containers in the function aren’t empty, and my du0/u0 aren’t empty.

Here are the values for each variable in there:

M = 96x96 Diagonal{Float64, Vector{Float64}}
f = 96x1 Matrix{Float64} (every 3rd element is non-zero)
C = 96x96 Diagonal{Float64, Vector{Float64}}
v = 96x1 Matrix{Float64}
K = 96x96 Matrix{Float64}
d = 96x1 Matrix{Float64}

so that

func = 96x1 Matrix{Float64}

So, I should be passing into the SecondOrderODEProblem() function:

96x1 matrix where every 3rd element is non zero
96x1 matrix of zeros
96x1 matrix of zeros

I’m not sure exactly what the problem is here, any help is greatly appreciated!

2

func is not a function. You might want to start by copying a tutorial demo:

function HH_acceleration(dv,v,u,p,t)
    x,y  = u
    dx,dy = dv
    dv[1] = -x - 2x*y
    dv[2] = y^2 - y -x^2
end
initial_positions = [0.0,0.1]
initial_velocities = [0.5,0.0]
prob = SecondOrderODEProblem(HH_acceleration,initial_velocities,initial_positions,tspan)
sol2 = solve(prob, KahanLi8(), dt=1/10);

etc.

3

Thanks yeah, silly mistake.

My code now looks like:

func(M, f, C, v, K, d) = inv(M) * (f - (C*v + K*d))
prob = SecondOrderODEProblem(func, v, d, (0.0,10.0))
sol = solve(prob, Tsit5(), dt = 1.0/10.0)

Which I arrived at after looking at the documented constructor:

SecondOrderODEProblem{isinplace}(f,du0,u0,tspan,callback=CallbackSet())

However, now it is throwing a long winded error:

ERROR: LoadError: MethodError: no method matching (::var"#func#84")(::Matrix{Float64}, ::Matrix{Float64}, ::Matrix{Float64}, ::SciMLBase.NullParameters, ::Float64)
Closest candidates are:
  (::var"#func#84")(::Any, ::Any, ::Any, ::Any, ::Any, ::Any) at c:\Users\grw40\code\Summer Research\MSD_matlab\msd_matlab\julia_port\ExplicitShellMSD.jl:86

I suppose there’s something about my function that doesn’t jive well.
I’m familiar with the syntax for the 2nd derivative function you’ve posted, but I thought the solver might do that for me.

I suppose I’ll have to explicitly define which each component acceleration is for every node within my acceleration functions? Are for loops the best way of going about this (in terms of performance)?

Thanks for the help!

4

That’s not one of the allowed function signatures. Did you mean func(v,u,p,t), with some of those things being parameters in p?

5

I don’t understand what you want to do, but I have made a SecondOrderODEProblem example that actually works.

Jupyter notebook: https://github.com/genkuroki/public/blob/main/0012/SecondOrderODEProblem%20example.ipynb

using DifferentialEquations
using LinearAlgebra
using Plots

function plot_2ndorder(sol; l1=:bottomright, l2=:bottomright, size=(720, 300), kwargs...)
    t = range(sol.prob.tspan...; length=400)
    n = Base.size(sol, 1) ÷ 2
    v = vcat((t -> sol(t)[1:n]').(t)...)
    d = vcat((t -> sol(t)[n+1:end]').(t)...)
    P = plot(t, d; label=permutedims(["\$d_{$i}\$" for i in 1:n]), legend=l1)
    Q = plot(t, v; label=permutedims(["\$v_{$i}\$" for i in 1:n]), legend=l2)
    plot(P, Q; size, kwargs...)
end

m = Float64[1, 2, 3, 4]
f = Float64[-10, -10, -10, -10]
c = Float64[1, 1, 1, 1]
K = Float64[
     2 -1  0  0
    -1  2 -1  0
     0  1  2 -1
     0  0 -1  2
]
p = (m, f, c, K)
v0 = Float64[0, 0, 0, 0]
d0 = Float64[3, 1, -1, -3]
tspan = (0.0, 10.0)

function g!(dv, v, d, p, t)
    m, f, c, K = p
    # The following is equivalent to dv .= M \ (f - (C * v + K * d)),
    # where K = Matrix, C = Diagonal(c), M = Diagonal(m)
    mul!(dv, K, d)
    @. dv = m \ (f - (c * v + dv))
    return
end

prob = SecondOrderODEProblem(g!, v0, d0, tspan, p)
sol = solve(prob)
plot_2ndorder(sol)

image
image1450×607 117 KB

6

Wow, idk why my brain can’t process the documentation. Thanks!

Incidentally, do you know of a good way to generate and plot a 3D planar mesh (think cloth sim, starts out as flat sheet and deforms in 3D)?

I tried Triangulate.jl, but could only find support for plotting in 2D. Tried Tetgen too, but it was pretty buggy (missing connections randomly), and was pretty hard to find good documentation on it - I don’t remember the C++ docs being too helpful, I think I remember the CLI switches being a source of frustration.

Anywho, thanks a lot, and when tf do you sleep? Have a good one.

P.S, you’re Stochastic Lifestyle guy right? Great blog, thanks for doin what you do!

1 Like
7

Hey, thanks for this, really helpful. Idk why my brain couldn’t process the docs, it’s fairly clear!

I’m ultimately trying to rewrite a cloth sim I’ve made into Julia, currently testing out implicit/explicit methods for performance and learnkng about Julia syntax. :slight_smile: Usually just handwrite the integration, but I’m trying to learn more about auto-diff/adaptive time stepping, and Julia is nice in that you get to see the source code for some world-class solvers for free!

I think, with some adjustment to account for n-dof systems, your example will be very helpful thanks!

8

Open an issue. Maybe we need another second order ODE tutorial.

9

Thanks, will do. To be clear, should I open an issue on the DifferentialEquations.jl github, or are you meaning there’s somewhere on the Julia docs site I should look to do this?

10

Yes, on DifferentialEquations.jl

11

Cool, thanks for your patience, have a good one.