Simplest thing to do would be to just use the interpolation of the derivative integ(t,Val{1}) (or its in-place form integ(du,t,Val{1})) and use a continuous callback to say when it’s 0.
One other thing that’s worth mentioning is that local extrema of the diffeq are roots of the input function. As such, you could first find all the roots, and then just saveat at those times.
ERROR: MethodError: no method matching sign(::Vector{Float64})
Closest candidates are:
sign(::Union{Static.StaticBool{N}, Static.StaticFloat64{N}, Static.StaticInt{N}}) where N at C:\Users\mihim\.julia\packages\Static\Ldb7F\src\Static.jl:404
sign(::Unsigned) at number.jl:163
sign(::Unitful.AbstractQuantity) at C:\Users\mihim\.julia\packages\Unitful\fbiNW\src\quantities.jl:453
...
How can I output the state varible when integ(t,Val{1})== 0?
How about extract the discrete 1d time series u from the ODE solver as a vector, use findmax to find the maximum value, form a local polynomial approximation for u(t) from a data points spanning the maximum (with polynomial degree matching the order of the ODE integrator), and then find the maximum of the polynomial over continuous t? Seems conceptually straightforward and requires only a few lines of code given Polynomials.jl.
I use a similar process to find Poincare intersections of ODEs. I think it took about six lines of code.
Yes, of course we can do that for normal periodic solutions. but for period-2 limit cycles this mothed would not work( period-2 limit cycle would have like two envelopes. That is the reason why I try to find the local extrema during the process.
Yes, this method will work, because findmax will find the global maximum of the discrete time series. For example, y(x) = \exp(-x^2) \cos 16x has many local maxima over [-1,1], but if you sample y discretely over this range and apply findmax, it finds the global maximum at x=0.
My mistake: I see you are looking for local extrema rather than the global. But it’s still simpler to do this with general-purpose Julia code on an extracted discrete time series than figuring out how to write the appropriate callback code for a library. Extract the discrete time series u, approximate du/dt with dudt = u[2:end]-u[1:end-1] (dividing by 1/\Delta t if desired), then look for indices where dudt goes from positive to negative, and do local polynomial approximation around those indices. Should be doable with ten or so lines of code.
using LinearAlgebra, Polynomials
"""
findlocalmaxima(t, u::Vector{T}, dudt::Vector{T}) where T <: Real
Compute and return the local maxima of a time series `u`. Inputs are discrete
samples of continuous variables with `u[n]` = u(`t[n]`) and similarly for dudt.
Outputs are vectors tmax, umax specifying the t and u values of local maxima of u(t).
Examples:
t = range(-5pi, 5pi, 64)
u = cos.(t)
dudt = -sin.(t)
tmax, umax = findlocalmaxima(t,u,dudt)
"""
function findlocalmaxima(t, u::Vector{T}, dudt::Vector{T}) where T <: Real
N = length(t)
@assert length(u) == N
@assert length(dudt) == N
# find indices of du/dt zero crossings downwards, indicating local maxima of u
signchange = [false; [dudt[i] >= 0 && dudt[i+1] < 0 for i in 1:N-1]]
maxlocations = findall(signchange)
Nmaxima = length(maxlocations)
umax = zeros(Nmaxima) # u values for local maxima of u(t)
tmax = zeros(Nmaxima) # t values for local maxima of u(t)
s = [-1; 0; 1] # auxiliary time variable for nbrhd of local max, s = (t-t[n])/Δt
for i in 1:Nmaxima
n = maxlocations[i] # index of local maximum of u data
upoly = fit(s, u[n .+ s]) # fit local polynomial u(s) to s,u data, order set by length(s)
tpoly = fit(s, t[n .+ s]) # fit local polynomail t(s) to s,t data (overkill for uniform time steps)
dudspoly = derivative(upoly) # compute du/ds(s)
smax = roots(dudspoly)[1] # solve du/ds(s) = 0 for s
umax[i] = upoly(smax) # evaluate polynomial u(s)
tmax[i] = tpoly(smax) # evaluate polynomial t(s)
end
return tmax, umax
end
The order of the local polynomial approximation is set by the size of s. I hardcoded s = [-1; 0; 1] to produce local quadratic approximation of u(t) around t=t[n] using data u[n .+ s] = u[n-1; n; n+1]. Production code should take the polynomial order as an input and set the vector ‘s’ accordingly.
Production code should have checks that u[n .+ s] doesn’t go out of bounds.
Maybe should be named interpolatelocalmaxima or something like that.