using DifferentialEquations: u is not updating

Question

I believe there is a bug in this code. For the sake of brevity though I will just write the function which defines the ODE

function clones(du,u,p,t)

    (Nmut,f) = p

    # average fitness
    phi = sum(f.*u)

    # constructing mutation kernel
    eps = 0.01
    Q = Qmatrix(Nmut,eps)

    # defining differential equations
    Nclones=2^Nmut;
    du = zeros(Nclones)

    ufQ = transpose(transpose(u.*f)*Q)
    du = ufQ .- phi*u

end

If the whole code is needed I can provide it, but it is messy and I'm not sure how to create a minimal example. I tried this when Nmut = 2 so I can compare to a hard-coded version. The output of du at the first time steps are identical. But this version does not seem to ever update u it stays at the u0 prescribed.

Does anyone have an idea why this might be the case? I can also provide the full script, but wanted to avoid that if someone could just see why u would not update.

EDIT:

maxdim=4;
for i in 1:maxdim
    du[i] = 0.0;
    for j in 1:maxdim
        du[i] += u[j].*w[j].*Q[j,i] 
    end
    du[i] -= u[i].*phi
end

The du is updated correctly if we use this version. Why would that be the case?

Answer 1

You are using the in-place form. With this, you have to update the values of du. In your script you are using

du = ufQ .- phi*u

That is replacing the name du with a new array, but not mutating the values in the original du array. A quick fix would be to use the mutating equals:

du .= ufQ .- phi*u

Notice that this is .=.

To understand what this means in a more example based format, think about this. We have an array:

a = [1,2,3,4]

Now we point a new variable to that same array

a2 = a

When we change a value of a we see that reflected in a2 since they point to the same memory:

a[1] = 5
println(a2) # [5,2,3,4]

But now if we replace a we notice nothing happens to a2 since they no longer refer to the same array

a = [1,2,3,4]
println(a2) # [5,2,3,4]

Packages like DifferentialEquations.jl utilize mutating forms so that way users can get rid of array allocations by repeatedly changing the values in cached arrays. As a consequence these means that you should be updating the values of du and not replacing its pointer.

If you feel more comfortable not using mutation, you can use the function syntax f(u,p,t), although there will be a performance hit for doing so if the state variables are (non-static) arrays.