Using du in system of ODEs

Specific Domains Numerics
, ,
1

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
2

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
3

Are you trying to define a differential-algebraic equation?

4

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())
plot_1
plot_2