Collocation for finding periodic orbits of ODEs

Every now and again I’m asked how to compute the periodic orbits of ODEs using a boundary value solver. Each time, I go looking for old code that does this and, each time, I can’t find it and end up rewriting the collocation code from scratch.

This time I thought I’d put my code here so that I have a better chance of finding it again in the future!

The basic idea is to use a Fourier differentiation matrix to approximate the derivatives along the orbit and use a nonlinear solver to ensure that those derivatives match the vector field. If you want to know more about these types of spectral methods, take a look at the excellent (and short!) introduction by Trefethen in “Spectral Methods in MATLAB”, SIAM 2000. If you want more detail then the magnum opus by Boyd “Chebyshev and Fourier Spectral Methods”, Dover 2001 (freely available on his personal website) is also very good.

Nowadays, my preference is for coding in Julia – it’s very clean and flexible. Here is the code (which could be better!).

# Released under the MIT expat license by David A.W. Barton (david.barton@bristol.ac.uk) 2020
using StaticArrays
using NLsolve
using OrdinaryDiffEq

"""
    duffing(u, p, t)

The vector field of the forced Duffing equation.
"""
duffing(u, p, t) = SVector(u[2], p.Γ*sin(p.ω*t) - 2p.ξ*u[2] - p.ωₙ^2*u[1] - p.β*u[1]^3)

"""
    fourier_diff([T=Float64,] N; order=1)

Create a Fourier differentiation matrix of the specified order with numerical type T on the
domain `x = LinRange{T}(0, 2π, N+1)[1:end-1]`.
"""
function fourier_diff(T::Type{<:Number}, N::Integer; order=1)
    D = zeros(T, N, N)
    n1 = (N - 1) ÷ 2
    n2 = N ÷ 2
    x = LinRange{T}(0, π, N+1)
    if order == 1
        for i in 2:N
            sgn = (one(T)/2 - iseven(i))
            D[i, 1] = iseven(N) ? sgn*cot(x[i]) : sgn*csc(x[i])
        end
    elseif order == 2
        D[1, 1] = iseven(N) ? -N^2*one(T)/12 - one(T)/6 : -N^2*one(T)/12 + one(T)/12
        for i in 2:N
            sgn = -(one(T)/2 - iseven(i))
            D[i, 1] = iseven(N) ? sgn*csc(x[i]).^2 : sgn*cot(x[i])*csc(x[i])
        end
    else
        error("Not implemented")
    end
    for j in 2:N
        D[1, j] = D[N, j-1]
        D[2:N, j] .= D[1:N-1, j-1]
    end
    return D
end
fourier_diff(N::Integer; kwargs...) = fourier_diff(Float64, N; kwargs...)

"""
    collocation_setup(u)

Return a data structure used internally by the `collocation!` function. `u` should be a
matrix with states down the columns and time across the rows (used for size/type information
only).
"""
function collocation_setup(u::AbstractMatrix)
    return (ndim=size(u, 1), nmesh=size(u, 2), Dt=-fourier_diff(eltype(u), size(u, 2))*2π)
end

"""
    collocation!(res, f, u, p, T, coll)

Calculate the residual of the collocation equations using a Fourier discretisation. Assumes
that a phase condition is not required (i.e., the equations are non-autonomous or the period
is known).

# Arguments
- `res`: residual (mutated)
- `f`: vector field function (expected to take the arguments (u, p, t))
- `u`: state variables along the orbit (vector)
- `p`: parameter vector passed to the vector field function
- `T`: period of oscillation
- `coll`: the output of `collocation_setup`

# Returns
- `res`: residual
"""
function collocation!(res, f, u, p, T, coll)
    # Matrix of derivatives along the orbit
    D = reshape(u, (coll.ndim, coll.nmesh))*coll.Dt
    ii = 1:coll.ndim
    for i in 1:coll.nmesh
        # Subtract the desired derivative from the actual derivative
        res[ii] .= D[ii] .- T.*f(u[ii], p, T*(i-1)/coll.nmesh)
        ii = ii .+ coll.ndim
    end
    return res
end

function example(; nmesh=20)
    p = (Γ=0.1, ω=1.0, ξ=0.05, ωₙ=1.0, β=0.1)
    # Do initial value simulation to get a reasonable starting point
    prob = ODEProblem(duffing, SVector(0.0, 0.0), (0.0, 100*2π/p.ω), p)
    odesol = solve(prob, Tsit5())
    # Refine using collocation
    t = range(0, 2π/p.ω, length=nmesh+1)[1:end-1]
    uvec = reinterpret(Float64, odesol(99*2π/p.ω .+ t).u)
    umat = reshape(uvec, (:, nmesh))
    coll = collocation_setup(umat)
    nlsol1 = nlsolve((res, u) -> collocation!(res, duffing, u, p, 2π/p.ω, coll), uvec)
    # Adjust the parameters slightly (actually quite a bit!) and correct
    p = (Γ=0.1, ω=1.1, ξ=0.05, ωₙ=1.0, β=0.1)
    nlsol2 = nlsolve((res, u) -> collocation!(res, duffing, u, p, 2π/p.ω, coll), uvec)
    return (nlsol1, nlsol2)
end

function plot_example()
    # Needs `using Plots` or similar
    nmesh = 20
    # The two solutions don't actually have the same period but normalize to [0, 2π]
    t = linspace(0, 2π, length=nmesh+1)[1:end-1]
    (sol1, sol2) = example()
    plot(t, sol1.zero[1:2:end])
    plot!(t, sol2.zero[1:2:end])
end

If you insist on using Matlab, the translation of the Julia code is below. Note that this uses the fourdif function by Reddy and Weideman to generate the Fourier differentiation matrix. (Also note that this can be put in a single file called fourier_collocation.m.)

function [nlsol1, nlsol2] = fourier_collocation()
% FOURIER_COLLOCATION Implement Fourier collocation for an arbitrary autonomous
% ODE. Assumes that the equations are non-autonomous or the period is known.

% Released under the MIT expat license by David A.W. Barton (david.barton@bristol.ac.uk) 2020

nmesh = 20;
p = struct('Gamma', 0.1, 'omega', 1.0, 'xi', 0.05, 'omegan', 1.0, 'beta', 0.1);
% Do initial value simulation to get a reasonable starting point
sol = ode45(@(t, u)duffing(t, u, p), [0, 100*2*pi/p.omega], [0, 0]);
% Refine using collocation
t = linspace(0, 2*pi/p.omega, nmesh+1);
t = t(1:end-1);
umat = deval(sol, 99*2*pi/p.omega + t);
uvec = umat(:);
coll = collocation_setup(umat);
nlsol1 = fsolve(@(u)collocation(@duffing, u, p, 2*pi/p.omega, coll), uvec)
% Adjust the parameters slightly (actually quite a bit!) and correct
p = struct('Gamma', 0.1, 'omega', 1.1, 'xi', 0.05, 'omegan', 1.0, 'beta', 0.1);
nlsol2 = fsolve(@(u)collocation(@duffing, u, p, 2*pi/p.omega, coll), uvec)

plot(t, nlsol1(1:2:end), 'b', t, nlsol2(1:2:end), 'r');

end

function du = duffing(t, u, p)
    du = [u(2); p.Gamma*sin(p.omega*t) - 2*p.xi*u(2) - p.omegan^2*u(1) - p.beta*u(1)^3];
end

function coll = collocation_setup(u)
    [~, D] = fourdif(size(u, 2), 1);
    coll = struct('ndim', size(u, 1), 'nmesh', size(u, 2), 'Dt', -D*2*pi);
end

function res = collocation(f, u, p, T, coll)
    % Matrix of derivatives along the orbit
    res = zeros(size(u));
    D = reshape(u, [coll.ndim, coll.nmesh])*coll.Dt;
    ii = 1:coll.ndim;
    for i = 1:coll.nmesh
        % Subtract the desired derivative from the actual derivative
        res(ii) = D(ii) - T*f(T*(i-1)/coll.nmesh, u(ii), p)';
        ii = ii + coll.ndim;
    end
end

Research Associate in Dynamics and Uncertainty

We are looking to recruit a passionate researcher to join a multi-institution project as a Research Associate, tenable for 1 year with the potential of an extension subject to satisfactory progress and funding availability.

The position is funded by an EPSRC-funded Programme Grant on Digital Twins for Improved Dynamic Design run in collaboration with the universities of Cambridge, Liverpool, Sheffield, Southampton, and Swansea. The overall aim of the programme grant is to create a robustly-validated virtual prediction tool called a “digital twin” for designing complex structures.

The Bristol component of this project is focussed on the interaction of numerical models and physical experiments in so-called hybrid tests where two substructures (one physical and one numerical) are coupled together in real-time. This novel approach to testing provides a great deal of flexibility in the design process since it enables poorly modelled (or externally supplied) parts of the structure to be tested physically while retaining the freedom to rapidly change and test different numerical models. Of particular interest within this testing framework is the investigation and exploitation of nonlinear behaviours and how uncertainties manifest and propagate through the system.

This position is ideal for a researcher with an interest in uncertainty quantification, nonlinear behaviour, and control. You should have good computational skills (MATLAB and/or Julia are commonly used) and experience of working with experimental data. You will be part of a vibrant Dynamics and Control research group, working alongside researchers dealing wide variety of application areas. In addition, there will be regular meetings of the full project team from each of the universities with an emphasis on cross-fertilisation of ideas and collaborative working.

Apply via the University of Bristol online portal.

PhD position in dynamical systems/nonlinear dynamics and Julia

A bit of a long shot, but if there is anyone who is looking to do a PhD in dynamical systems/nonlinear dynamics, would like to develop Julia-based dynamics software, and is based in the UK, I’d love to hear from you. I’m particularly interested in stochastic dynamics and links to machine learning.

I’m part of a small group researchers in dynamics at the University of Bristol, UK based in the Department of Engineering Mathematics and there is the opportunity for funded PhD studentships (ca. £15k a year plus tuition fees) starting in January. These are competitively awarded (i.e., I’m not guaranteed any for this project) and unfortunately restricted to people who have been resident in the UK for a minimum of 3 years and have leave to stay (see https://epsrc.ukri.org/skills/students/help/eligibility/ – note that the 10% rule mentioned has already been allocated this year, so you do need to be UK-based).

The deadline for application is the end of September 2019. (It’s not an entirely strict deadline but there are some internal processes I will need to complete before the hard deadline.)

Post-doc position available

I have a post-doctoral research associate position available to work on control-based continuation (nonlinear dynamics in experiments). The position will run until May 2020 with a possible extension to August 2020 (subject to EPSRC approval).

For details see the University of Bristol jobs website. The deadline for applications is 24 February 2019 with a provisional interview date of 7 March.

The text of the advert is below –

We seek a highly motivated Research Associate who is interested in working as part of a team at the interface between Engineering and Applied Mathematics to investigate new methods for exploring the nonlinear behaviour of engineered systems. The post will run until 31 May 2020, funded by an EPSRC grant with the possibility of an extension subject to funds and EPSRC permission.

Modern test methods for investigating the dynamics of engineered structures are inadequate for dealing with the presence of significant nonlinearity since they have largely been developed under the assumption of linear behaviour. In contrast, control-based continuation (CBC), a versatile non-parametric identification method, has been developed with nonlinearity in mind from the beginning. It has been demonstrated on simple experiments but now advances in underlying methodology are required to apply CBC to real-world experiments which have higher levels of measurement noise and many degrees of freedom. The versatility of CBC is such that, with these advances, it will also become relevant for researchers studying nonlinear systems in both engineering and other fields, such as in the biological sciences.

We are seeking a Research Associate to drive this research forward alongside other researchers (both PhD students and other post-doctoral staff) who are working on closely related problems. Support will be readily available from the investigators David Barton, Simon Neild and Djamel Rezgui. More widely, you will be part of the Dynamics and Control research group and the Applied Nonlinear Mathematics research group both of which carry out cutting-edge research in a wide range of application areas.

CBC presently draws on a wide range of underlying areas including, but not limited to, dynamical systems and bifurcation theory, control theory, system identification, and machine learning. Applicants are expected to have experience in at least one of these areas in addition to a first degree and preferably a PhD in Applied Mathematics/Physics/Engineering (or a closely related discipline).

Possible initial avenues of research include

  • Improving the robustness of CBC in the presence of noise using surrogate models. Gaussian processes have previously been investigated and may be useful.
  • Investigating the scaling up of CBC to many degree-of-freedom systems. Ideas from numerical continuation of PDE systems could yield insights.
  • Implementation of CBC on existing aerospace experiments for dynamic testing and wind tunnel testing.

New PhD scholarship opportunity in robotics/machine learning

There is the opportunity for fully-funded PhD scholarships starting September 2019 as part of the next University of Bristol funding competition. The deadline for applications is January 2019 (the precise date is to be announced).

Funding can be awarded to students of any nationality, though the chances of funding are likely higher for UK nationals (and others eligible for EPSRC doctoral funding) and Chinese nationals (via the CSC funding programme) since more funding is available through those routes.

I am particular interested in recruiting students for a PhD opportunity in tactile robotics and machine learning (though do also get in touch if you are considering nonlinear dynamics and control more generally).

A short project description is below.

Present approaches to tactile sensing and control require large amounts of data to train machine learning algorithms, or other statistical methods, to transform low-level sensory data into high-level information such as contact position, angle and force. Once a suitable model is learnt from data it is then used to within a control policy to complete the desired robotic manipulation task. While this approach is effective, it is far from efficient. This project will investigate the use of online learning combined with a high-level objective function to minimise the amount of prior training required. A local interaction model can be learnt from online sensor readings and the known movements between them and, as such, a robot manipulator can learn how to interact with its surroundings as it is carrying out useful tasks. This project has the opportunity to make use of extensive experimental facilities in conjunction with the Bristol Robotics Laboratory.

A more detailed version is also available.

Funding available for PhD positions – Oct 2018 start

There is funding available (competitively awarded across my department) for PhD students to start September/October 2018. There are a variety of funding sources including: EPSRC, China Scholarship Council (CSC), and University Scholarships. These all provide funding for fees and living costs.

I am particularly interested in topics around computational dynamics with links to machine learning and/or uncertainty quantification, largely from an engineering point of view but other areas might be considered.

I’m also interested in experiment-based dynamics and the real-time link with computational dynamics.

If you are considering a PhD in any of these areas, get in touch with me at david.barton@bristol.ac.uk.

(The image above is borrowed from Mike Henderson’s Multifaro page – a nice example of computational dynamics in action!)

Postdoc position now available

3 year postdoc position available for 1 June 2017 start!

We seek a highly motivated Research Associate who is interested in working as part of a team at the interface between Engineering and Applied Mathematics to investigate new methods for exploring the nonlinear behaviour of engineered systems and to develop numerical continuation techniques for physical experiments.

Modern test methods for investigating the dynamics of engineered structures are inadequate for dealing with the presence of significant nonlinearity since they have largely been developed under the assumption of linear behaviour. In contrast, control-based continuation (CBC), a versatile non-parametric identification method, has been developed with nonlinearity in mind from the beginning. It has been demonstrated on simple experiments but now advances in underlying methodology are required to apply CBC to real-world experiments which have higher levels of measurement noise and many degrees of freedom. The versatility of CBC is such that, with these advances, it will also become relevant for researchers studying nonlinear systems in both engineering and other fields, such as in the biological sciences.

We are seeking a Research Associate to drive this research forward alongside researchers working on closely related problems from the Departments of Engineering Mathematics, Mechanical Engineering and Aerospace Engineering. Support will be readily available from the investigators David Barton, Simon Neild and Djamel Rezgui. More widely, you will be part of the Dynamics and Control research group and the Applied Nonlinear Mathematics research group both of which carry out cutting-edge research in a wide range of application areas.

CBC presently draws on a wide range of underlying areas including, but not limited to, dynamical systems and bifurcation theory, system identification, control theory and machine learning. Applicants are expected to have experience in at least one of these areas in addition to a first degree and preferably a PhD in Applied Mathematics/Physics/Engineering (or a closely related discipline).

Possible initial avenues of research include

  • Improving the robustness of CBC in the presence of noise using surrogate models. Gaussian processes have previously been investigated and may be useful.
  • Investigating the scaling up of CBC to many degree-of-freedom systems. Ideas from numerical continuation of PDE systems could yield insights.
  • Implementation of CBC on existing aerospace experiments for dynamic testing and wind tunnel testing.

The post is available from 1 June 2017 with funding available for up to three years.

To apply visit http://www.bristol.ac.uk/jobs/find/details.html?nPostingId=5644&nPostingTargetId=21245

Please direct any questions to David Barton, david.barton@bristol.ac.uk.

Studentships available for immediate start

We have some EPSRC funded studentships available for immediate start (competitively awarded) in my department/school. The studentships pay an annual stipend of ~£14k plus tuition fees, though you must meet EPSRC eligibility requirements (there are no awards available for overseas students).

The deadline is very tight though, you need to have an application completed for 10 October.

Applications coming in after this time will be considered for a September 2017 start; but please do still get in touch.

If you are interested in nonlinear dynamics (from a numerical/computational/experimental viewpoint), please get in touch.  I’m particularly interested in stochastic dynamics and/or links with machine learning.

Email me at david.barton@bristol.ac.uk with your current CV and a short description of your (academic) interests. General enquires along the lines of “I don’t really know what I want to do” are also welcome.

Nonlinear vibrations, localisation and energy transfer

I recently had the privilege of being invited to the 6th conference on Nonlinear Vibrations, Localization and Energy Transfer in Liege, 4-8 July 2016. The conference was superbly organised by Gaetan Kerschen and Jean-Philippe Noël in remarkable surroundings and with very good food (and beer — a trip to the Val-Dieu abbey to taste good Belgian beer was a highlight of the social part of the conference!).

The conference was a good mix of cutting-edge research combined with tutorial lectures for PhD students (and above!) on a range of topics. The lectures on nonlinear system identification were particularly useful for me — it’s an area that I’m keen to learn more about and import into my own work.

I had the pleasure of presenting on “Control-based Continuation – A new experimental approach for testing nonlinear dynamic systems” (with co-authors Ludovic Renson and Simon Neild); hopefully I was able to impart some of the enthusiasm we have for working on nonlinear experiments and understanding their nonlinear dynamics and bifurcations via numerical continuation!

Other talks of particular interest were by Simon Peter (Stuttgart) and Peter Bruns (Hannover) on using phase-locked loops to extract backbone (nonlinear normal modes) from experiments. Their approaches are highly complementary to my own. Cyril Touzé also gave a fascinating talk on turbulence in metallic plates, along with convincing experiments that show energy transfers from low-frequency modes to higher-frequency modes of vibration. I was also pleased to see that software developments are slowly propagating through into applied communities; Bruno Cochelin presented an interesting take on using Automatic Differentiation for high-order numerical continuation methods, I’m curious as to how far this can be applied.

Overall a thoroughly productive week!