Why is the argument of vectorized in-place function is marked as not used?

New to Julia
, , ,
1

Following the tutorial, I want to use OrdinaryDiffEq to solve an equation:

using OrdinaryDiffEq

function A(t, p)
    ω, ε = p
    return [[0.0 1.0]; [-(ω + ε * cos(t)) 0.0]]
end

function mathieu!(du, u, p, t)
    du = A(t, p) * u
end

u0 = [1.0, 1.0]
p = [1.0, 1.0]
tspan = (0.0, 20 * π)
prob = ODEProblem(mathieu!, u0, tspan, p)
sol = solve(prob, Tsit5())

Such definition doesn’t produce the correct result, and the linter in VS Code warns about du:

An argument is included in a function signature but not used within its body.Julia(UnusedFunctionArgument)

However, if I define the system like in the tutorial

function mathieu!(du, u, p, t)
    du[1] = u[2]
    du[2] = # ...
end

everything works fine.

Why does my definition not work? If the u vector had 1000 of components, would I have to write every one of them by hand?

2

You need to use .= instead of = to reuse the memory in du—using = will just create a new local variable called du.

3

That is, assignment operators are different for scalars and vectors, and only in in-place functions?

4

Not just in in-place functions or for vectors. In general, x = … says “attach the label x to ….” If x is an object that already exists, this is different from changing/mutating the object x!

What .= says instead is “do x[i] = … for every i. Unlike x = …, x[i] = ... really does change x!

3 Likes
5

While du .= A(t, p) * u would work here, use mul!(du, A(t, p), u) to avoid allocation in that operation altogether. You will need to using LinearAlgebra to access that function.

You should be able to write this more simply (and likely with fewer allocations) as return [0 1; -(ω + ε * cos(t)) 0]. Or, if you don’t mind it being a little less “clean” you could write out the operation manually as

function mathieu!(du, u, p, t)
    du[1] = u[2]
    du[2] = -(p[1] + p[2] * cos(t)) * u[1]
    return du
end

Since you’re working with tiny vectors and matrices, you might eventually check out StaticArrays.jl. They can be considerably faster at small sizes. Although it may take a little work to adapt this to it, so only try that once you have this working and if you think a little extra performance is useful.

1 Like
6

I think you’re misunderstanding what assignment does. See the manual section on “assignment vs mutation”.

5 Likes