Accessing solutions from earlier times in Callback functions

Specific Domains Numerics
,
1

Hi everyone, I’ve only recently started using Julia and I was wondering if it’s possible to access the values obtained at earlier times by a solver in a callback function.

For example, I’m trying to identify maxima (as a function of time) in the solution of a simple oscillator and store them in a vector, but I don’t know how/if it’s possible to access the values of u at previous times to compare them to the value at the current time. I’ve written out the logic I’m trying to achieve in the code below, but I don’t know how to define u_prev or u_prev2 in the condition function. I want for the affect! function to store the current time when u_prev is a local maximum.

using DifferentialEquations, Plots

function condition(u, t, integrator)
	# u_prev = ?    # Value of u[1] at the previous time
	# u_prev2 = ?   # Value of u[1] at the time before the previous time
	u_curr = u[1]
	if (u_curr < u_prev) && (u_prev2 < u_prev)
		return 0
	else
		return 1
	end
end
function affect!(integrator)
	push!(peaks, integrator.t)   # Store current time in peaks vector
end
cb = ContinuousCallback(condition, affect!)

# Initialise problem
u0 = ones(2)
p = (2.0, 1.0)
peaks = zeros(0)
function f!(du,u,p,t)
	a, b = p
	du[1] = a * u[2]
	du[2] = -b * u[1]
end
prob = ODEProblem(f!, u0, (0.0, 50.0), p)
sol = solve(prob, Tsit5(), maxiters = 1e3, callback = cb, abstol=1e-6, reltol=1e-5)

I was thinking that I might be able to access previous values of u through the integrator, but I haven’t been able to figure this out. Similarly, when I try to display the value of integrator.du, only the word “nothing” is outputted.

Any help would be greatly appreciated, and thank you in advance for your help!

1 Like
2

Since u'(t) = 0 and u''(t) < 0 at the maxima, you can write a continuous callback which evaluates those terms, and saves the current state at time t, I think? Your condition instead seems to be trying to find maxima discretely from u_{n-1}, u_n, and u_{n+1}, which is (afaik!) not something the callback interface can do. If you really want to do that, then I would not run it online (with the ODE solve) but analyze the solution array after it is complete.

3

I think that a continuous callback function that calculates when u_1'(t)=0 and u_1''(t)<0 would be ideal, but I haven’t been able to find out how to evaluate the derivatives of u(t) within a callback function. I don’t know how/if it’s possible to explicitly pass the derivatives to the callback function, and I don’t think they’re accessible through the integrator.

Do you know how I can access u_1'(t) and u_1''(t) within the callback function?

4

integrator(t,Val{i}) for the ith derivative.

5

Thanks for the help!

After testing, I’ve found that I can store the times of maxima in u[1] using the functions:

condition(u, t, integrator) = integrator(t, Val{1})[1]
function affect!(integrator)
	if (integrator(integrator.t, Val{2})[1] < 0.0)
		push!(peaks, integrator.t)   # Store current time in peaks vector
	end
end
cb = ContinuousCallback(condition, affect!)

I’ve used the time stored in the integrator (integrator.t) in the affect! function since t is not passed to it explicitly, but do you know exactly what time this will correspond to? Ideally I would want the time at which the condition is satisfied, but I wanted to check whether integrator.t would reference the time at the end of the current interval in this case (as I would expect it to if I displayed it from the condition function).

Also, when using these functions in my original code two values of time are outputted to the peaks vector for each maxima in u[1]. These times are very similar (for example 0.6755106670700617 and 0.6755108454503309 for the first maxima). Is this due to floating point errors in the observed derivatives over short time intervals?