# Numerical continuation and periodic orbits

Apologies if this question is poorly worded, I’m very new to numerical continuation and nonlinear frequency analysis. I have a nonlinear system of ODEs, and I’d like to numerically find periodic solutions for a changing parameter. Ideally, I’d like to recreate Figure 1 in \verb|[1]| for the Duffing equation. I haven’t included that figure here, I don’t think it’s necessary. Basically, given a nonlinear system, I want to find a periodic (steady state) response for a range of (single) parameter inputs. Looking at the Duffing Equation, which is basically a forced spring-mass damper with a nonlinearity thrown in, I’d like to keep all parameters constant except for \omega, and then numerically find a one-period-long time-domain solution to the nonlinear dynamics for each value of \omega.

#### Duffing Equation

\ddot{z} + c \dot{z} + k z + \alpha z^3 = A \cos \left(\omega t \right)

#### Duffing Equation (redefined)

• We can get rid of t by rewriting the dynamics as a fourth-order autonomous system

\begin{align} \dot{z}_1 &= z_2 \\ \dot{z}_2 &= -k z_1 - c z_2 -\alpha z_1^3 + A z_4 \\ \dot{z}_3 &= z_3 + \omega z_4 - z_3 \left(z_3^2 + z_4^2 \right) \\ \dot{z}_4 &= -\omega z_3 + z_4 - z_4 \left( z_3^2 + z_4^2 \right) \end{align}

Note that the solution of z_4 is equivalent to \cos (\omega t), as shown in the \dot{z_2} expression.

For many values of \omega (and all other parameters held constant), I would like to find a periodic solution for the steady state response of the system. I’m trying to accomplish this with the code below.

using Setfield
using LinearAlgebra
using BifurcationKit
using ModelingToolkit
using DifferentialEquations

@variables t z[1:4](t)
@parameters c k α A ω
@derivatives D'~t

duffing, f, J = let

eqs = [
D(z[1]) ~ z[2],
D(z[2]) ~ -k*z[1] - c*z[2] - α*z[1]^3 + A*z[4],
D(z[3]) ~ z[3] + ω*z[4] - z[3]*(z[3]^2 + z[4]^2),
D(z[4]) ~ -ω*z[3] + z[4] - z[4] * (z[3]^2 + z[4]^2)
]

# Build ModelingToolkit system
duffing = ODESystem(eqs, t, z, [ω, A, c, k, α]; name=:Duffing)

# Build function for dynamics
func = eval(first(generate_function(duffing)))
f = (x,u,t=0)->func(x, (u), t)

# Build function for Jacobian
jac = eval(first(generate_jacobian(duffing)))
J = (x,u,t=0)->jac(x,(u), t)

duffing, f, J

end

# Time slices
t = collect(0.0:0.01:10.0)

# Guess for orbit period
T = 1.0

# Guess for periodic orbit
orbitguess = vcat(sin.(t), cos.(t), ones(length(t) * 2), T)

# Construct Periodic Orbit Trapezoidal Problem
pb = PeriodicOrbitTrapProblem(f, J, 0.0, 0.0, 10)

newton(pb, orbitguess, (2.5, 2.5, 0.2, 1.0, 0.05), NewtonPar())



This code errors on the newton call, because pb seems to think the number of states, N, is zero. Am I setting up PeriodicOrbitTrapProblem incorrectly?

Thanks in advance for the help.

And thanks for the package @rveltz!

### Reference

\verb|[1]| Nguyen, Duc H., Mark H. Lowenberg, and Simon Neild. “Nonlinear Frequency Response Analysis to Inform Aircraft Control Law Design.” In AIAA Scitech 2020 Forum, p. 0605. 2020.

1 Like

Hi,

Thank you for trying out BifurcationKit. I have not coded non-autonomous systems, it would be an easy task because it amounts to removing the phase conditions.

pb = PeriodicOrbitTrapProblem(f, J, 0.0, 0.0, 10)

Note that you can use keyword arguments to call PeriodicOrbitTrapProblem if it helps. According to the docs, you have to pass vectors ϕ, xπ of size M * 4 where M is the number of time sections. These vectors specify a phase condition as explained here

For example, you can pass xπ = zeros(2N) and for ϕ, you pass the time derivative of your orbitguess.

The “funny” thing is that you need this phase condition because you removed time from the equations (but you had no other choice for now)

2 Likes

Thank you for this quick reply! I’ve reviewed the docs again, but I’m not sure I fully understand.

Should ϕ be of size M * N, and xπ be of size M * N?

Or is xπ of size 2*N?

I’m trying the following code for a working example, but I’m getting the error DimensionMismatch("dot product arguments have lengths 4 and 40")

using Setfield
using LinearAlgebra
using BifurcationKit
using ModelingToolkit
using DifferentialEquations

@variables t z[1:4](t)
@parameters c k α A ω
@derivatives D'~t

duffing, f, J = let

eqs = [
D(z[1]) ~ z[2],
D(z[2]) ~ -k*z[1] - c*z[2] - α*z[1]^3 + A*z[4],
D(z[3]) ~ z[3] + ω*z[4] - z[3]*(z[3]^2 + z[4]^2),
D(z[4]) ~ -ω*z[3] + z[4] - z[4] * (z[3]^2 + z[4]^2)
]

# Build ModelingToolkit system
sys = ODESystem(eqs, t, z, [ω, A, c, k, α])

# Build function for dynamics
func = eval(first(generate_function(sys)))
f = (x,u,t=0)->func(Array(x), (u), t)

# Build function for Jacobian
jac = eval(first(generate_jacobian(sys)))
J = (x,u,t=0)->jac(Array(x),(u), t)

duffing, f, J

end

# The orbit guess zᵧ is of size M*N + 1, where zᵧ[end] = 2π (guess at period)
γ  = collect(range(0.0, stop=2π, length=10))                        # "time" vector
zᵧ = vcat(sin.(γ)...,  cos.(γ)...,  cos.(γ)...,  cos.(γ)..., 2π) 	# orbit guess
żᵧ = vcat(cos.(γ)..., -sin.(γ)..., -sin.(γ)..., -sin.(γ)...)        # time derivative of orbit guess

M = length(γ)
N = 4 # number of states in model
pb = PeriodicOrbitTrapProblem(
F=f, J=J,
ϕ=żᵧ, xπ=zeros(N * M), M=M, N=4
)

# Here, [2.5, 2.5, 0.2, 1.0, 0.05] are the desired parameter values
newton(pb, zᵧ, [2.5, 2.5, 0.2, 1.0, 0.05], NewtonPar())



Thanks for the help, I solved the issue by checking out the dev branch!

julia> ]dev BifurcationKit

γ  = collect(range(0.0, stop=2π, length=100))
zᵧ = vcat(sin.(γ)...,  cos.(γ)...,  cos.(γ)...,  cos.(γ)..., 2π)
żᵧ = vcat(cos.(γ)..., -sin.(γ)..., -sin.(γ)..., -sin.(γ)...)

M = length(γ)
N = 4 # number of states in model
pb = PeriodicOrbitTrapProblem(
F=f, J=J,
ϕ=żᵧ, xπ=zeros(N * M), M=M, N=4
)

newton(pb, zᵧ, [2.5, 2.5, 0.2, 1.0, 0.05], NewtonPar())



There seems to be a bug.

You can try

pb = PeriodicOrbitTrapProblem(
F=f, J=J,
ϕ=żᵧ, xπ=zeros(N * M), M=M, N = N
)


while I fix it.

1 Like

Actuellement, it was another constructor

pb = PeriodicOrbitTrapProblem(f, J, żᵧ, zeros(N * M), M)


See also here for how to extract the vector field and the jacobian from your pb: https://github.com/rveltz/BifurcationKit.jl/issues/27

2 Likes

Did you get it to work?

1 Like

Thanks for checking – I did get the code to run with your suggestions (the kwarg N=4 was the key part which made the previous code I posted run, and I’m using the latest version of master).

I’m trying to replicate the following plot from [1] in my previous post.

In this plot, x is the solution to x1.

My code produces the following plot, which is periodic, but isn’t quite as smooth. I’m sure I’m implementing something incorrectly, but I haven’t found the issue so far. I’ve tried a variety of time slices between M=31 and M=351. The code below should be completely self-contained. Can newton plot a periodic, steady state solution for nonlinear dynamics as in Fig 1?


import Pkg
Pkg.add("Plots")
Pkg.develop("BifurcationKit")
Pkg.add("ModelingToolkit")

using Plots
using BifurcationKit
using ModelingToolkit

# Construct system, system function, system Jacobian
duffing, f, J = let

t, z, D, ω, c, k, α, A, ω = let

@variables t z[1:4](t)
@parameters c k α A ω
@derivatives D'~t

t, z, D, ω, c, k, α, A, ω

end;

eqs = [
D(z[1]) ~ z[2],
D(z[2]) ~ -k*z[1] - c*z[2] - α*z[1]^3 + A*z[4],
D(z[3]) ~ z[3] + ω*z[4] - z[3]*(z[3]^2 + z[4]^2),
D(z[4]) ~ -ω*z[3] + z[4] - z[4] * (z[3]^2 + z[4]^2)
]

# Build ModelingToolkit system
sys = ODESystem(eqs, t, z, [ω, A, c, k, α])

# Build function for dynamics
func = eval(first(generate_function(sys)))
f = (x,u,t=0.0)->func(Array(x), (u), t)

# Build function for Jacobian
jac = eval(first(generate_jacobian(sys)))
J = (x,u,t=0.0)->jac(Array(x),(u), t)

sys, f, J

end

γ, sol, histories, converged = let

γ  = collect(range(0.1, stop=2π, length=31))
zᵧ = vcat(cos.(γ)...,  -sin.(γ)...,  cos.(γ)...,  cos.(γ)..., 2π)
żᵧ = vcat(-sin.(γ)...,  cos.(γ)..., -sin.(γ)..., -sin.(γ)...)

N = 4
M = length(γ)

pb = PeriodicOrbitTrapProblem(
F=f, J=J,
ϕ=żᵧ, xπ=zeros(N * M), M=M, N=4
)

output, history, converged, numIter = newton(
pb, zᵧ, [0.24, 2.5, 0.2, 1.0, 0.005],
NewtonPar(tol=1e-12, maxIter=100)
)

z₁ = output[1:M]
z₂ = output[M+1:2M]
z₃ = output[2M+1:3M]
z₄ = output[3M+1:4M]

sol = hcat(z₁, z₂, z₃, z₄)

γ, sol, history, converged

end

# Plot the solution (and the forcing cosine)
plot(γ, hcat(sol[:,1], 2.5*sol[:,4]);
title="Periodic Orbit within Autonomous Duffing Dynamics",
xlabel="Time (seconds)", ylabel="Response",
labels=["z(t)" "Forcing Cosine"])



why do you repeat omega? It seems there are mistakes in the parameters. I would suggest using named tuples.

Also, for the period to be 2pi, you need omega = 1.

As a good test for what you are doing, you should plot the residual:

pb(zᵧ, [1.00, 2.5, 0.2, 1.0, 1.00])[1:end-1] |> plot


It shows that the last two equations, for which you have analytical expression, are not behaving as intended. The guess is wrong (4th component) and it should be ordered in the other way (this is not well documented, so … thank you for spotting!). Basically:

zᵧ = vcat(vec(hcat(cos.(γ),  -sin.(γ),  cos.(γ),  sin.(γ))'), 2pi)


To get the solution, it is easier to call:

sol = BifurcationKit.extractTimeSlices(pb, output)
plot(sol[3,:])


You could also use a shooting method, maybe easier to set up from an ODEProblem.

2 Likes

Ah excellent! Thank you so much for the detailed replies, this has been really helpful.

So should the initial guess order be: x4_guess, x3_guess, x2_guess, x1_guess, period_guess? Or the reverse order, x1_guess, x2_guess x3_guess, x4_guess, period_guess?

As for the extra omega – that was a typo! Thank you for catching it.

Thank you again!

The latter. I updated the docs accordingly. Please keep us posted.

Given what you wrote, it would be easier to just do

function F(x, par)
@unpack c, k, α, A, ω = par
[z[2],
-k*z[1] - c*z[2] - α*z[1]^3 + A*z[4],
z[3] + ω*z[4] - z[3]*(z[3]^2 + z[4]^2),
-ω*z[3] + z[4] - z[4] * (z[3]^2 + z[4]^2)]
end

using Parameters, ForwardDiff, LinearAlgebra

function Fd(z,par)
@unpack ω, A, c, k, α = par
[z[2]*1,
(-k*z[1] - c*z[2] - α*z[1]^3 + A*z[4])*1,
z[3] + ω*z[4] - z[3] * (z[3]^2 + z[4]^2),
-ω*z[3] +   z[4] - z[4] * (z[3]^2 + z[4]^2)]
end

Jd(x,p) = ForwardDiff.jacobian(z->Fd(z,p),x)

γ  = LinRange(0,2pi,101)[1:end-1]
zᵧ = reduce(hcat, [[0cos(t),  -0sin(t),  cos(t),  sin(t)] for t in γ]) |> vec; zᵧ = vcat(zᵧ, 2pi)
żᵧ = reduce(hcat, [[-sin(t),  -cos(t), -sin(t), cos(t)] for t in γ]) |> vec
N = 4
M = length(γ)

par = (ω=1.00, A=0.1, c=0.2, k=1.0, α=1.00)

pb = PeriodicOrbitTrapProblem(
F=Fd, J=Jd,
ϕ=żᵧ, xπ=zeros(N * M), M=M, N=4
)

res = pb(zᵧ, par)[1:end-1]
res = reshape(res, N, M)
plot(res[3:4,:]')

output, history, converged, numIter = newton(
pb, zᵧ, par,
NewtonPar(tol=1e-12, maxIter=100, verbose=true),
)

sol = BifurcationKit.extractTimeSlices(pb, output)
plot(sol[3,:])

2 Likes