CODEWITHSUNDEEP
juliadifferential-equationsdifferentialequations.jl
winkmal
winkmal

Reputation: 632

How to read corresponding value at specified time in DiscreteCallback?

Similar to this question, I am trying to solve this ODE with a time-dependent input parameter. It consists of a series of discrete callbacks. At certain times, a parameter is changed (not a state!). Times and values are stored in a nx2 Array. But I can't get the affect function to find the corresponding parameter value at the specified time. In the given examples, the value assigned to u[1] is usually constant. Consider this MWE (with a very Matlab-like approach), which works correctly without the callback:

using DifferentialEquations
using Plots

function odm2prod(dx, x, params, t)
    k_1, f_1, V_liq, X_in, Y_in, q_in = params

    rho_1 = k_1*x[1]
    q_prod = 0.52*f_1*x[1]
    # Differential Equations
    dx[1] = q_in/V_liq*(X_in - x[1]) - rho_1
    dx[2] = q_in/V_liq*(Y_in - x[2])
end

x0      = [3.15, 1.5]
tspan   = (0.0, 7.0)
params  = [0.22, 43, 155, 249, 58, 0]
prob    = ODEProblem(odm2prod, x0, tspan, params)

input   = [1.0 60; 1.1 0; 2.0 60; 2.3 0; 4.0 430; 4.05 0]
dosetimes = input[:,1]
function affect!(integrator)
    ind_t = findall(integrator.t == dosetimes)
    integrator.p[6] = input[ind_t, 2]
end
cb = PresetTimeCallback(dosetimes, affect!)
sol = solve(prob, Tsit5(), callback=cb, saveat=1/12)

plot(sol, vars=[1, 2])

It does not work. The error originates at line 22, since comparing a vector to a scalar seems not to be defined in Julia, or there is a special syntax I am unaware of.

I know that it is possible to use time-dependent parameters in Julia, but I suppose that would only work for continuous functions, not discrete changes!? I haven taken a look at the help for interpolate, but I am not sure how to use it for my specific case.

Could someone tell me how to get this to work, please? Should probably need just a few lines of code. Also, I do not necessarily want dosetimes as part of sol.t, unless they coincide.

Upvotes: 3

Views: 536

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

You are using findall wrong, the documentation says

findall(f::Function, A)

Return a vector I of the indices or keys of A where f(A[I]) returns true.

Then you have to take into account that the result of a search for "all" is a list. As you expect it to only have one element, use the first one only

function affect!(integrator)
    ind_t = findall(t -> t==integrator.t, dosetimes)
    integrator.p[6] = input[ind_t[1], 2]
end

and you get the plot

plot of the ode solution

Upvotes: 3

Related Questions