Continuation of nonlinear normal modes with BifurcationKit.jl

[KaioBenedetti]

Dear all,

I am new to Julia, and came to notice this package for the continuation of autonomous dynamical systems. I tried to apply it to the computation of normal modes of second-order dynamical systems. Briefly, I am interested in the free vibrations of nonlinear mechanical systems, which are modeled as second-order autonomous ODES. However, these vibrations can be quite complex due to their nonlinearity. Specifically, there are solutions with Auto that I wish to mimic:
COMPUTATION OF NONLINEAR NORMAL MODES, PART II: NUMERICAL CONTINUATION IN AUTO
I have already tried following the examples in the package, with the trapezoidal method and simple shooting, but with no success.

Edit:
As an example, consider the Duffing oscillator, which is a second-order autonomous conservative ODE:

where the H function is the Hamiltonian (total energy). The relation Amplitude vs frequency of vibration is called the backbone curve. If b = 0, a > 0, then the system is a linear oscillator and the backbone is a vertical line at omega = sqrt(a). For any other real parameters, this line bends towards higher or lower omegas. This is the curve I am interested in.

The problem is, for constant a, b, the continuation problem does not depend on any parameter, but on the initial energy itself, i.e., the initial condition. If I encode this system, this would result in:

> using BifurcationKit, Setfield, Parameters
> function F!(f, x, p, t=0)
    @unpack a, b = p
    f[1] = x[2]
    f[2] = -a*x[1] -b*(x[1]^3)) 
    return f
end
F(x,p,t = 0) = F!(similar(x),x,p,t)

Here, F is the vector field of the first-order equivalent system, with phase-space variables x[1] = x, x[2] = dx/dt. If i follow the provided reference, then this vector field would be rewritten as

> function F2!(f, x, p, t=0)
@unpack a, b, T, λ = p
    f[1] = T*(x[2]) + λ*(a*x[1] + b*(x[1]^3))
    f[2] = T*(-a*x[1] -b*(x[1]^3)) + λ*(x[2])
    return f
end
F2(x,p,t = 0) = F2!(similar(x),x,p,t)

and now the system has two additional parameters, the period T and the unfolding parameter lambda. F2 should be continued against lambda by a collocation method, similar to how Auto proceeds. There are at least 2 branches for F2, one for trivial solutions, x[1,2] = 0, and another bifurcated from lambda = 0. This bifurcated branch is the backbone of the original system.

I tried following the tutorials of BifurcationKit.jl for F2, applying the trapezoidal continuation method and shooting, with no success. The original F has no explicitly continuation parameters, so I do not know if the package can deal with it directly. Kerschen also defines an algorithm for the computation of backbones directly, based on shooting and a Moore-Penrose inverse:

Peeters M, Viguié R, Sérandour G, Kerschen G, Golinval JC: Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques. Mech. Syst. Signal Process. 2009, 23 (1) :195–216. Redirecting

Any help would be deeply appreciated!
Kaio Benedetti

[cmarcotte]

Hi Kaio! It might be easier for people with experience using BifurcationKit.jl to help you if you provide a minimal example of what you are trying to accomplish so that we can diagnose what is not working.

[KaioBenedetti]

Thanks for the reply! I have edited the original post, explaining the problem in more detail and containing the vector field of a typical example.

[cmarcotte]

I can certainly get the trivial solution using BifurcationKit based on what you’ve described, and even the bifurcation at \lambda=0, but without information about your starting parameters and starting solution I’m not sure what more help I can offer. If you include a complete, runnable, example of the script you are working on, then perhaps I can understand better. My own script is below, based on the pp2 example.

using Revise, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

# define the sup norm
norminf(x) = norm(x, Inf)

# function to record information from a solution
recordFromSolution(x, p) = (u1 = x[1], u2 = x[2])

####################################################################################################

function F!(f, x, p, t=0)
	@unpack a, b, T, λ = p
    f[1] = T*(x[2]) + λ*(a*x[1] + b*(x[1]^3))
    f[2] = T*(-a*x[1] -b*(x[1]^3)) + λ*(x[2])
    return f
end
F(x,p) = F!(similar(x),x,p,0.0)

jet  = BK.getJet(F; matrixfree=false)

# parameters of the model
par_F = (a = 1., b = 0.5, T = 1., λ = -0.01)

# initial condition
z0 = zeros(2)

# continuation options
opts_br = ContinuationPar(pMin = -0.1, pMax = 1.0, dsmax = 0.01,
	# options to detect bifurcations
	detectBifurcation = 3, nInversion = 8, maxBisectionSteps = 25,
	# number of eigenvalues
	nev = 2,
	# maximal number of continuation steps
	maxSteps = 1000,
	# parameter theta, see `? continuation`. Setting this to a non
	# default value helps passing the transcritical bifurcation
	theta = 0.3)
	
diagram = bifurcationdiagram(jet...,
	# initial point and parameter
	z0, par_F,
	# specify the continuation parameter
	(@lens _.λ),
	# very important parameter. This specifies the maximum amount of recursion
	# when computing the bifurcation diagram. It means we allow computing branches of branches of branches
	# at most in the present case.
	3,
	(args...) -> setproperties(opts_br; ds = 0.001, dsmax = 0.01, nInversion = 8, detectBifurcation = 3);
	# δp = -0.01,
	recordFromSolution = recordFromSolution,
	verbosity = 1, plot = true)

scene = plot(diagram; code = (), title="$(size(diagram)) branches", legend = false)

But things go off the rails when I try to select the Hopf point:
brH = BK.getBranch(diagram, (n,m)).γ fails, for all n,m \in \{0,1,2,3\} (though I don’t know where else the point is in the diagram structure). Perhaps a question for @rveltz?
By diagram structure is below:

julia> diagram
[Bifurcation diagram]
 ┌─ From 0-th bifurcation point.
 ├─ Children number: 0
 └─ Root (recursion level 1)
      ┌─ Branch number of points: 90
      ├─ Branch of Equilibrium
      ├─ Type of vectors: Vector{Float64}
      ├─ Parameter λ starts at -0.01, ends at 1.0
      └─ Special points:

 (ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`)

- #  1,  hopf at λ ≈ +0.00000022 ∈ (-0.00000033, +0.00000022), |δp|=6e-07, [converged], δ = ( 2,  2), step =   5, eigenelements in eig[  6], ind_ev =   2

Sorry I could not be more helpful!

[KaioBenedetti]

Many thanks again for your time! I have constructed a code example from yours:

using Revise, Parameters, Setfield, Plots, LinearAlgebra
using BifurcationKit
const BK = BifurcationKit

# define the sup norm
norminf(x) = norm(x, Inf)

# function to record information from a solution
recordFromSolution(x, p) = (u1 = x[1], u2 = x[2])

####################################################################################################

function F!(f, x, p, t=0)
	@unpack a, b, T, λ = p
    f[1] = T*(x[2]) + λ*(a*x[1] + b*(x[1]^3))
    f[2] = T*(-a*x[1] -b*(x[1]^3)) + λ*(x[2])
    return f
end
F(x,p) = F!(similar(x),x,p,0.0)

# Defining the Jacobian
function dF(x,p)
	@unpack a, b, T, λ = p
    M = zeros(2,2)
    M[1,1] = T*0 + λ*(a + 3*b*(x[1]^2))
    M[1,2] = T*1 + λ*0
    M[2,1] = T*(-a -b*3*(x[1]^2)) + λ*0
    M[2,2] = T*0 + λ*1
    return M
 end

jet  = BK.getJet(F; matrixfree=false)

# parameters of the model
P = @with_kw (a=1., b=0.5, T=2*π, λ=-0.01)
par_F = P()

# initial condition
z0 = [0.001;0.001]

# continuation options
opt_newton = NewtonPar(verbose = true, tol = 1e-9, maxIter = 25)
opts_br = ContinuationPar(pMin = -0.1, pMax = 1.0, dsmax = 1., ds = 0.001,
	# options to detect bifurcations
	detectBifurcation = 3, nInversion = 8, maxBisectionSteps = 25,
	# number of eigenvalues
	nev = 2,
	# maximal number of continuation steps
	maxSteps = 1000,
	# parameters for the Newton method
	newtonOptions = opt_newton,
	# parameter theta, see `? continuation`. Setting this to a non
	# default value helps passing the transcritical bifurcation
	theta = 0.3)

# Computing the trivial branch (equilibrium points)
br, = continuation(F, dF,	# here I don't known how to use the jet variable
 	z0, par_F,
	(@lens _.λ),
	opts_br, verbosity = 1,
	recordFromSolution = recordFromSolution,
	plot = true,
	normC = norminf)

# index of the Hopf point in br.specialpoint
ind_hopf = 1

# number of time slices
M = 51

# Hopf bifurcation orbit guess
p_hopf, Th, orbitguess, hopfpt, vec_hopf = guessFromHopf(
	# branch and index of bifurcation point
	br, ind_hopf,
	# Newton eigensolver
	opts_br.newtonOptions.eigsolver,
	# number ot time slices and amplitude guess
	M, 1)

# reshaping orbit initial guess as a vector
orbitguess = reduce(vcat, orbitguess)
orbitguess = vcat(vec(orbitguess), Th) |> vec
plot(orbitguess[1:2:end-1])
plot!(orbitguess[2:2:end-1])

# defining the periodic orbit trapezoidal problem
poTrap = PeriodicOrbitTrapProblem(
	F,    				# pass the vector field
	dF, 				# pass the jacobian of the vector field
	real.(vec_hopf),	# used to set ϕ, see the phase constraint
	hopfpt.u,           # used to set uhopf, see the phase constraint
	M,					# number of time slices
	2)			        # phase-space dimension

# arguments for period orbit continuation
args_po = (	recordFromSolution = (x, p) -> begin
		xtt = BK.getPeriodicOrbit(p.prob, x, @set par_F.T = p.p)
		return (max = maximum(xtt[1,:]),
				min = minimum(xtt[1,:]),
				period = getPeriod(p.prob, x, @set par_F.T = p.p))
	end,
	normC = norminf)

# continuation parameters for periodic orbits
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.03, ds= 0.01,
	pMax = 30.0, maxSteps = 100,
	newtonOptions = opt_newton, nev = 2, precisionStability = 1e-8, detectBifurcation = 3)

# here I am trying to continuate from the hopf point to the backbone, controling parameter T
br_po, = continuation(poTrap,
	orbitguess, (@set par_F.T = par_F.T + 0.01), (@lens _.T),
	opts_po_cont;
	verbosity = 2,	plot = true,
	args_po)

I have change some steps, namely the initial period definition, which should be 2*pi for a = 1, and the rest of the code is my attempt using the trapezoidal method. The initial continuation computes the trivial branch, and the trapezoidal method should have computed the backbone. However, I am not sure about the second procedure, if it indeed is the correct way. The authors of the first reference told that they indeed continuated in the parameter T, which is what I am trying.

[rveltz]

Hi,

I looked at the paper you refer. If you speak about (9a-9b), it seems that the problem is more difficult than what you state. You look for periodic functions solutions of (9a-9b). Auto07p has a dedicated structure for that and we dont have it in BifurcationKit.jl

You could use a BifurcationKit.PeriodicOrbitTrapProblem but it is tailored to solving perioidic solutions to dx/dt = F(x). I think (9a-9b) is more general than that unless I am wrong. In case it is more general than that, you can probably sett up a TwoPointBVProblem from DifferentialEquations.jl and follow the solutions with BifurcationKit.

The second paper using Shooting is in the same vain. It is different from what BK can do. It can of course find the solutions (z,T) by shooting but it cannot parametrize it as function of the variable T itself. Indeed, the period is an unknown, not a parameter.

The easiest way to proceed would be to copy the structure ShootingProblem and pass the period as a parameter like

function (sh::MyShootingPb)(x::AbstractVector, par)
	T = par.period
	M = getM(sh)
	N = div(length(x) - 1, M)

	# extract the orbit guess and reshape it into a matrix as it's more convenient to handle it
	xc = getTimeSlices(sh, x)

	# variable to hold the computed result
	out = similar(x)
	outc = getTimeSlices(sh, out)
		for ii in 1:M
			ip1 = (ii == M) ? 1 : ii+1
			# we can use views but Sundials will complain
			outc[:, ii] .= sh.flow(xc[:, ii], par, sh.ds[ii] * T) .- xc[:, ip1]
		end

	# add constraint
	out[end] = @views sh.section(getTimeSlice(sh, xc, 1), T)

	return out
end

[KaioBenedetti]

Thanks for the reply!

I see. I thought it should be possible to apply the trapezoidal method or even the collocation method to continue these solutions because, in the end, it is equation (3) with periodicity condition (4) that is solved. Indeed, for such problems coupled with a phase-condition (x[2] = 0, in this example), the solution is continued in the product space RxT, where R is the phase-space and T is the period space.
Nevertheless, is it possible to pass a TwoPointBVProblem to the continuation procedures? I am thinking of defining a shooting procedure by hand for those problems. If that is the case, it should be possible to apply the tools implemented in BifurcationKit.jl, right?

[rveltz]

@KaioBenedetti

Nevertheless, is it possible to pass a TwoPointBVProblem to the continuation procedures?

It should be. It is just a functional to be evaluated and you can surely AD it to get the jacobian.

I am thinking of defining a shooting procedure by hand for those problems. If that is the case, it should be possible to apply the tools implemented in BifurcationKit.jl , right?

I don’t know. Your periodic orbits are not isolated, there are foliated by the energy.

[KaioBenedetti]

Thanks for the response.

Your periodic orbits are not isolated, there are foliated by the energy.

What do you mean by this? I know that the periodic orbits are not uniquely defined, thus the continuation must be supplemented by a phase condition.