# Differential equation with periodic orbit: Can Diffeq.jl determine the period?

**URL:** <https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982>\
**Category:** General Usage\
**Created:** [December 23, 2023, 6:38pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982 "2023-12-23T18:38:43Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![HMegh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hmegh/32/216684_2.png) [@HMegh](https://discourse.julialang.org/u/HMegh)\
**Post date:** [December 23, 2023, 6:38pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/1 "2023-12-23T18:38:43Z")

</div>

Hi, I am working on a problem where the differential equation \mathbf{x}'(t)=F(\mathbf{x}(t)) admits a periodic solution. For example:

\begin{cases} x'(t) = -y(t) \\ y'(t)=x(t)\end{cases},\qquad x(0)=1,\ y(0)=0. 

The solution in this case is (x(t),y(t))= (\cos(t),\sin(t)) which has a period of 2\pi. The following code solves the ODE and more importantly, the solution `sol` provides a good parametrization of the unit circle.

```julia-auto
using DifferentialEquations
f(u,p,t) = [-u[2],u[1]]
tspan= (0.0,2π)
u0=[1.0,0.0]
prob=ODEProblem(f,u0,tspan)
sol=solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

using Plots
plot(sol,linewidth=5,idxs=(1,2),label="")

```

My question: Say we did not know that the period is 2\pi. Is there a way to determine it?

I tried a naive approach (see below) and I was wondering if there is a better way to do it, maybe using DifferentialEquations.jl directly

* * *

## Naive approach

Solve the ODE on a _large_ interval, then find a solution `guess` to \mathbf{x}(t)=\mathbf{x}\_0. This solution is a multiple of the period, then we find the largest integer k such that `guess/k` is a solution too, which gives us the period

```julia-auto

using DifferentialEquations
f(u,p,t) = [-u[2],u[1]]
tspan= (0.0,2π)
u0=[1.0,0.0]
prob=ODEProblem(f,u0,tspan)
sol=solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

#naive approach
big_tspan=(0.0,1000.0)
big_prob=ODEProblem(f,u0,big_tspan)
big_sol=solve(big_prob,Tsit5(),reltol=1e-8,abstol=1e-8)

using LinearAlgebra
#g(t)=0 when 𝐱(t)=𝐱(0)
g(t)=norm(big_sol(t)-u0)^2
#g'(t)= 2*sol'(t)*(sol(t)-u0)=2*f(sol(t))* (sol(t)-u0)
#The ODE is autonomous so f(x,p,t)=f(x,0,0)
gp(t)= 2*dot(f(big_sol(t),0.0,0.0), (big_sol(t)-u0))

#Newton's method
guess=big_tspan[2]/2 #initial guess
while norm(g(guess))>1e-8
    guess=guess .- g(guess)/gp(guess)
end
#guess is a multiple of the period: guess= period*k (assume k<1000)
k_list=range(1,1000)
possible_periods= guess ./ k_list
k=findlast(g.(possible_periods) .< 1e-8)

period=guess/k

```

This returns a good approximation of the period: `period-2\pi =8e-7`.

---

<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:** [December 23, 2023, 7:49pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/2 "2023-12-23T19:49:28Z")

</div>

A good way to find periodic orbits is via a bifurcation tool like BifurcationKit.jl. Have you tried that? Pure simulation is one way but that might be a bit easier.

---

<div class="post-metadata">

**Author:** ![HMegh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hmegh/32/216684_2.png) [@HMegh](https://discourse.julialang.org/u/HMegh)\
**Post date:** [December 23, 2023, 10:09pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/3 "2023-12-23T22:09:21Z")

</div>

Thanks for the answer, I will take a look at it.

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [December 23, 2023, 10:50pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/4 "2023-12-23T22:50:47Z")

</div>

I you already know that you look for a periodic solution, why don’t you linearize your system and compute the eigensolutions ?

---

<div class="post-metadata">

**Author:** ![HMegh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hmegh/32/216684_2.png) [@HMegh](https://discourse.julialang.org/u/HMegh)\
**Post date:** [December 24, 2023, 8:33pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/5 "2023-12-24T20:33:02Z")

</div>

Can you expand more? would a ‘local’ linearization be enough to find the period?

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [December 26, 2023, 7:27pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/7 "2023-12-26T19:27:16Z")

</div>

Looking at your equations, you can see that your system is linear, conservative, and leads to periodic orbits.

One way to find the natural frequencies, and therefore periods of periodic solutions is to linearize the underlying system of equations (which is straightforward here, but I’ll show how to do it with MTK), and compute the eigenvalues of the state matrix of the state-space model of the system.

```julia
using ModelingToolkit
using DifferentialEquations, ControlSystemsBase

@variables begin
    t
    x(t)
    y(t)
end

@parameters begin
    p
end
D = Differential(t)

eqs = [
    D(x) ~ -y
    D(y) ~ x
]

@named model = ODESystem(eqs, t, [x, y], [p])
sys = structural_simplify(model)

tspan = (0.0, 2π)
@nonamespace initial_conditions = Dict([x => 1.0, y => 0.0])

prob = ODEProblem(sys, initial_conditions, tspan,[p=>0])
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

(; A, B, C, D), simplified_sys = linearize(sys, [],[x,y], op = merge(initial_conditions,Dict([p=>0])))

state_sys = ss(A,B,C,D)

dampreport(state_sys)

Wn, zeta, ps = damp(state_sys)

```

This gives 2 times 1 rad/s periods, which correspond to `2pi` periods.

@ChrisRackauckas I do not know why, maybe it is a bug, but the `linearization` does not works without a dummy parameter (with some `NullParameter` instead).

---

<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:** [December 26, 2023, 8:58pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/8 "2023-12-26T20:58:10Z")

</div>

> [@BambOoxX](#):
>
> @ChrisRackauckas I do not know why, maybe it is a bug, but the `linearization` does not works without a dummy parameter (with some `NullParameter` instead).

That’s worth an issue. I presume it just was tested with parameters each time.

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [December 28, 2023, 4:16pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/9 "2023-12-28T16:16:59Z")

</div>

For the record, the problem with a parameter-less model is tracked in [this issue](https://github.com/SciML/ModelingToolkit.jl/issues/2400)

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [December 28, 2023, 5:40pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/10 "2023-12-28T17:40:17Z")

</div>

Note that if the system is nonlinear, linearization is not very useful for finding periodic orbits. If you’re sitting right at a Hopf bifurcation, the linearized system will have periodic orbits, but the actual system may not. If the system has a periodic orbit that’s small and close to being annihilated in a Hopf bifurcation, you can expand to cubic order and extract the period from the normal form coefficients, see [Hopf bifurcation - Wikipedia](https://en.wikipedia.org/wiki/Hopf_bifurcation#Supercritical_and_subcritical_Hopf_bifurcations). Otherwise, you need more sophisticated methods, and BifurcationKit.jl is probably a good place to start.

---

<div class="post-metadata">

**Author:** ![BambOoxX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bambooxx/32/22179_2.png) [@BambOoxX](https://discourse.julialang.org/u/BambOoxX)\
**Post date:** [December 28, 2023, 10:09pm UTC](https://discourse.julialang.org/t/differential-equation-with-periodic-orbit-can-diffeq-jl-determine-the-period/107982/11 "2023-12-28T22:09:37Z")

</div>

Yes, linearization is the very first step here. Bifurcation analysis is surely needed for nonlinear systems.