# Using du in system of ODEs

**URL:** <https://discourse.julialang.org/t/using-du-in-system-of-odes/114590>\
**Category:** Numerics\
**Tags:** question, ode, differentialequation\
**Created:** [May 22, 2024, 7:17pm UTC](https://discourse.julialang.org/t/using-du-in-system-of-odes/114590 "2024-05-22T19:17:13Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![cgreysongaito](https://avatars.discourse-cdn.com/v4/letter/c/d9b06d/32.png) [@cgreysongaito](https://discourse.julialang.org/u/cgreysongaito)\
**Post date:** [May 22, 2024, 7:17pm UTC](https://discourse.julialang.org/t/using-du-in-system-of-odes/114590/1 "2024-05-22T19:17:14Z")

</div>

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?

```julia
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

```

---

<div class="post-metadata">

**Author:** ![JonasWickman](https://avatars.discourse-cdn.com/v4/letter/j/9de0a6/32.png) [@JonasWickman](https://discourse.julialang.org/u/JonasWickman)\
**Post date:** [May 22, 2024, 7:39pm UTC](https://discourse.julialang.org/t/using-du-in-system-of-odes/114590/2 "2024-05-22T19:39:25Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 23, 2024, 8:11pm UTC](https://discourse.julialang.org/t/using-du-in-system-of-odes/114590/3 "2024-05-23T20:11:57Z")

</div>

Are you trying to define a differential-algebraic equation?

> **[Differential Algebraic Equations · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/dae_example/)**
>
> Documentation for DifferentialEquations.jl.

---

<div class="post-metadata">

**Author:** ![Pugolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pugolas/32/209242_2.png) [@Pugolas](https://discourse.julialang.org/u/Pugolas)\
**Post date:** [May 29, 2024, 10:42am UTC](https://discourse.julialang.org/t/using-du-in-system-of-odes/114590/4 "2024-05-29T10:42:20Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/c/e/ceb668fa33ec247150e362b5f770c1cb142cb7b9.png)  
 ![plot_2](https://global.discourse-cdn.com/julialang/original/3X/1/0/10208b7f94c5b1f98f005d044a894304b60d0547.png)