System of Two Differential equations

I am trying to implement a system of two ODEs. The following works:

function diffusion!(du, u, p, t)
α, β = p
c, ce = u
du[1] = -αc
du[2] = -β
ce
du = [dc, dce]
end

u0 = [3.0, 7.0]
tspan = (0., 2.)
p = [1.0, 1.]

prob = ODEProblem(diffusion!, u0, tspan, p)
sol = solve(prob, Tsit5())

But the does not: notice the difference in the diffusion! function.

function diffusion!(du, u, p, t)
α, β = p
c, ce = u
dc = -αc
dce = -β
ce
du = [dc, dce]
end

The second version leads to a solution always equal to the initial conditions, while the first case leads to the proper solution. I do not understand why the line

du = [dc, dce]

does not do what I expect. Thanks. Gordon.

Instead updating the contents of du, with the line du = [dc, dce] you are creating a new vector and setting du to point (binding) to that new vector. The underlying DiffEq code is still using the old one and so does not see any change.

An alternative syntax that would work is broadcasting du .= [dc, dce] since this effectively expands out into the syntax in your first example where you mutate the underlying vector.

2 Likes

Is here a typo here [first example]… just before end, you define du = [dc,dce], but dc and dce have not been defined…I assume that your example should be:

function diffusion!(du, u, p, t)
α, β = p
c, ce = u
du[1] = -α c
du[2] = -β ce
#du = [dc, dce]
end

?

Thank you!! That is what I was missing. Coming from the world of C++ and Python, it will take me a while to get used to this approach using broadcasting.

Yes, you are correct. I made an error. I do have the answer I need though, using broadcasting. I appreciate the time you took to help out. Gordon.