Differentialequations.jl - modify `u` in `f` possible without callback?

General Usage
, ,
1

Dear all,

I’m wondering if it is possible to set the value of some states to a known value from within the rhs-function f? E.g. by setting u to the desired value and du to 0. Without modfication of u it seems difficult to set a du that achieves the desired effect because it would depend on the length of the current time step.

See below a MWE, where I try to achieve this behavior. But the suggested code doesn’t update u:

using DifferentialEquations
using Plots

# check if it is possible to update state vector in f-function
u0 = [0;0;0;0]
tspan = (0,10)
p = []
function f(du, u, p, t)
    du .= [1; 2; 0.5; -1]
end
function f2(du, u, p, t)
    if t > 2 # update state vector
        u[3] = 12 # this has no effect on u[3]
        u[4] = 14 # this has no effect on u[4]
        du .= [1; 2; 0; 0]
    else
        du .= [1; 2; 0.5; -1]
    end
end

ode = ODEProblem(f,u0,tspan,p)
sol = solve(ode, Tsit5())
ode2 = ODEProblem(f2,u0,tspan,p)
sol2 = solve(ode2, Tsit5())

plot(plot(sol), plot(sol2); layout = (2,1))

Above code does not modify the solution to 12 and 14 for the 3rd and 4th states:
plot_6

Thanks for any suggestions.

2

As far as I know this isn’t possible, but it seems a perfect case for a callback? What’s your reason for not wanting to use one?

3

I don’t think there’s any integrator that supports that, in any language, even at a mathematical level. It doesn’t quite make sense because time doesn’t always go forwards: what if a step is rejected? This is one of the quintessential cases for using a DiscreteCallback: use PresetTimeCallback at 2.0 to do this change.

4

Dear Thomas and Chris,

Thank you for your replies! That’s what I feared/expected. My main reason for not wanting to use a callback is mostly lazyness in refactoring a larger codebase. :sweat_smile: At least good to have certainty that it’s not doable like I intended a first.

Background:
I have re-implemented an existing code solving water mass balance across a forest ecosystem, see figure below. There are intermittent storage compartments such as snow cover or intercepted rain on leaves, that can disappear (represented by 0 mm of water in the state vector u). I am now adding the transport of concentrations with these waters fluxes. As a concentration is undefined when there is no volume I need to come up with ways to handle certain edge-cases.
Updating u in f would have been a quick fix e.g. for the first snow that falls, where the concentration of the storage will just be equal to the incoming flux without any mixing.

I’m afraid I’m not gonna be able to avoid major code refactors to address these concentration related issues. I have ideas how to go about this… it’s just more work. :smiley:

LWFBrook90jl_overview_titleLogo
LWFBrook90jl_overview_titleLogo1154×1300 163 KB

5

The solution suggested by Chris with the PresetTimeCallback would look like this for the MWE:

using DifferentialEquations
using Plots
u0 = [0;0;0;0]
tspan = (0,10)
p = []
function f(du, u, p, t)
    du .= [1; 2; 0.5; -1]
end
affect!(integrator) = integrator.u[[3, 4]] .= [12, 14]
cb = PresetTimeCallback(2.0,affect!)
ode3 = ODEProblem(f, u0, tspan, p;
    callback = cb)
sol3 = solve(ode3, Tsit5())
plot(sol3)

plot_4