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

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