# How to Update Non-Time-Dependent States in Mixed ODE Problems?

**URL:** <https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594>\
**Category:** General Usage\
**Tags:** question, differentialequation\
**Created:** [October 22, 2024, 2:13pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594 "2024-10-22T14:13:11Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [October 22, 2024, 2:13pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/1 "2024-10-22T14:13:11Z")

</div>

Hello, I’d like to ask for your help. In an ODE problem, after I complete the calculations for a given time step ( t ), I need to update a certain state. This state is not time-dependent but is recalculated based on another variable ( u ), and then used for subsequent calculations.

I have a mixed ODE problem described as follows:

```julia
# where s is in m^3 and q_in, q_out are in m^3/s
ds/dt = q_in - q_out

dq_in = new_q_in - q_in 

# where new_q_in is calculated by a function func(q_out), and q_gen is a forcing data.
# We can assume that the function is: func = (q_out, q_gen) -> q_out + q_gen
new_q_in = func(q_out, q_gen)

```

As you can see, these two balance equations are different: one updates the state ( s ) (in m³), and the other updates the flow rate ( q\_{\text{in}} ) (in m³/s, which is time-dependent).

I tried to construct an ODE problem based on this formula. The code is as follows:

```julia
using DataInterpolations
using OrdinaryDiffEq

# q_gen
input_arr = [3, 4, 4, 5, 6, 10, 23, 24, 38, 40, 59, 63, 49, 32, 22, 12, 9, 4]
input_itp = LinearInterpolation(input_arr[1,:], collect(1:length(input_arr)))

# ode function
function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river + q_in) / (p[1] + 1)
    du[1] = q_in - q_out
    du[2] = q_out + input_itp(t) - q_in
end

# build problem and solve
prob = DiscreteProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
sol = solve(prob, FunctionMap())

```

In the code, I built a `DiscreteProblem` to solve this problem. However, mathematically, this formula should be continuous, but when constructing an `ODEProblem`, the results are different:

```julia
prob2 = ODEProblem(discharge2!, [0.1, 3], (1, length(input_arr)), [0.3])
sol2 = solve(prob2, Tsit5(), saveat=1.0)

```

The reason lies in the equation `du[2] = q_out + input_itp(t) - q_in`, which doesn’t hold in the `ODEProblem`.

Therefore, I attempted to replace ( u[2] ) by using a `DiscreteCallback`:

```julia
function condition(u, t, integrator)
    return true
end

function affect!(integrator)
    q_in = integrator.u[2]
    q_out = (integrator.uprev[1] + q_in) / (p[1] + 1)
    integrator.u[2] = q_out + input_itp(integrator.t) - q_in
end

cb = DiscreteCallback(condition, affect!)
prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [p])
sol = solve(prob, Tsit5(), dt=1.0, callback=cb)
sol_u = Array(sol)

```

I am not sure if this approach is correct, as the results are still different from earlier (perhaps there’s a mistake in one of the steps).

---

<div class="post-metadata">

**Author:** ![janneshb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/janneshb/32/36035_2.png) [@janneshb](https://discourse.julialang.org/u/janneshb)\
**Post date:** [October 22, 2024, 2:43pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/3 "2024-10-22T14:43:30Z")

</div>

Hi,  
sorry for my earlier answer, I didn’t get your problem at first and took a closer look now.

Why are you trying to implement this as a `DiscreteProblem`?  
Nothing about your problem seems to be discrete…  
In a discrete problem, your function should define a discrete state evolution. I.e. u\_{k+1} = f(u\_n, p, t\_{k+1}).  
This is completely different from the ordinary differential equation that you defined, so of course you get very different results from the `DiscreteProblem` vs the `ODEProblem`.  
Your problem seems to be well defined and I don’t see any issue with taking your first code block and replacing the `DiscreteProblem` part with

```julia
u0 = [0.1, 3]
prob = ODEProblem(discharge!, u0, (1, length(input_arr)), [0.3])
sol = solve(prob, Tsit5())

```

What you are trying to achieve with the discrete callback? As far as I can tell, you don’t need any discrete dynamics here…

---

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [October 22, 2024, 3:09pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/4 "2024-10-22T15:09:39Z")

</div>

Thank you for your answer.

Indeed, this formula is not a discrete problem, but the issue is that u[2], which is q\_in, is not updated like a common state in the ODE.

In the `discharge!` function,

```julia

du[2] = q_out + input_itp(t) - q_in

```

This calculation formula is theoretically problematic. The unit of `q_in` is m3/s, so the formula should be:

```julia

dq_in/dt = (new_q_in - q_in) /dt

```

But in the code, `dt` is ignored, which is why the calculation results differ between DiscreteProblem (dt=1) and ODEProblem(dt!=1).

I also believe that this calculation formula doesn’t need discrete dynamics, but when using DifferentialEquations.jl, I don’t know how to update the state variable q\_in. So I considered using a discrete callback, which seems to support updating q\_in.

---

<div class="post-metadata">

**Author:** ![janneshb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/janneshb/32/36035_2.png) [@janneshb](https://discourse.julialang.org/u/janneshb)\
**Post date:** [October 22, 2024, 6:44pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/5 "2024-10-22T18:44:17Z")

</div>

I honestly don’t really understand your problem. What’s a theoretical sound formulation of this problem, not considering any implementation difficulties?

As far as I can tell, you have an ODE

```julia
ds = q_in - q_out

```

where `q_in` is computed as

```julia
q_in = q_out + q_gen

```

and `q_out` is computed as

```julia
q_out = (s + q_in) / (p[1] + 1)

```

`q_gen` is an input to your system that you compute via linear interpolation of a vector.

Why do you need any discrete dynamics here?

---

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [October 23, 2024, 2:09am UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/6 "2024-10-23T02:09:23Z")

</div>

I apologize for not expressing my problem well. The ODE problem you described is somewhat different from my idea:

First, I have an ODE problem:

bash

```julia
ds/dt = q_in - q_out # the unit of s is m³, q_in and q_out are m³/s

```

Then calculate q\_out:

```julia

q_out = s / (p+1) # here I replace s+q_in with s, because the units of these two variables are different

```

Finally, I need to obtain the next q\_in for calculation based on q\_out combined with q\_gen, rather than the q\_in for this calculation, which is provided by the state u:

`q_in(t+dt) = func(q_out(t), q_gen(t))`

After obtaining the q\_in for the next calculation, I set u[2] = q\_in(t+dt) to use it for the next calculation.

I constructed a discrete loop calculation to help you understand this problem:

```julia

input_arr = [3 4 4 5 6 10 23 24 38 40 59 63 49 32 22 12 9 4]
s_river = zeros(length(input_arr))
q_in = zeros(length(input_arr))
q_out = zeros(length(input_arr))
s_river0, q_in0 = 0.1, 3
p = [0.3]

for i in eachindex(input_arr)
    q_out[i] = (s_river0) / (p[1] + 1)
    s_river[i] = s_river0 + q_in0 - q_out[i]
    q_in[i] = q_out[i] + input_arr[i]
    s_river0, q_in0 = s_river[i], q_in[i]
    @info [q_out[i] q_in[i] s_river[i]]
end

```

Then I solved this problem using both ODEProblem and DiscreteProblem, and the results were consistent with the discrete loop calculation, thus confirming your statement:

```julia

function discharge2!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river) / (p[1] + 1)
    du[1] = q_in - q_out
    du[2] = q_out + input_itp(t) - q_in
end

prob = ODEProblem(discharge2!, [0.1, 3], (1, length(input_arr)), [0.3])
sol3 = solve(prob, Tsit5(), saveat=1.0)

function discharge2(u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river) / (p[1] + 1)
    [s_river + q_in - q_out, q_out + input_itp(t)]
end

prob = DiscreteProblem(discharge2, [0.1, 3], (1, length(input_arr)), [0.3])
sol4 = solve(prob, FunctionMap(), saveat=1.0)

```

The results are plotted as follows:

 ![test_routeflux](https://global.discourse-cdn.com/julialang/original/3X/6/e/6e84710ca22ddbb6ca333e17f39ff1b1c06ac939.png)  
My problem seems to be resolved, but I’m not quite clear why `du[2] = q_out + input_itp(t) - q_in` can update without considering units. Does du refer to du/dt or just du?

---

<div class="post-metadata">

**Author:** ![janneshb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/janneshb/32/36035_2.png) [@janneshb](https://discourse.julialang.org/u/janneshb)\
**Post date:** [October 23, 2024, 1:26pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/7 "2024-10-23T13:26:20Z")

</div>

I don’t really understand the issue you have with units… If your problem is formulated well, all units work out automatically.

But this looks better now:

> [@chooron](#):
>
> q\_in(t+dt) = func(q\_out(t), q\_gen(t))

However, this

> [@chooron](#):
>
> `du[2] = q_out + input_itp(t) - q_in`

doesn’t make a lot of sense, in my opinion. `du` is a time derivative (in math you would write `du/dt`), so `du[2]` would denote the _rate of change_ of `q_in`, not the next `q_in(t+dt)`. For that you indeed need a discrete problem.

Discrete and continuous problems are completely different. You cannot use the same `discharge` function for both of them.  
Also, `u(t + dt) = s_river + q_in - q_out` doesn’t seem correct if you’re are assuming sample times other than `T = 1`.  
So, I think the fact that your results make sense are a bit of a coincidence.

For this problem:

> [@chooron](#):
>
> `ds/dt = q_in - q_out`

> [@chooron](#):
>
> q\_in(t+dt) = func(q\_out(t), q\_gen(t))

You do need some hybrid dynamics. I would propose the following:

```julia
using DataInterpolations
using OrdinaryDiffEq

input_arr = [3, 4, 4, 5, 6, 10, 23, 24, 38, 40, 59, 63, 49, 32, 22, 12, 9, 4]

function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river + q_in) / (p[1] + 1)
           
    du[1] = q_in - q_out
    du[2] = 0 # q_in is a discrete state, so it will not be changed by the ODE
end

function update_q_in!(integrator)
    q_out = (integrator.u[1] + integrator.u[2]) / (integrator.p[1] + 1)
    integrator.u[2] = q_out + input_arr[Int(round(integrator.t))]
end
cb = PeriodicCallback(update_q_in!, 1.0)
prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
sol = solve(prob, Tsit5(), dt=1.0, callback=cb)

```

This will update `q_in` once a second based on the current `q_out` and `q_gen` (from the array).  
I don’t understand your problem on a deeper physical level, so I don’t know if this makes sense, but it is a good way to implement such hybrid systems.

Generally, don’t mix continuous and discrete states. If your state evolution is described by an ODE (`du[2] = ...`) then it probably isn’t also governed by some discrete rule.  
Hope this helps 🙂

---

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [October 23, 2024, 2:11pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/8 "2024-10-23T14:11:08Z")

</div>

> [@janneshb](#):
>
> However, this
> 
> > [@chooron](#):
> >
> > `du[2] = q_out + input_itp(t) - q_in`
> 
> doesn’t make a lot of sense, in my opinion. `du` is a time derivative (in math you would write `du/dt`), so `du[2]` would denote the _rate of change_ of `q_in`, not the next `q_in(t+dt)`. For that you indeed need a discrete problem.

This is exactly what I wanted to express. du[2] represents the rate of change of q\_in, so it cannot be represented by q\_out + input\_itp(t) - q\_in. That’s why I think this might be a discrete problem.

In your code

> [@janneshb](#):
>
> ```julia
> using DataInterpolations
> using OrdinaryDiffEq
> 
> input_arr = [3, 4, 4, 5, 6, 10, 23, 24, 38, 40, 59, 63, 49, 32, 22, 12, 9, 4]
> 
> function discharge!(du, u, p, t)
> s_river, q_in = u[1], u[2]
> q_out = (s_river + q_in) / (p[1] + 1)
>            
> du[1] = q_in - q_out
> du[2] = 0 # q_in is a discrete state, so it will not be changed by the ODE
> end
> 
> function update_q_in!(integrator)
> q_out = (integrator.u[1] + integrator.u[2]) / (integrator.p[1] + 1)
> integrator.u[2] = q_out + input_arr[Int(round(integrator.t))]
> end
> cb = PeriodicCallback(update_q_in!, 1.0)
> prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
> sol = solve(prob, Tsit5(), dt=1.0, callback=cb)
> 
> ```

I’m using DifferentialEquations.jl version 7.14.0, so I’m using DiscreteCallback. The code is as follows:

```julia
function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river) / (p[1] + 1)
    du[1] = q_in - q_out
end

function condition(u, t, integrator)
    return true
end

function affect!(integrator)
    q_out = (integrator.u[1]) / (integrator.p[1] + 1)
    integrator.u[2] = q_out + input_itp(integrator.t)
end

cb = DiscreteCallback(condition, affect!)
prob = ODEProblem(discharge!, [0.1, 3], (timeidx[1], timeidx[end]), [0.3])
sol1 = solve(prob, Tsit5(), dt=1.0, callback=cb)

```

I hope this code should be correct. Thank you for your suggestions！ 😁

---

<div class="post-metadata">

**Author:** ![janneshb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/janneshb/32/36035_2.png) [@janneshb](https://discourse.julialang.org/u/janneshb)\
**Post date:** [October 23, 2024, 2:20pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/9 "2024-10-23T14:20:38Z")

</div>

> [@chooron](#):
>
> ```julia
> function condition(u, t, integrator)
> return true
> end
> 
> ```

This makes no sense, though. If your condition is always true, your discrete time dynamics depend on the step size of the continuous time integrator. That’s why I used a `PeriodicCallback` with sampling time `1.0`. This is crucial. If you solve your `ODEProblem` with a fixed step size `dt=1.0` you are missing out on a lot of nice integration methods provided by `DifferentialEquations.jl`.

---

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [October 23, 2024, 2:25pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/10 "2024-10-23T14:25:01Z")

</div>

I just find that `PeriodicCallback ` is provided in [DiffEqCallbacks.jl](https://docs.sciml.ai/DiffEqCallbacks/stable/), I will try it.

---

<div class="post-metadata">

**Author:** ![janneshb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/janneshb/32/36035_2.png) [@janneshb](https://discourse.julialang.org/u/janneshb)\
**Post date:** [October 23, 2024, 2:30pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/11 "2024-10-23T14:30:57Z")

</div>

You can use `PeriodicCallback` if you use `using DifferentialEquations`.

The problem with your implementation is, that you did not define a sampling time for your discrete callback. It will _always_ be executed. For that to work you need to define a fixed step size `dt = 1.0`, like you did.  
But this is something you don’t usually want for continuous time systems.  
So maybe think about using a periodic callback or adjust your `condition()` function to only return `true` once every `dt` seconds.

---

<div class="post-metadata">

**Author:** ![chooron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chooron/32/208300_2.png) [@chooron](https://discourse.julialang.org/u/chooron)\
**Post date:** [October 23, 2024, 2:56pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/12 "2024-10-23T14:56:17Z")

</div>

I have implemented the method you provided, but it seems to have discrepancies with the calculation results of loop discretization:

```julia
# your method
function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river + q_in) / (p[1] + 1)
           
    du[1] = q_in - q_out
    du[2] = 0 # q_in is a discrete state, so it will not be changed by the ODE
end

function update_q_in!(integrator)
    # there may be uprev that we need to use the previous value of s_river
    q_out = (integrator.uprev[1] + integrator.u[2]) / (integrator.p[1] + 1)
    integrator.u[2] = q_out + input_itp(integrator.t)
end
cb = PeriodicCallback(update_q_in!,1.0)
prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
sol1 = solve(prob, Tsit5(), dt=1.0, callback=cb)

# Use a for loop to calculate
input_arr = [3 4 4 5 6 10 23 24 38 40 59 63 49 32 22 12 9 4]
s_river = zeros(length(input_arr))
q_in = zeros(length(input_arr))
q_out = zeros(length(input_arr))
s_river0, q_in0 = 0.1, 3
p = [0.3]
for i in eachindex(input_arr)
    q_out[i] = (s_river0) / (p[1] + 1)
    s_river[i] = s_river0 + q_in0 - q_out[i]
    q_in[i] = q_out[i] + input_arr[i]
    s_river0, q_in0 = s_river[i], q_in[i]
    @info [q_out[i] q_in[i] s_river[i]]
end

# plot results
plot(sol1)
plot!(s_river, label="s_river, u[1]")
plot!(q_in, label="q_in, u[2]")

```

 ![routeflux](https://global.discourse-cdn.com/julialang/original/3X/6/6/666653e978a1fb635b81479d08eb7a546178835a.png)  
It can be seen that there is a significant discrepancy between these two calculation methods. The problem may lie in the fact that u[2] (q\_in) is only updated at dt intervals, so u[2] appears to be in a step-like pattern.

However, in reality, u[2] should be continuously updated based on the constantly changing q\_out and q\_gen, so the result should be continuous. Perhaps I should consider using other callbacks.

---

<div class="post-metadata">

**Author:** ![janneshb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/janneshb/32/36035_2.png) [@janneshb](https://discourse.julialang.org/u/janneshb)\
**Post date:** [October 23, 2024, 3:07pm UTC](https://discourse.julialang.org/t/how-to-update-non-time-dependent-states-in-mixed-ode-problems/121594/13 "2024-10-23T15:07:58Z")

</div>

If`q_in` should be continuously updated then it is in fact _not_ a discrete variable.  
I feel like the problem here is not an implementation issue, but a lack of understanding of how the underlying physical system works. Maybe it’s just on my side though.

To summarize: only use `du[...] = ` if you are writing out an ODE, i.e.

\frac{du}{dt} = ...

Use a discrete or periodic callback if you have something like

u\_{k+1} = ...

where the you know the time that passes between k and k+1.

If you just have a continuous variable like

```julia
q_in = q_out + inout_itp(t)

```

then you can simply use that in your `discharge!`function, but `q_in` should _not_ be a state of the system.

Hope this helps.