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

General Usage
1

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.

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


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.

1 Like
2

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.

1 Like
3

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

4

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

5

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

7

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.

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

1 Like
8

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

9

For the record, the problem with a parameter-less model is tracked in this issue

1 Like
10

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. Otherwise, you need more sophisticated methods, and BifurcationKit.jl is probably a good place to start.

11

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