# Numerical continuation and periodic orbits

**URL:** <https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580>\
**Category:** Modelling & Simulations\
**Tags:** question, package\
**Created:** [December 10, 2020, 3:08am UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580 "2020-12-10T03:08:56Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [December 10, 2020, 3:08am UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/1 "2020-12-10T03:08:56Z")

</div>

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.

```julia
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](https://discourse.julialang.org/u/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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 10, 2020, 7:41am UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/2 "2020-12-10T07:41:29Z")

</div>

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](https://rveltz.github.io/BifurcationKit.jl/dev/library/#BifurcationKit.PeriodicOrbitTrapProblem), 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](https://rveltz.github.io/BifurcationKit.jl/dev/periodicOrbitFD/)

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)

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [December 10, 2020, 1:56pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/3 "2020-12-10T13:56:25Z")

</div>

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

> [@rveltz](#):
>
> pass vectors ϕ, xπ of size M \* 4

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

> [@rveltz](#):
>
> you can pass `xπ = zeros(2N)`

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")`

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

```

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [December 10, 2020, 3:06pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/4 "2020-12-10T15:06:07Z")

</div>

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

`julia> ]dev BifurcationKit`

```julia
γ = 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())

```

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 10, 2020, 4:14pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/5 "2020-12-10T16:14:32Z")

</div>

There seems to be a bug.

You can try

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

```

while I fix it.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 10, 2020, 4:18pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/6 "2020-12-10T16:18:51Z")

</div>

Actuellement, it was another constructor

```julia
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](https://github.com/rveltz/BifurcationKit.jl/issues/27)

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 11, 2020, 7:59pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/7 "2020-12-11T19:59:57Z")

</div>

Did you get it to work?

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [December 17, 2020, 8:58pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/8 "2020-12-17T20:58:43Z")

</div>

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.

 ![Screen Shot 2020-12-17 at 10.53.44 AM](https://global.discourse-cdn.com/julialang/original/3X/5/9/59fbad0276208b6424c9f98161eb0ee356e16af7.jpeg)

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

```julia

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"])

```

 ![Screen Shot 2020-12-17 at 3.54.00 PM](https://global.discourse-cdn.com/julialang/original/3X/c/c/cc2e90440362a6e8d0a63a7351f255c41a6d928d.jpeg)

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 17, 2020, 9:14pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/9 "2020-12-17T21:14:22Z")

</div>

> [@cadojo](#):
>
> `t, z, D, ω, c, k, α, A, ω`

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:

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

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

```

To get the solution, it is easier to call:

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

```

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

---

<div class="post-metadata">

**Author:** ![cadojo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cadojo/32/25328_2.png) [@cadojo](https://discourse.julialang.org/u/cadojo)\
**Post date:** [December 18, 2020, 6:15pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/10 "2020-12-18T18:15:34Z")

</div>

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

> [@rveltz](#):
>
> The guess is wrong (4th component) and it should be ordered in the other way

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!

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 18, 2020, 8:18pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/11 "2020-12-18T20:18:45Z")

</div>

The latter. I updated the [docs](https://rveltz.github.io/BifurcationKit.jl/dev/library/#BifurcationKit.PeriodicOrbitTrapProblem) accordingly. Please keep us posted.

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

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

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 19, 2020, 3:08pm UTC](https://discourse.julialang.org/t/numerical-continuation-and-periodic-orbits/51580/12 "2020-12-19T15:08:07Z")

</div>

```julia

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,:])

```