# Using du in system of ODEs

I am numerically solving a system of odes. One of the equations uses two other odes within the system of odes. However, I get different results depending on the order of the equations in the function. For MWEtest1!, ωₘ, ωₑ, λ all land on a periodic orbit. For MWEtest2!, ωₘ, ωₑ, λ continue to expand (as well as cycle). Why does the difference in the order of the equations change the outcome?

``````using DifferentialEquations
using Parameters
using PyPlot

@with_kw mutable struct MWESystemPar
νₘ = 3.0
νₑ = 3.0
αₘ = 0.025
αₑ = 0.025
δₘ = 0.01
δₑ = 0.01
β = 0.02
ϕ₀ = 0.95
ϕ₁ = 1.0
aem = 0.1
amm = 0.1
end

function MWEtest1!(du, u, p, t, )
@unpack νₘ, νₑ, αₘ, αₑ, δₘ, δₑ, β, ϕ₀, ϕ₁, aem, amm  = p
Aₘ, Aₑ, ωₘ, ωₑ, Kₘ, Kₑ, λ = u
du[1] = Aₘ * αₘ
du[2] = Aₑ * αₑ
du[3] = ωₘ * (-ϕ₀ + ϕ₁ * λ - αₘ)
du[4] = ωₑ * (-ϕ₀ + ϕ₁ * λ - αₘ)
du[5] = Kₘ * (((1 - ωₘ) / νₘ) - δₘ - aem * (1 / (νₘ * (1 - amm))) )
du[6] = Kₑ * (aem * (Kₘ / Kₑ) * (1 / (νₘ * (1 - amm))) + ((1 - ωₑ) / νₑ) - δₑ )
du[7] = λ * ((1 / ((Kₘ / (νₘ * Aₘ)) + (Kₑ / (νₑ * Aₑ))) )* ((Kₘ / (νₘ * Aₘ))*((du[5]/Kₘ) - αₘ ) + (Kₑ / (νₑ * Aₑ))*((du[6] / Kₑ) - αₑ ) ) - β)
return
end

function MWEtest2!(du, u, p, t, )
@unpack νₘ, νₑ, αₘ, αₑ, δₘ, δₑ, β, ϕ₀, ϕ₁, aem, amm  = p
Aₘ, Aₑ, ωₘ, ωₑ, λ, Kₘ, Kₑ = u
du[1] = Aₘ * αₘ
du[2] = Aₑ * αₑ
du[3] = ωₘ * (-ϕ₀ + ϕ₁ * λ - αₘ)
du[4] = ωₑ * (-ϕ₀ + ϕ₁ * λ - αₘ)
du[5] = λ * ((1 / ((Kₘ / (νₘ * Aₘ)) + (Kₑ / (νₑ * Aₑ))) )* ((Kₘ / (νₘ * Aₘ))*((du[6]/Kₘ) - αₘ ) + (Kₑ / (νₑ * Aₑ))*((du[7] / Kₑ) - αₑ ) ) - β)
du[6] = Kₘ * (((1 - ωₘ) / νₘ) - δₘ - aem * (1 / (νₘ * (1 - amm))) )
du[7] = Kₑ * (aem * (Kₘ / Kₑ) * (1 / (νₘ * (1 - amm))) + ((1 - ωₑ) / νₑ) - δₑ )
return
end

let
u0 = [1.0, 1.0, 0.8,0.8,100, 100, 0.8]
t_span =(0.0, 500.0)
prob = ODEProblem(MWEtest1!, u0, t_span, MWESystemPar())
sol = solve(prob, reltol = 1e-8)
test = figure()
subplot(4,2,1)
plot(sol.t, sol[1,:])
xlabel("Time")
ylabel("Aₘ")
subplot(4,2,2)
plot(sol.t, sol[2,:])
xlabel("Time")
ylabel("Aₑ")
subplot(4,2,3)
plot(sol.t, sol[3,:])
xlabel("Time")
ylabel("ωₘ")
subplot(4,2,4)
plot(sol.t, sol[4,:])
xlabel("Time")
ylabel("ωₑ")
subplot(4,2,5)
plot(sol.t, sol[5,:])
xlabel("Time")
ylabel("Kₘ")
subplot(4,2,6)
plot(sol.t, sol[6,:])
xlabel("Time")
ylabel("Kₑ")
subplot(4,2,7)
plot(sol.t, sol[7,:])
xlabel("Time")
ylabel("λ")
tight_layout()
return test
end

let
u0 = [1.0, 1.0, 0.8,0.8, 0.8, 100, 100]
t_span =(0.0, 500.0)
prob = ODEProblem(MWEtest2!, u0, t_span, MWESystemPar())
sol = solve(prob, reltol = 1e-8)
test = figure()
subplot(4,2,1)
plot(sol.t, sol[1,:])
ylabel("Aₘ")
subplot(4,2,2)
plot(sol.t, sol[2,:])
ylabel("Aₑ")
subplot(4,2,3)
plot(sol.t, sol[3,:])
ylabel("ωₘ")
subplot(4,2,4)
plot(sol.t, sol[4,:])
ylabel("ωₑ")
subplot(4,2,5)
plot(sol.t, sol[5,:])
ylabel("λ")
subplot(4,2,6)
plot(sol.t, sol[6,:])
ylabel("Kₘ")
subplot(4,2,7)
plot(sol.t, sol[7,:])
ylabel("Kₑ")
tight_layout()
return test
end
``````

For `MWEtest2!` you are using `du[6]` and `du[7]` before you have assigned any value to them, and they could contain basically anything depending on exactly how the specific solver may or may not use `du` for internal storage between steps.

1 Like

Are you trying to define a differential-algebraic equation?

I noticed that you are using the default solver, which may be related to this? Your image is noticeably oscillating, and when specified as the rigid solver QNDF(), it seems that it will not change due to changes in order.
`sol = solve(prob, reltol = 1e-8, QNDF())`