Problem about using DiscreteCallback function in solving ODE problem

General Usage
,
1

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

# load variables

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[1] = γ/(2*π) * u[2] * p[3] - γ/(2*π) * u[3] * p[2]

  du[2] = γ/(2*π) * u[3] * p[1] - γ/(2*π) * u[1] * p[3]

  du[3] = γ/(2*π) * u[1] * p[2] - γ/(2*π) * u[2] * p[1]

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)
2

should be (u,t,integrator)

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

1 Like
3

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:
 [1] to_index(::Float64) at .\indices.jl:297
 [2] to_index(::Array{Float64,2}, ::Float64) at .\indices.jl:274
 [3] to_indices at .\indices.jl:325 [inlined]
 [4] to_indices at .\multidimensional.jl:711 [inlined]
 [5] to_indices at .\indices.jl:321 [inlined]
 [6] getindex(::Array{Float64,2}, ::Function, ::Float64) at .\abstractarray.jl:980
 [7] control_loop!
...
...
D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.4\REPL\src\REPL.jl:86
 [20] run_backend(::REPL.REPLBackend) at C:\Users\maxiao\.julia\packages\Revise\AMRie\src\Revise.jl:1023
 [21] top-level scope at none:0
4

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

5

Yes it works now. Thanks a lot!

1 Like