Terminate condition in ContinuousCallback in DifferentialEquations.jl

General Usage
, , ,
1 

 using DifferentialEquations
 u0 = [1.,0.]
harmonic! = @ode_def HarmonicOscillator begin
          dv = -x
          dx = v
       end
tspan = (0.0,10.0)
 prob = ODEProblem(harmonic!,u0,tspan)
function terminate_affect!(integrator)
           terminate!(integrator)
end
function terminate_condition(u,t,integrator)
           u[2] < 0.5
end

terminate_cb = ContinuousCallback(terminate_condition,terminate_affect!)
julia> sol = solve(prob,callback=terminate_cb)
retcode: Terminated
Interpolation: automatic order switching interpolation
t: 8-element Vector{Float64}:
 0.0
 0.0009990009990009992
 0.010989010989010992
 0.07985922249873038
 0.2403882280626971
 0.48125583199780264
 0.7872724750204353
 0.7872724750204353
u: 8-element Vector{Vector{Float64}}:
 [1.0, 0.0]
 [0.9999995009985435, 0.000999000832833342]
 [0.9999396214263468, 0.010988789821184269]
 [0.9968129466114029, 0.07977436592568168]
 [0.9712456180254766, 0.23807970964822978]
 [0.8864142948207253, 0.4628927349710369]
 [0.7057801319307302, 0.7084309070362929]
 [0.7057801319307302, 0.7084309070362929]

why is u[2] < 0.5 condition not satisfied here? Am I missing something?

2

Certainly definition of terminate_affect! function is missing.

3

@zdenek_hurak updated the code. Still i get same output .

4

I guess my way of writing condition was wrong

function terminate_condition(u,t,integrator)
                 u[2] - 0.5
       end
terminate_condition (generic function with 1 method)

this seems to give expected result

5

A minor issue first: you can just easily write

terminate_affect!(integrator) = terminate!(integrator)

Now for the major issue: why aren’t you happy with the result you get? If I run the code with the terminate_condition relaxed, that is, I set

function terminate_condition(u,t,integrator)
           u[2] < -10e6
end

and I plot the results (but just the points that were obtained by the solver, nothing interpolated afterwards), I get this:
discourse

Note that the ninth component of the simulated sequence of u[2]s (in the figure labelled as x) is already below 0.5. Hence the eighth u[2] was the last one that satisfied the condition.

For convenience the whole code here:

using DifferentialEquations
u0 = [1.,0.]
harmonic! = @ode_def HarmonicOscillator begin
          dv = -x
          dx = v
      end
tspan = (0.0,10.0)
prob = ODEProblem(harmonic!,u0,tspan)
function terminate_condition(u,t,integrator)
           u[2] < -10e6
end
terminate_affect!(integrator) = terminate!(integrator)
terminate_cb = ContinuousCallback(terminate_condition,terminate_affect!)
sol = solve(prob,callback=terminate_cb,dense=true)

using Plots
scatter(sol,denseplot=false)

savefig("discourse.png")
6

I was expecting the result to be terminated at the 4th redpoint, as that would be first time x is crossing 0.5. Am i right ?

7

I see. Well, the documentation for callbacks says that

ContinuousCallback is applied when a continuous condition function hits zero.

I interpret it that it detects change. And although u[2] actually started while not satisfying the condition, nothing changed because “zero was not hit” or even better “zero was not crossed”.

1 Like
8

I guess thats why u[2] - 0.5
condition worked. Thanks for pointong out the definition

9

Yes, ContinuousCallback has a continuous definition to its conditions, and it will apply the callback at the point at which the condition hits zero. DiscreteCallback is a true/false callback and applies at the end of each step based on true/false. Choose the one that is appropriate for your situation.

1 Like
10

@ChrisRackauckas , may I use this as an opportunity to give a comment on the terminology, in particular the meaning of the continuity and discreteness of the callbacks? Essentially I am asking for confirmation that I understand it correctly. The manual gives

The ContinuousCallback is applied when a continuous condition function hits zero. This type of callback implements what is known in other problem solving environments as an Event. A DiscreteCallback is applied when its condition function is true .

I think I have more or less understood the continuous case – a solver is checking if a function that is expected to evolve continuously hits (or perhaps crosses?) zero. This hitting can happen even within an integration step and the algorithm must be able to detect it and adjust the step accordingly.

But the discrete case is only starting to be clear to me after the definition that you give here – that the condition is only evaluated at the end of the step. Obviously this must be computationally cheaper than the former.

I am now thinking that perhaps the lucid explanation you give here should be incorporated to the part of the manual that I quoted above. How about extending the last sentence as

A DiscreteCallback is applied when its condition function evaluated at the end of the integration step is true .

?

11

That would be a good change.

12

PR submitted.

13

@zdenek_hurak , @ChrisRackauckas
Any directions on how to terminate integration when the derivative of a particular state is nearly zero … say du[1] ~ 0 in the above case?
I am a little confused because there could be multiple time instances when the derivative becomes zero, but I want the time when steady-state for u[1] is reached(which could be different from when u[2] reaches the steady-state in multiscale systems)

14

Look at terminate and the SteadyStateCallback.

1 Like