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

**URL:** <https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628>\
**Category:** General Usage\
**Tags:** question, differentialequation\
**Created:** [February 14, 2023, 3:49pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628 "2023-02-14T15:49:30Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![YiboXia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiboxia/32/50734_2.png) [@YiboXia](https://discourse.julialang.org/u/YiboXia)\
**Post date:** [February 14, 2023, 3:49pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/1 "2023-02-14T15:49:30Z")

</div>

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:

```julia
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!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 14, 2023, 3:56pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/2 "2023-02-14T15:56:02Z")

</div>

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](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#SciMLBase.ContinuousCallback) to track when the derivative is zero.

---

<div class="post-metadata">

**Author:** ![YiboXia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiboxia/32/50734_2.png) [@YiboXia](https://discourse.julialang.org/u/YiboXia)\
**Post date:** [February 14, 2023, 3:58pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/3 "2023-02-14T15:58:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![YiboXia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiboxia/32/50734_2.png) [@YiboXia](https://discourse.julialang.org/u/YiboXia)\
**Post date:** [February 14, 2023, 3:59pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/4 "2023-02-14T15:59:52Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 14, 2023, 4:14pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/5 "2023-02-14T16:14:28Z")

</div>

> [@YiboXia](#):
>
> can I use `callback` to obtain the `du` from the last step and output some variables during the `step!(integer)`?

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 14, 2023, 6:07pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/6 "2023-02-14T18:07:42Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![YiboXia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiboxia/32/50734_2.png) [@YiboXia](https://discourse.julialang.org/u/YiboXia)\
**Post date:** [February 15, 2023, 10:13am UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/7 "2023-02-15T10:13:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![YiboXia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiboxia/32/50734_2.png) [@YiboXia](https://discourse.julialang.org/u/YiboXia)\
**Post date:** [February 15, 2023, 1:24pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/8 "2023-02-15T13:24:59Z")

</div>

> [@ChrisRackauckas](#):
>
> integ(t,Val{1})

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

```julia
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:

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

```

But error comes like this:

```julia
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`?

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [February 15, 2023, 1:57pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/9 "2023-02-15T13:57:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![YiboXia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiboxia/32/50734_2.png) [@YiboXia](https://discourse.julialang.org/u/YiboXia)\
**Post date:** [February 15, 2023, 2:52pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/10 "2023-02-15T14:52:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [February 15, 2023, 5:59pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/11 "2023-02-15T17:59:18Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [February 15, 2023, 8:47pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/12 "2023-02-15T20:47:24Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![YiboXia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiboxia/32/50734_2.png) [@YiboXia](https://discourse.julialang.org/u/YiboXia)\
**Post date:** [February 17, 2023, 11:36am UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/13 "2023-02-17T11:36:39Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [February 17, 2023, 2:15pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/14 "2023-02-17T14:15:20Z")

</div>

Here’s a working function

```julia
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
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

```

EDIT Some comments:

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

---

<div class="post-metadata">

**Author:** ![YiboXia](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yiboxia/32/50734_2.png) [@YiboXia](https://discourse.julialang.org/u/YiboXia)\
**Post date:** [May 8, 2023, 1:19pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/15 "2023-05-08T13:19:32Z")

</div>

oops…I forgot to press the reply button until now. Sorry and really thanks for your code!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 8, 2023, 3:48pm UTC](https://discourse.julialang.org/t/how-to-find-the-local-extrema-of-the-solutions-when-solving-a-differential-equation/94628/16 "2023-05-08T15:48:40Z")

</div>

This function also exists in EasyModelAnalysis.jl, a helper library for SciML tools.

[https://docs.sciml.ai/EasyModelAnalysis/stable/api/basic\_queries/#Extrema](https://docs.sciml.ai/EasyModelAnalysis/stable/api/basic_queries/#Extrema)

See the tutorial:

[https://docs.sciml.ai/EasyModelAnalysis/stable/getting\_started/](https://docs.sciml.ai/EasyModelAnalysis/stable/getting_started/)