# How to find the local extrema of the solutions when solving a differential equation？

Hello guys,

I got a problem when I want to calculate the local extrema in an ODE. I did the following steps:

1. define the ODE problem;
2. define a integrator integ;
3. take 2000 successful step ahead to insure the trajectory is on the steady state, for example step!(integ, 2000).
4. start a loop:
maxi = []
minm = []
@time for i in 1:1600000
temp = get_du(integ)
step!(integ)
if temp[1] < 0.0 && temp[1] * get_du(interds)[1] < 0.0
push!(maxi,integ.u[1])
elseif temp[1] < 0.0 && temp[1] * get_du(interds)[1] > 0.0
push!(minm,integ.u[1])
else
continue
end
end


but it seems very slow then running the loop and also I will get wrong results. So I ask for help for this problem.

many thanks!

I think the best way to do this is to add a variable to your equation that tracks the derivative and add a Event Handling and Callback Functions · DifferentialEquations.jl to track when the derivative is zero.

2 Likes

@ChrisRackauckas can I use callback to obtain the du from the last step and output some variables during the step!(integer)?

I also think about the callback, but how can I get du from the last step?

Yes.

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.

2 Likes

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.

1 Like

For the equilibrium points, yes, but I am trying to calculate the local extrema of periodic solusions.

Hi Chris, still have some questions, I apply the callback like this:

maxi = []
minm = []
function condition(u,t,integrator)
integrator(t,Val{1})
end
function affect!(integrator,maxi)
maxi = vcat(maxi, integrator.u')
end
cb = ContinuousCallback(condition,affect!)
diffeq = (alg=Tsit5(), adaptive=false, dt=0.001, callback=cb)


and then use the step! like following:

integ = integrator(ds; diffeq)
step!(interds, 2000)


But error comes like this:

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.

julia> x = range(-1, 1, 201);

julia> y = cos.(16x) .* exp.(-x.^2);

julia> (ymax, i) = findmax(y)
(1.0, 101)

julia> x[i]
0.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.

dudt = u[2:end]-u[1:end-1]

Yes you are right, this can be a solution to my problem, I would try that!

Still, I want to explore the potentials of callbacks and see what It can do during the Iterations.

1 Like

Here’s a working function

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


and usage

julia> t = range(-5pi, 5pi, 256);
julia> u = cos.(t);
julia> dudt = -sin.(t);
julia> tmax, umax = findlocalmaxima(t, u, dudt);

julia> tmax
5-element Vector{Float64}:
-12.566370614359172
-6.283185307179586
-6.938893903907228e-18
6.283185307179586
12.566370614359172

julia> umax
5-element Vector{Float64}:
0.999994607368687
0.999994607368687
0.999994607368687
0.999994607368687
0.9999946073686868


• 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.