# Problem about using DiscreteCallback function in solving ODE problem

Hello everyone,

I am a PostDoc in Department of Radiology, University of Minnesota. I want to use Julia to solve an ODE problem, where the user-provided parameter `p` is changed dynamically. I refered to several topics in this forum(https://discourse.julialang.org/t/diffeqs-hybrid-continuous-discrete-system-periodic-callback/23791/2), and I tried to use the function `DiscreteCallback` to dynamically update `integrator.p`. However, in my code the parameter p is not updated successfully. It is always equal to the initial value. I attached my code below. Is there something wrong in it?
Thanks a lot for your kind help!

Xiaodong

``````using DifferentialEquations

using MATLAB

using Plots

mf = MatFile("C:/Users/maxiao/.julia/dev/BlochSimWithDiffEq/test/testdata_blochsimu.mat")

Brt = get_variable(mf, "Brt")

Δt = get_variable(mf, "dt")

Brt = Brt[:,1000,:]

B = [0.0, 0.0, 0.0]

# Define the Bloch equation

function BlochEq!(du, u, p, t)

du = γ/(2*π) * u * p - γ/(2*π) * u * p

du = γ/(2*π) * u * p - γ/(2*π) * u * p

du = γ/(2*π) * u * p - γ/(2*π) * u * p

end

# Constant parameters

γ = 2.675e8

# Define callback function for discrete control

Nt = 736

tstops = collect(0:Δt:(Nt-1)*Δt)

function condition_control_loop(t,u,integrator)

(t in tstops)

end

function control_loop!(integrator)

nt = round(integrator.t/Δt+1)

integrator.p = Brt[:,nt]

end

cb = DiscreteCallback(condition_control_loop, control_loop!)

# Initial conditions

u0 = [0.0, 0.0, 1.0]

prob = ODEProblem(BlochEq!, u0, (0.0, (Nt-1)*Δt), B)

sol  = solve(prob, Tsit5(), callback = cb, tstops=tstops)
``````

should be `(u,t,integrator)`

BTW, you might want look into using https://diffeq.sciml.ai/latest/features/callback_library/#PeriodicCallback-1

1 Like

Thank you Chris! That is amazingly fast.
Actually I tried to change the control condition function to `(u,t,integrator)`, but it failed. I attached the error message below. It is weired since `(u,t,integrator)` should be the correct order based on the document.
I will look into the link you suggested.

``````julia> sol  = solve(prob, Tsit5(), callback = cb, tstops=tstops)
ERROR: ArgumentError: invalid index: 2.0 of type Float64
Stacktrace:
 to_index(::Float64) at .\indices.jl:297
 to_index(::Array{Float64,2}, ::Float64) at .\indices.jl:274
 to_indices at .\indices.jl:325 [inlined]
 to_indices at .\multidimensional.jl:711 [inlined]
 to_indices at .\indices.jl:321 [inlined]
 getindex(::Array{Float64,2}, ::Function, ::Float64) at .\abstractarray.jl:980
 control_loop!
...
...
D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\REPL\src\REPL.jl:86
 run_backend(::REPL.REPLBackend) at C:\Users\maxiao\.julia\packages\Revise\AMRie\src\Revise.jl:1023
 top-level scope at none:0
``````

that’s because this is a Float64, not an integer. Try `nt = round(Int,integrator.t/Δt+1)`

Yes it works now. Thanks a lot!

1 Like