GitHub

Invariant Models

This package is an improved version of FMA. Current capabilities are

  • Identifies invariant foliations in autonomous, periodic and quasi-periodic systems from data.
  • Calculates invariant manifolds in autonomous, periodic and quasi-periodic systems from discrete-time systems and vector fields.
  • Calculates invariant foliations in autonomous, periodic and quasi-periodic systems from discrete-time systems and vector fields.
  • Calculates instantaneous frequencies and damping ratios on two-dimensional invariant manifolds.

Background

The system in question is assumed to obey the discrete-time map

\[\begin{aligned} \boldsymbol{x}_{k+1} &=\boldsymbol{F}\left(\boldsymbol{x}_{k},\boldsymbol{\theta}_{k}\right)\\\boldsymbol{\theta}_{k+1} &= \boldsymbol{\theta}_{k}+\boldsymbol{\omega} \Delta t, \end{aligned} \tag{MAP}\]

where it is assumed that $\Delta t$ is the sampling period, that is the duration of time between two samples $\boldsymbol{x}_k$ and $\boldsymbol{x}_{k+1}$. or the vector field

\[\begin{aligned} \dot{\boldsymbol{x}} &=\boldsymbol{F}\left(\boldsymbol{x},\boldsymbol{\theta}\right)\\\dot{\boldsymbol{\theta}} &= \boldsymbol{\omega}, \end{aligned} \tag{ODE}\]

It may be the case that variable $\boldsymbol{\theta}$ is absent and therefore the system is autonomous.

Invariance

There are four ways to connect a low-order model $\boldsymbol{R}$ to $\boldsymbol{F}$. The figure below shows the four combinations. Only invariant foliations and invariant manifolds produce meaningful reduced order models. Only invariant foliations and autoencoders can be fitted to data. The intersection is invariant foliations.

Therefore,

  • when a system of equations is given, invariant manifolds are the most appropriate (foliations are still possible),
  • when data is given, only invariant foliations are appropriate.

Autoencoders such as SSMLearn do not enforce invariance and therefore generate spurious results as shown in Data-Driven Reduced Order Models Using Invariant Foliations, Manifolds and Autoencoders.

Invariant Foliations

The invariance equation for foliations is

\[\boldsymbol{R}\left(\boldsymbol{U}\left(\boldsymbol{x},\boldsymbol{\theta}\right),\boldsymbol{\theta}\right) = \boldsymbol{U}\left(\boldsymbol{F}\left(\boldsymbol{x},\boldsymbol{\theta}\right),\boldsymbol{\theta}+\boldsymbol{\omega}\right) \tag{FOIL}\]

Invariant Manifolds

The invariance equation for manifolds is

\[\boldsymbol{S}\left(\boldsymbol{V}\left(\boldsymbol{x},\boldsymbol{\theta}\right),\boldsymbol{\theta}\right) = \boldsymbol{V}\left(\boldsymbol{F}\left(\boldsymbol{x},\boldsymbol{\theta}\right),\boldsymbol{\theta}+\boldsymbol{\omega}\right) \tag{MAN}\]

Case Study of a Minimal Forced System

We consider a forced version of the Shaw-Pierre example, that originally appeared in [].

\[\begin{aligned} \dot{x}_1 &= x_3 \\ \dot{x}_2 &= x_4 \\ \dot{x}_3 &= - c x_3 - k x_1 - \kappa y_1^3 + k (x_2 - x_1) + c (x_4 - x_3) + p \cos(\theta+0.1) \\ \dot{x}_4 &= - c x_4 - k x_2 - k (x_2 - x_1) - c (x_4 - x_3) + p \cos(\theta) \\ \dot{\theta} &= \omega \\ \end{aligned}\]

The parameters are $k=1$, $\kappa = 0.2$, $c = 2^{-5}$. The forcing frequency is $\omega = 1.2$ and the the forcing amplitude is $p=0.25$.

Without forcing the natural frequencies are approximately $\omega_1 = 1$, $\omega_2 = \sqrt{3}$. The damping ratios are $\zeta_1 = 0.01561$ and $\zeta_2 = 0.02705$. The spectral quotient for the first mode is $\beth_1 = 3$ and for the second mode is $\beth_2 = 1$. When forcing is turned on these values change slightly due to nonlinearity.

Set-up

First we import the required packages

using InvariantModels
using CairoMakie

Consider the Shaw-Pierre vector field

NDIM = 4

function shawpierre!(x, y, p, t)
    k = 1.0
    kappa = 0.2
    c = 2^-5

    x[1] = y[3]
    x[2] = y[4]
    x[3] = - c*y[3] - k*y[1] - kappa*y[1]^3 + k*(y[2]-y[1]) + c*(y[4]-y[3]) + p * cos(t+0.1)
    x[4] = - c*y[4] - k*y[2] - k*(y[2]-y[1]) - c*(y[4]-y[3]) + p * cos(t)
    return x
end

function shawpierre(z, p, t)
    dz = zero(z)
    return shawpierre!(dz, z, p, t)
end
shawpierre (generic function with 1 method)

Set up some system parameters

Amplitude = 0.25        # forcing amplitude
omega_ode = 1.2         # forcing frequency

# parameters of the numerical method
fourier_order = 7       # fourier harmonics to be resolved
ODE_order = 7           # polynomial order of the calculations
SEL = [1 2]             # which invariant vector bundle to use
dispMaxAmp = 1.0        # maximum amplitude to display
1.0

Invariant Manifolds of Vector Fields

Create a polynomial MPode, XPode out of our vector field

MPode = QPPolynomial(NDIM, NDIM, fourier_order, 0, ODE_order)
XPode = fromFunction(MPode, (x, t) -> shawpierre(x, Amplitude, t))
4×330×15 Array{Float64, 3}:
[:, :, 1] =
 0.0        0.0       1.0       0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0        1.0       0.0       0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.248751   0.03125  -0.0625    1.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.25      -0.0625    0.03125  -2.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

[:, :, 2] =
 0.0        0.0       1.0       0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0        1.0       0.0       0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.217094   0.03125  -0.0625    1.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.228386  -0.0625    0.03125  -2.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

[:, :, 3] =
 0.0        0.0       1.0       0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0        1.0       0.0       0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.147899   0.03125  -0.0625    1.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.167283  -0.0625    0.03125  -2.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

;;; … 

[:, :, 13] =
 0.0         0.0       1.0       0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0         1.0       0.0       0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.100605    0.03125  -0.0625    1.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0772542  -0.0625    0.03125  -2.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

[:, :, 14] =
 0.0        0.0       1.0       0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0        1.0       0.0       0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.184995   0.03125  -0.0625    1.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.167283  -0.0625    0.03125  -2.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

[:, :, 15] =
 0.0        0.0       1.0       0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0        1.0       0.0       0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.237397   0.03125  -0.0625    1.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.228386  -0.0625    0.03125  -2.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0

The invariant manifold can now be calculated using QPODETorusManifold:

MK, XK, MSn, XSn, MWn, XWn, _, _, _, _, XWre, oeva = QPODETorusManifold(MPode, XPode, omega_ode, SEL, threshold=0.1, resonance=false)
(QPConstant{4, 7, ℝ}(Euclidean(4, 15; field=ℝ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), [-0.5851607327172489 -0.5032137404631323 … -0.446152015418224 -0.5648704815203789; -0.6041692786018199 -0.5146117408994948 … -0.472707892894138 -0.5893509885342703; 0.09125378349006542 0.3707150321377271 … -0.46382979776303934 -0.20563127531679973; 0.11021647846213264 0.39545080631049584 … -0.46470386262474167 -0.19388036864429362], InvariantModels.QPFourierPolynomial{2, 2, 7, 0, 7}([0 0 … 6 7; 0 1 … 1 0], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], Euclidean(2, 36, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; … ;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im], InvariantModels.QPFourierPolynomial{4, 2, 7, 0, 7}([0 0 … 6 7; 0 1 … 1 0], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], Euclidean(4, 36, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im -1.8376913549748457e-9 + 1.051594954146569e-9im … -1.5439941936356905e-12 + 9.248693889784205e-13im 2.1316402941299244e-13 - 2.607237887585272e-13im; 0.0 + 0.0im 1.932254633285949e-9 - 1.3080859575918551e-9im … 9.28598075803996e-13 - 7.233615874021984e-13im 7.232125210589336e-14 - 4.846538232505142e-13im; 0.0 + 0.0im -1.7694128944588512e-9 - 3.491877011805131e-9im … 2.4502228536544e-13 + 6.298273112677297e-13im -2.7468655152228985e-13 - 1.2067134099140594e-13im; 0.0 + 0.0im 2.1360102001730727e-9 + 3.442459872016412e-9im … 1.1387881679423064e-12 + 1.544295728981521e-12im -6.836996712822771e-13 - 1.496625619954442e-13im;;; 0.0 + 0.0im 1.0666824568035874e-8 + 1.0553868518826516e-8im … -1.7178021622532267e-5 + 4.196099527654553e-5im 1.9116224228103998e-5 + 3.148537663326983e-5im; 0.0 + 0.0im -2.5990556540578604e-10 - 2.0010935692654432e-10im … 8.435720043748689e-6 - 1.9456626081984035e-5im 9.523903200527766e-6 + 1.5589254396108572e-5im; 0.0 + 0.0im 8.728431277508471e-8 - 8.841851459625687e-8im … 8.64854412913753e-5 + 3.1757347853003205e-5im -1.3890440260468833e-6 - 3.880919251619603e-6im; 0.0 + 0.0im -9.379409229527458e-10 + 1.6619851978925422e-9im … -4.0207309281285285e-5 - 1.549498818746899e-5im -6.267234591758321e-7 - 1.9612164264237774e-6im;;; 0.0 + 0.0im -1.9772078347982622e-11 + 1.6744078229450918e-12im … 9.437730299039082e-13 - 7.765029331492843e-12im -4.610300696040772e-13 - 8.830437724455137e-13im; 0.0 + 0.0im 1.2977996962095417e-11 - 1.0142833758526422e-13im … 6.523221125096276e-13 - 7.469251232101085e-12im -1.2837213160291325e-13 - 1.535719619682079e-12im; 0.0 + 0.0im -7.166246536587105e-12 - 4.897808484575253e-11im … -7.589431799024196e-12 - 8.400183543664812e-13im 1.0778801958441185e-12 - 4.658009491791358e-13im; 0.0 + 0.0im -4.3604716720044987e-13 + 2.42051567975758e-11im … -7.79662003452865e-12 - 1.0305199858418843e-12im 1.8323445316441114e-12 + 3.8574508501966636e-14im;;; … ;;; 0.0 + 0.0im 1.8371897445235535e-11 + 5.344609012645188e-12im … 4.522739776233169e-13 - 5.868851997447412e-13im -5.802433581221375e-13 - 2.0123049344654438e-14im; 0.0 + 0.0im -1.12018256149926e-11 - 3.948741559273708e-12im … -6.534048607434035e-13 + 3.1896835749798125e-13im 5.553992959649155e-13 + 7.855372158807943e-14im; 0.0 + 0.0im 1.2415359019986196e-11 - 4.339243934295115e-11im … 5.483907444619001e-13 + 1.949715070906039e-12im -9.806575609203639e-14 + 5.803727421003136e-13im; 0.0 + 0.0im -6.412022755411212e-12 + 2.1392791784228224e-11im … 1.1935533985904015e-12 + 1.6059672558836076e-12im 9.759212319616923e-14 - 1.0427727223766795e-12im;;; 0.0 + 0.0im -4.0094552977912924e-8 + 2.4992752885677616e-8im … 1.99135390458561e-11 - 3.0303696861617204e-11im -4.891805743787622e-13 + 5.972302047859109e-13im; 0.0 + 0.0im 1.3107042754998609e-9 - 5.541178826697017e-10im … -4.358618592218882e-12 + 5.144782069119434e-12im 2.5210343637805417e-13 - 3.251492896919229e-13im; 0.0 + 0.0im -1.5423608299572943e-7 - 2.4839651701616223e-7im … 3.351993024295125e-10 + 2.0763615292777148e-10im -5.097471537290631e-12 - 3.949468075768805e-12im; 0.0 + 0.0im 2.855061474882248e-9 + 7.684832130582808e-9im … -9.74051170077051e-12 - 9.707888307085478e-12im 5.623430076274752e-13 + 4.880809204846179e-13im;;; 0.0 + 0.0im 1.6658031077518076e-9 + 1.028527422930394e-9im … -9.52226377470191e-13 - 2.8932690655971893e-13im 2.922659710814459e-14 + 8.890432921208027e-15im; 0.0 + 0.0im -1.7885297022407263e-9 - 1.2788788227361565e-9im … 9.26752907145415e-13 + 4.084780818668961e-13im 1.639664102834293e-14 + 1.19525279699603e-14im; 0.0 + 0.0im 1.6141528229650275e-9 - 2.6876239695338146e-9im … -3.437988289901927e-13 + 7.563112312846452e-13im -4.430977242938976e-15 - 1.4268074753611958e-15im; 0.0 + 0.0im -2.101211729627376e-9 + 3.198689745637149e-9im … 6.203731886133029e-13 - 1.6967305375035845e-12im -6.006597913402906e-15 + 3.1701784149461015e-14im], InvariantModels.QPFourierPolynomial{4, 4, 7, 0, 7}([0 0 … 6 7; 0 0 … 1 0; 0 0 … 0 0; 0 1 … 0 0], [5 4 3 2; 12 9 7 6; … ; 329 326 325 324; 330 329 328 327], [1 1 1 1; 2 2 2 2; … ; 209 209 209 209; 210 210 210 210], [1 1 1 1; 1 1 1 2; … ; 6 2 1 1; 7 1 1 1], Euclidean(4, 330, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; … ;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im], QPPolynomial{4, 4, 7, 1, 1, ℂ}([0 0 0 1; 0 0 1 0; 0 1 0 0; 1 0 0 0], [0 0 … 0 1; 0 0 … 1 0; 0 0 … 0 0; 0 1 … 0 0], [2, 3, 4, 5], [4 3 2 1], [1 1 1 1], [1 1 1 1], Array{Int64, 3}(undef, 0, 4, 4), Array{Int64, 3}(undef, 0, 4, 4), Array{Int64, 3}(undef, 0, 4, 4), Euclidean(4, 4, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.003023892931199202 - 0.34393748976902505im 0.003023892931199089 + 0.34393748976902555im -0.006364175040882533 + 0.48525095019578185im -0.006364175040882914 - 0.48525095019578113im; -0.011930478993420549 + 0.3527819212586303im -0.01193047899342003 - 0.3527819212586307im -0.012655352114915922 + 0.5306030873648908im -0.01265535211491539 - 0.5306030873648894im; -0.6450530663936729 - 0.0029433904442031896im -0.6450530663936727 + 0.0029433904442033565im 0.497514687772557 - 0.006621968681592415im 0.49751468777255725 + 0.006621968681592858im; 0.5912975949551774 - 0.005552120385582407im 0.5912975949551771 + 0.005552120385582606im 0.4836616902495777 - 0.009437763581131997im 0.48366169024957806 + 0.009437763581132406im;;; -0.008511893956247385 - 0.3535079535658547im -0.008511893956247407 + 0.35350795356585446im -0.006034735471503913 + 0.48183616848092725im -0.006034735471504322 - 0.4818361684809278im; -0.01845528815664635 + 0.3461842706003635im -0.01845528815664656 - 0.346184270600364im -0.029356335987524733 + 0.5180670493895639im -0.029356335987525444 - 0.5180670493895637im; -0.6348652094468953 - 0.006177255420232403im -0.6348652094468951 + 0.006177255420232746im 0.495769296317179 - 0.014144481584118814im 0.4957692963171796 + 0.01414448158411863im; 0.5956309967068604 - 0.009737575514572147im 0.5956309967068606 + 0.00973757551457216im 0.5019127461518936 - 0.03146151541730548im 0.5019127461518929 + 0.03146151541730459im;;; -0.008138987663617082 - 0.3671391151429722im -0.00813898766361732 + 0.3671391151429724im -0.005280592495098128 + 0.477341058334339im -0.005280592495097246 - 0.4773410583343393im; -0.017995252861316983 + 0.336829167879279im -0.017995252861316844 - 0.33682916787927963im -0.030980733899207625 + 0.4974957937801633im -0.030980733899208076 - 0.4974957937801632im; -0.6210876464532612 - 0.004434858231903509im -0.6210876464532608 + 0.0044348582319036125im 0.49166181612743165 - 0.013477185155667528im 0.4916618161274318 + 0.01347718515566628im; 0.6014025897883122 - 0.009048990088989625im 0.6014025897883118 + 0.009048990088989553im 0.5312831603209808 - 0.03344614529637526im 0.5312831603209808 + 0.033446145296376045im;;; … ;;; 0.024828300696341565 - 0.3693314705776267im 0.02482830069634197 + 0.3693314705776269im -0.004668710268486181 + 0.4773911422435022im -0.004668710268486386 - 0.477391142243502im; 0.001664387226165991 + 0.33577804332505595im 0.001664387226165916 - 0.335778043325056im 0.013154588528966692 + 0.4894249368919463im 0.013154588528966692 - 0.4894249368919463im; -0.6193094448617277 + 0.006554960933426945im -0.6193094448617271 - 0.006554960933426337im 0.4880863020199717 + 0.008798957310583203im 0.4880863020199714 - 0.008798957310584386im; 0.60112786069006 + 0.004605437612981434im 0.6011278606900599 - 0.004605437612982686im 0.5417887190232843 + 0.02511499972798819im 0.5417887190232847 - 0.025114999727987815im;;; 0.02807016050049592 - 0.3560094095159763im 0.028070160500495572 + 0.3560094095159765im -0.005323199287426482 + 0.48187882825723954im -0.005323199287425896 - 0.48187882825724im; 0.003082931472690057 + 0.34503902788163665im 0.00308293147269058 - 0.3450390278816365im 0.01910817786687246 + 0.509203679708621im 0.019108177866872206 - 0.5092036797086212im; -0.6326126675056183 + 0.0069633683219010286im -0.6326126675056181 - 0.006963368321900895im 0.49182587839845193 + 0.010660742544823024im 0.49182587839845204 - 0.010660742544823638im; 0.59530365587296 + 0.005101453834684994im 0.5953036558729601 - 0.005101453834685426im 0.5134509272960314 + 0.03282705310110336im 0.5134509272960307 - 0.032827053101103586im;;; 0.01879050679822032 - 0.34503764003874693im 0.01879050679822039 + 0.34503764003874726im -0.00604713541480837 + 0.4852651902080783im -0.006047135414808978 - 0.4852651902080789im; -0.0027358251480628535 + 0.35229504071742235im -0.002735825148062982 - 0.35229504071742285im 0.008070016292990393 + 0.5268124954557836im 0.008070016292989613 - 0.5268124954557827im; -0.643991802924056 + 0.0030235209742055357im -0.6439918029240566 - 0.0030235209742056116im 0.4958228721699905 + 0.004096816026594389im 0.49582287216999 - 0.004096816026594762im; 0.5911493628823044 + 0.0007440251608407484im 0.5911493628823045 - 0.0007440251608402798im 0.48859656431868626 + 0.01804919301137592im 0.488596564318687 - 0.01804919301137582im], [-0.009000302656132372 -0.6862484749213057 0.004276430394466084 0.4864010626399133; -0.01789737059752039 -0.750386082388464 -0.016872245198102154 -0.49890897760399266; 0.7035920189277657 0.009364877919117649 -0.9122427949442852 0.004162582685551378; 0.684000921951247 0.013347013254908632 0.8362210781842051 0.007851883949219417;;; -0.008534404749135334 -0.6814192443076157 -0.01203763587440648 0.49993574233959; -0.041516128495139866 -0.7326574474653328 -0.026099718808632953 -0.4895784905632721; 0.7011236626599215 0.020003317688996946 -0.8978349894786348 0.008735958393535657; 0.7098118127359315 0.044493301795962434 0.8423494337126469 0.013771011357339933;;; -0.007467885523933311 -0.6750621985739491 -0.01151026673787472 0.5192131159128487; -0.04381337405253141 -0.7035652987874754 -0.025449130654807688 -0.47634837741772074; 0.6953148084684015 0.019059618029758457 -0.8783505730365877 0.006271836658759714; 0.7513478507863703 0.04729999228723489 0.8505116989249345 0.012797204509628983;;; … ;;; -0.006602553380484434 -0.675133027917544 0.035112519575444055 0.5223135747020793; 0.018603397505102277 -0.6921513835161869 0.002353798988283846 -0.47486186281739484; 0.6902582679251744 -0.012443604763369713 -0.8758358162292075 -0.009270114652877337; 0.7662049543834749 -0.03551797323431719 0.8501231733082072 -0.006513072332942585;;; -0.007528140627493286 -0.6814795743418439 0.03969720167779063 0.503473335269931; 0.027023044291568047 -0.7201227498542172 0.004359923500545272 -0.4879588727782389; 0.695546827557152 -0.015076566691856348 -0.8946494141154668 -0.009847689920632004; 0.7261292649950893 -0.04642446370832267 0.8418865038658259 -0.007214545200831477;;; -0.00855194091712929 -0.6862686133398245 0.026573789557907134 0.48795691007200254; 0.011412726490118309 -0.7450253759011835 -0.0038690410286716893 -0.49822042453936094; 0.7011994303575825 -0.005793772787357927 -0.9107419417523027 -0.004275904367840944; 0.6909798877883842 -0.025525413546577988 0.83601144637637 -0.0010522104732076053], ComplexF64[-0.04687545134990123 - 5.453620181640978im, -0.04687545134990123 + 5.453620181640978im, -0.04687545134990034 - 6.653620181640956im, -0.04687545134990034 + 6.653620181640956im, -0.046875186380086475 - 6.546404839125937im, -0.046875186380086475 + 6.546404839125937im, -0.04687518638008603 - 7.746404839125932im, -0.04687518638008603 + 7.746404839125932im, -0.04687517363425853 - 4.253595039088093im, -0.04687517363425853 + 4.253595039088093im  …  -0.015624826371528289 - 2.574945434876184im, -0.015624826371528289 + 2.574945434876184im, -0.015624826339016962 - 4.974945436050899im, -0.015624826339016962 + 4.974945436050899im, -0.01562482633901685 - 3.7749454360508903im, -0.01562482633901685 + 3.7749454360508903im, -0.015623791867367176 - 6.175008711448712im, -0.015623791867367176 + 6.175008711448712im, -0.015623791867366732 - 7.375008711448703im, -0.015623791867366732 + 7.375008711448703im])

The frequencies and damping ratios are extracted from the reduced order model using ODEFrequencyDamping

That_ode, Rhat_ode, rho_ode, gamma_ode = ODEFrequencyDamping(MWn, XWn, MSn, XSn, dispMaxAmp)
odeAmp = range(0, dispMaxAmp, length=1000)
odeFreq = abs.(That_ode.(odeAmp))
odeDamp = -Rhat_ode.(odeAmp) ./ odeFreq
1000-element Vector{Float64}:
 0.015242921599632966
 0.0152429223704787
 0.015242924683044052
 0.015242928537413769
 0.015242933933729084
 0.015242940872187756
 0.015242949353043926
 0.015242959376608292
 0.015242970943248021
 0.015242984053386623
 ⋮
 0.04292747374427692
 0.04298695677903298
 0.043046327778356964
 0.043105586767119154
 0.043164733772742096
 0.04322376882516649
 0.04328269195681742
 0.04334150320257098
 0.04340020259972138

The results are then plotted

fig = Figure(size=(1200, 400))
axFreq = Axis(fig[1, 3])
axDamp = Axis(fig[1, 4])
axFreq.xlabel = "Frequency"
axFreq.ylabel = "Amplitude"
axDamp.xlabel = "Damping ratio"
axDamp.ylabel = "Amplitude"

lines!(axFreq, odeFreq, odeAmp, label="ODE O($(ODE_order)) A=$(Amplitude)", linestyle=:dash, linewidth=2)
lines!(axDamp, odeDamp, odeAmp, label="ODE O($(ODE_order)) A=$(Amplitude)", linestyle=:dash, linewidth=2)
CairoMakie.Screen{SVG}

Invariant Manifolds of Maps

First we set up some additional parameters, such as the sampling period $\Delta t$ = Tstep

Tstep = 0.8                 # sampling period
omega = omega_ode * Tstep   # shift angle
0.96

We now create a discrete-time map from the vector field by Taylor expanding an ODE solver about the origin using mapFromODE. We take 500 time steps on the interval $[\theta,\theta + \Delta t]$. The resulting map MP, XPmap is dependent on the phase variable $\theta \in [0,2\pi)$

MP = QPPolynomial(NDIM, NDIM, fourier_order, 0, ODE_order)
XPmap = zero(MP)
mapFromODE(MP, XPmap, shawpierre!, Amplitude, omega_ode, Tstep / 500, Tstep)

The invariant manifold of the map is calculated using QPMAPTorusManifold

MK, XK, MSn, XSn, MWn, XWn, MSd, XSd = QPMAPTorusManifold(MP, XPmap, omega, SEL, threshold = 0.1, resonance = false)
(QPConstant{4, 7, ℝ}(Euclidean(4, 15; field=ℝ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), [-0.5851607959895752 -0.503213852519958 … -0.446152085568514 -0.5648705525668267; -0.6041694556951728 -0.5146118759348656 … -0.47270795577275115 -0.5893510958069279; 0.09125379732965186 0.37071517278053323 … -0.4638298362044751 -0.2056313971558085; 0.11021643017553637 0.39545079914101283 … -0.46470406279621485 -0.19388044333377438], InvariantModels.QPFourierPolynomial{2, 2, 7, 0, 7}([0 0 … 6 7; 0 1 … 1 0], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], Euclidean(2, 36, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; … ;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im], InvariantModels.QPFourierPolynomial{4, 2, 7, 0, 7}([0 0 … 6 7; 0 1 … 1 0], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], Euclidean(4, 36, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im 1.6096564545911838e-10 - 1.8650185095601463e-9im … -3.5585206074910105e-11 + 4.679332879911928e-12im -3.7133476634440594e-12 - 3.299639439844445e-12im; 0.0 + 0.0im 1.257416242045188e-8 + 4.511137062991128e-9im … 2.3821942333749533e-11 - 1.5657058136401012e-11im -9.315050205156803e-12 - 2.8572168253363406e-12im; 0.0 + 0.0im -1.2458189120153121e-8 + 1.273404228583685e-8im … 2.9641285396208623e-11 - 1.5544085189864997e-11im -4.6764449310494455e-12 + 5.5423814120749225e-12im; 0.0 + 0.0im 7.4411904216884545e-9 - 5.81683795980518e-9im … 2.4451483872350513e-11 + 4.160385904592029e-11im -2.6086964975285358e-12 + 1.1082161281661599e-11im;;; 0.0 + 0.0im 6.861334857878946e-9 + 4.199367108945106e-8im … -1.7487397833512394e-5 + 4.207154064487529e-5im 1.9107622796159105e-5 + 3.150705698556681e-5im; 0.0 + 0.0im -2.711534316341898e-9 + 1.7331867195538535e-8im … 8.572739819274295e-6 - 1.9509974981510097e-5im 9.51941003107807e-6 + 1.5600302614286446e-5im; 0.0 + 0.0im 9.817559141260878e-8 - 7.437339069438706e-8im … 8.682486294436047e-5 + 3.2329364973738496e-5im -1.3869108579079906e-6 - 3.878891172519823e-6im; 0.0 + 0.0im 5.427302364824983e-9 + 1.4431214279610586e-9im … -4.0325746350290447e-5 - 1.576829949034412e-5im -6.262677763649995e-7 - 1.9625281123109286e-6im;;; 0.0 + 0.0im 1.9579043532440248e-10 + 3.65644643949092e-10im … -1.0251103897379674e-11 + 3.2533968010698394e-12im -6.479022710317625e-12 + 2.726210353162191e-12im; 0.0 + 0.0im 8.154140792933999e-11 + 2.425983129850051e-10im … -1.0036323777665484e-11 + 1.1872842106814974e-11im -1.094567060799761e-11 - 1.4847824505484044e-12im; 0.0 + 0.0im -3.3901421416019435e-10 + 2.272988738722241e-10im … 2.318023505561409e-11 + 1.942496505715805e-11im -3.373880055651523e-12 - 7.761636969226134e-12im; 0.0 + 0.0im -1.660233299415831e-10 + 3.2752789769369e-11im … -7.852087215131832e-12 + 5.692108433009568e-12im 3.40161507212711e-12 - 1.2901018340165223e-11im;;; … ;;; 0.0 + 0.0im 3.096212306655086e-11 + 1.7424885837807284e-11im … -1.1005097528282161e-11 - 4.856247203335042e-12im -2.298141573370562e-11 + 6.549746409545569e-11im; 0.0 + 0.0im -8.543529677531697e-13 + 7.213359082656068e-12im … -2.087414998052218e-11 + 2.0566956117483358e-11im 5.656845208960112e-11 - 6.54806272033856e-11im; 0.0 + 0.0im 2.2560111187933568e-11 - 4.9471997496152025e-11im … -1.0498893665064553e-10 - 3.619890943147108e-11im 5.852730492489044e-11 + 1.61045575187345e-11im; 0.0 + 0.0im -2.1446339155512725e-13 + 1.1552715345371666e-11im … 4.649521198782412e-11 + 3.004200217694576e-11im -1.172447842756743e-10 - 8.700956405943766e-11im;;; 0.0 + 0.0im -1.1204056271175495e-7 + 5.987111450814777e-8im … 7.617467694197232e-9 - 1.2678059379803262e-8im -6.22581425059412e-10 + 8.212338826975854e-10im; 0.0 + 0.0im 8.479677010524265e-8 - 5.344046468317667e-8im … -5.296013433764207e-10 + 8.914694071940781e-10im 3.8176888017394806e-10 + 2.9034312688446874e-9im; 0.0 + 0.0im -1.0208813106454502e-7 - 1.1433967171430315e-7im … 6.379886570139672e-9 + 3.150798444533642e-8im 3.575844840471151e-10 - 4.02308578333711e-10im; 0.0 + 0.0im -8.721226864950777e-8 - 1.3205579791205493e-7im … 6.150489393470725e-10 - 1.1838681157400594e-9im 3.787159495947737e-9 - 7.864015773799591e-10im;;; 0.0 + 0.0im 2.4574208047818933e-9 + 1.27875416457231e-9im … 1.0937686755065814e-10 - 1.1479050607367259e-10im 4.1697942391586164e-12 - 1.5519490849474802e-11im; 0.0 + 0.0im -9.047632958878595e-10 - 8.624519893287102e-11im … -7.27513595783608e-11 + 8.890014470911479e-11im 2.554238759555487e-12 - 3.668079204654547e-12im; 0.0 + 0.0im 1.7155147163807613e-9 - 2.069240560008937e-9im … -2.7845484757654775e-10 - 2.553722804828445e-10im -7.62177985399976e-12 + 6.105005498328489e-12im; 0.0 + 0.0im -1.926095182900996e-9 + 2.1075921432648457e-9im … 1.6375957172499253e-10 + 1.2558654176089229e-10im 5.604847142951536e-12 + 4.776514499935169e-12im], InvariantModels.QPFourierPolynomial{4, 4, 7, 0, 7}([0 0 … 6 7; 0 0 … 1 0; 0 0 … 0 0; 0 1 … 0 0], [5 4 3 2; 12 9 7 6; … ; 329 326 325 324; 330 329 328 327], [1 1 1 1; 2 2 2 2; … ; 209 209 209 209; 210 210 210 210], [1 1 1 1; 1 1 1 2; … ; 6 2 1 1; 7 1 1 1], Euclidean(4, 330, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im 0.0 + 0.0im … 2.3285055313438408e-12 + 1.7172495042994646e-12im -1.4956848286201933e-13 - 2.6778074180150707e-13im; 0.0 + 0.0im 0.0 + 0.0im … 2.2093741186915324e-12 + 7.536940594329681e-13im -1.79820864592217e-13 - 1.6681144761768238e-13im; 0.0 + 0.0im 0.0 + 0.0im … 3.936003234238607e-12 - 8.41122357713394e-12im -8.515237304892955e-13 + 6.517922202173927e-13im; 0.0 + 0.0im 0.0 + 0.0im … 1.8453235437882974e-12 + 4.621362526745371e-13im -1.527989556237306e-13 - 1.2666644743356996e-13im;;; 0.0 + 0.0im 0.0 + 0.0im … 4.84300751416392e-10 - 1.4661195091379724e-11im -6.637611174242484e-11 - 3.832725705341978e-11im; 0.0 + 0.0im 0.0 + 0.0im … 4.95025603359833e-10 - 1.5114323029071058e-10im -7.733343601856426e-11 - 2.0655866741477422e-11im; 0.0 + 0.0im 0.0 + 0.0im … 3.8999047062266025e-10 + 4.5259234376281066e-11im -4.909648581927023e-11 - 3.8760349459945674e-11im; 0.0 + 0.0im 0.0 + 0.0im … 1.4460596160650998e-11 + 7.710780247258699e-10im 5.821788828949061e-11 - 9.734602721755295e-11im;;; 0.0 + 0.0im 0.0 + 0.0im … 3.323640168765593e-15 + 1.37013869558294e-14im -8.492515535269642e-16 - 1.6095721051744446e-15im; 0.0 + 0.0im 0.0 + 0.0im … 6.199036231395893e-15 + 6.71107049879305e-15im -1.3246509014992684e-15 - 7.66811706015617e-16im; 0.0 + 0.0im 0.0 + 0.0im … 1.0586623897729446e-13 - 4.7568019935386394e-14im -1.3101621539717817e-14 - 3.5723904474437597e-16im; 0.0 + 0.0im 0.0 + 0.0im … 2.228031467095478e-16 + 9.873226293835473e-15im -1.9510566026496376e-16 - 7.005935018050507e-16im;;; … ;;; 0.0 + 0.0im 0.0 + 0.0im … -8.276771099316657e-14 - 8.033970653844721e-14im 5.1601073621943056e-15 + 2.0977412806926803e-14im; 0.0 + 0.0im 0.0 + 0.0im … -1.1460627667137634e-13 - 6.574414028579685e-14im 1.1376239763958226e-14 + 2.121227034025142e-14im; 0.0 + 0.0im 0.0 + 0.0im … -6.600763886498896e-14 - 7.815258291143808e-14im 2.4561769852701826e-15 + 1.8661165965625e-14im; 0.0 + 0.0im 0.0 + 0.0im … 3.0284525726578e-13 - 1.143385061485608e-13im -5.340089587127108e-14 - 8.219034730247659e-15im;;; 0.0 + 0.0im 0.0 + 0.0im … 2.344423142439024e-10 - 4.8903061311881177e-11im -2.6449969523527067e-11 - 7.506371133891471e-12im; 0.0 + 0.0im 0.0 + 0.0im … 1.8013541276113166e-10 - 9.188275242711251e-11im -2.284722924682227e-11 - 3.4139156794216207e-13im; 0.0 + 0.0im 0.0 + 0.0im … -1.7655370972672825e-10 - 4.465609558321994e-10im -8.16611158133065e-12 + 5.8113119700794515e-11im; 0.0 + 0.0im 0.0 + 0.0im … 1.301967527682106e-10 - 9.259779993780517e-11im -1.7853595084766018e-11 + 2.3381669142474515e-12im;;; 0.0 + 0.0im 0.0 + 0.0im … 3.925158583949803e-12 - 7.001652314820555e-12im -1.176663076072344e-12 + 6.860224433939179e-13im; 0.0 + 0.0im 0.0 + 0.0im … 2.094942832197633e-12 - 8.779387803621222e-12im -1.0419684712956307e-12 + 1.0720025440352656e-12im; 0.0 + 0.0im 0.0 + 0.0im … 3.998079916447401e-12 - 5.1921821165808045e-12im -1.034604114688052e-12 + 4.140305326699184e-13im; 0.0 + 0.0im 0.0 + 0.0im … 1.6464622458055132e-11 + 1.1252157194159592e-11im -1.282661660554124e-12 - 2.8189252141344484e-12im], [-0.009000552376276246 -0.6862483185122276 -0.004276407352488767 -0.48640067761565947; -0.017897254651445697 -0.7503865740391369 0.016872272758789334 0.49890932043992986; 0.7035922896796005 0.009364870304002947 0.9122431731000222 -0.004162502641731281; 0.6840002447183959 0.01334678083875924 -0.8362207122536642 -0.007851826177318635;;; -0.008534596625593775 -0.6814190925977878 0.012037558058138342 -0.49993567495520186; -0.04151652881382364 -0.7326579434629639 0.026099504148112254 0.4895784347472987; 0.7011236029152482 0.02000421919688411 0.8978355177434907 -0.008736030312877348; 0.7098116707172549 0.04449317738274341 -0.8423492678160096 -0.013771156786870542;;; -0.007467959735652206 -0.6750622332291538 0.011510588830342125 -0.5192133230764211; -0.04381407966808927 -0.7035652443688287 0.025449446288709952 0.4763480049673945; 0.6953149371364757 0.019059841316049306 0.8783508225111686 -0.006271869904731265; 0.7513477806817049 0.04730036943509936 -0.8505117098738074 -0.012797345094037613;;; … ;;; -0.006602667849535517 -0.6751328653782404 -0.03511271790461955 -0.5223139021277927; 0.018603359885669295 -0.6921508493822353 -0.002354022390841644 0.47486155382407413; 0.6902582512373165 -0.012443437854453555 0.8758357272319904 0.009270269309865293; 0.766205509726251 -0.035518628580773454 -0.8501229850030501 0.006513447369276166;;; -0.007528149818173058 -0.6814793436035967 -0.0396975845313865 -0.5034732632575389; 0.027023444680149444 -0.7201226486455146 -0.0043602414245211924 0.4879587390731162; 0.695546257867568 -0.015076956935466183 0.8946493403431712 0.00984772100919135; 0.7261299611145323 -0.046425055258419036 -0.8418867258560513 0.007214518624960564;;; -0.008551767004856978 -0.6862683172949282 -0.026574080270754093 -0.4879567714208438; 0.011412749522674076 -0.745025583309213 0.0038689823128509595 0.4982203216952609; 0.7011993319284525 -0.0057940196089069015 0.9107423743698615 0.004275731074284006; 0.6909797349226072 -0.025526334684940463 -0.8360110324009611 0.0010519848455631403])

The frequencies and damping ratios are extracted from the reduced order model using MAPFrequencyDamping

That, Rhat_r, rho, gamma = MAPFrequencyDamping(MWn, XWn, MSn, XSn, dispMaxAmp)
mapAmp = range(0, dispMaxAmp, length=1000)
mapFreq = abs.(That.(mapAmp)) ./ Tstep
mapDamp = -log.(abs.(Rhat_r.(mapAmp))) ./ abs.(That.(mapAmp))
1000-element Vector{Float64}:
 0.015242844555505817
 0.01524284531148357
 0.015242847579443537
 0.015242851359465955
 0.01524285665168461
 0.015242863456286546
 0.015242871773512343
 0.01524288160365584
 0.015242892947064546
 0.015242905804139496
 ⋮
 0.04083879342727106
 0.040893794263591465
 0.040948688166818785
 0.041003475142387356
 0.04105815519823288
 0.04111272834475912
 0.04116719459480775
 0.04122155396362614
 0.04127580646883644

The results are plotted

lines!(axFreq, mapFreq, mapAmp, label="MAP O($(ODE_order)) A=$(Amplitude)", linestyle=:dashdot, linewidth=2)
lines!(axDamp, mapDamp, mapAmp, label="MAP O($(ODE_order)) A=$(Amplitude)", linestyle=:dashdot, linewidth=2)
CairoMakie.Screen{SVG}

Invariant Foliations of Maps

Using QPGraphStyleFoliations, two invariant foliations are calculated and the invariant manifold defined by the zero level-set of the second foliation is extracted.

MRf, XRf, MW, XW, MRt, XRt, MUt, XUt, MSt, XSt, MVt, XVt = QPGraphStyleFoliations(MP, XPmap, omega, SEL; dataScale=1, resonance=false, threshold=0.1)
(InvariantModels.QPFourierPolynomial{2, 2, 7, 0, 7}([0 0 … 6 7; 0 1 … 1 0], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], Euclidean(2, 36, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; … ;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im], InvariantModels.QPFourierPolynomial{4, 2, 7, 0, 7}([0 0 … 6 7; 0 1 … 1 0], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], Euclidean(4, 36, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im 1.609655638692909e-10 - 1.8650185960650235e-9im … 1.0290558369103696e-9 - 4.204263442287593e-9im -2.384331883517654e-10 + 7.199827224961264e-10im; 0.0 + 0.0im 1.257416235245072e-8 + 4.511136983425144e-9im … 2.4298816296180224e-9 - 2.535820368995087e-9im -4.4080111368391514e-10 + 1.1736066023472129e-9im; 0.0 + 0.0im -1.245818918630391e-8 + 1.2734042338688093e-8im … 4.666808202934532e-10 - 1.901456464727019e-9im -7.799567025574025e-10 - 1.0508409506706679e-10im; 0.0 + 0.0im 7.441190325469126e-9 - 5.8168378825521616e-9im … -1.002825353280026e-8 + 1.5568268369240198e-9im -7.959713849005924e-10 - 3.9251883785039293e-10im;;; 0.0 + 0.0im 6.86133492241066e-9 + 4.199367109708384e-8im … -1.7901542403742208e-5 + 4.0852223947663156e-5im 1.8652323726665304e-5 + 3.0837363224179226e-5im; 0.0 + 0.0im -2.711534276558906e-9 + 1.733186720640947e-8im … 8.637304783226125e-6 - 1.847992185581663e-5im 9.340793503718753e-6 + 1.5250959899605255e-5im; 0.0 + 0.0im 9.817559143735751e-8 - 7.437339078019805e-8im … 8.295443348263511e-5 + 3.366947495676803e-5im -9.161488077842688e-7 - 4.153052144347647e-6im; 0.0 + 0.0im 5.427302390730187e-9 + 1.4431213359050662e-9im … -3.887791643125982e-5 - 1.6509828632766455e-5im -4.551524719217775e-7 - 1.9850519360720224e-6im;;; 0.0 + 0.0im 1.9579048366536336e-10 + 3.656446717046676e-10im … 1.680323984654669e-11 - 5.49554231756108e-12im -8.399707288127474e-12 + 1.104991702853182e-11im; 0.0 + 0.0im 8.154149073347391e-11 + 2.425983631763377e-10im … 1.0446331113763386e-10 - 2.5993054945393095e-10im -1.762258246352099e-11 + 1.1425589720680598e-11im; 0.0 + 0.0im -3.390141743772027e-10 + 2.2729876261862517e-10im … -4.3720920239570996e-11 + 3.471637783796989e-11im -1.371446144169229e-11 - 1.1626111096268096e-11im; 0.0 + 0.0im -1.660232268990086e-10 + 3.275270048893411e-11im … 3.47203829119864e-10 + 1.7749492081028715e-10im -6.293225061813378e-12 - 2.2471654127537753e-11im;;; … ;;; 0.0 + 0.0im 3.0962061310395114e-11 + 1.7424920532276804e-11im … -1.1835391208963098e-7 - 4.202064839468382e-8im 1.5523579748598905e-7 + 7.153226922257133e-8im; 0.0 + 0.0im -8.544392413340416e-13 + 7.213423614369374e-12im … -2.4607535751155917e-8 - 9.62267102376045e-9im -5.032565542772564e-8 - 1.243915066578639e-8im; 0.0 + 0.0im 2.256013339239406e-11 - 4.947189364403993e-11im … -9.142073719544414e-8 + 1.682652766635009e-7im 3.722681093456883e-7 - 6.217485941075916e-7im; 0.0 + 0.0im -2.1439909113828438e-13 + 1.1552807632660588e-11im … 4.749359690165787e-8 - 1.100801742111024e-7im -2.1853735976882665e-8 + 1.0387899231582906e-7im;;; 0.0 + 0.0im -1.1204056278322555e-7 + 5.987111450884166e-8im … 8.043066013477332e-9 - 1.2384897361438994e-8im -8.262235326359048e-10 + 7.756923777544639e-10im; 0.0 + 0.0im 8.479677006152761e-8 - 5.344046465102646e-8im … -2.717621265604913e-10 + 2.3299495530997415e-9im 4.124195171824434e-10 + 3.0185768041250934e-9im; 0.0 + 0.0im -1.0208813105413667e-7 - 1.1433967164283255e-7im … 8.047075084892903e-9 + 3.041168309584435e-8im 8.076171100854934e-10 - 6.011122659980952e-10im; 0.0 + 0.0im -8.72122686414124e-8 - 1.3205579783086987e-7im … -1.6291328706554152e-9 - 1.3945012781563408e-9im 3.337415826621432e-9 - 6.746798660443002e-10im;;; 0.0 + 0.0im 2.4574208916915393e-9 + 1.2787540850063265e-9im … -6.241104884583714e-6 - 6.8706656831050205e-6im 1.4502089651544031e-5 + 1.0493276339384446e-5im; 0.0 + 0.0im -9.047631900697274e-10 - 8.624524380438494e-11im … 1.4948712051246226e-6 + 5.5970891203389314e-6im -1.0037500127692299e-5 - 8.249953905392769e-6im; 0.0 + 0.0im 1.7155146224743972e-9 - 2.069240641714413e-9im … -1.1653767063828496e-5 + 1.3415115388806707e-5im 1.932658677922084e-5 - 2.6156222493853716e-5im; 0.0 + 0.0im -1.9260952860592186e-9 + 2.1075920663009472e-9im … 9.397992945074158e-6 - 4.067820124174221e-6im -1.3992685206945333e-5 + 1.9007633898723696e-5im], QPPolynomial{2, 2, 7, 0, 7, ℝ}([0 0 … 6 7; 0 1 … 1 0], [0 0 … 6 7; 0 1 … 1 0], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  27, 28, 29, 30, 31, 32, 33, 34, 35, 36], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], [6 5; 9 8; … ; 35 34; 36 35;;; 5 4; 8 7; … ; 34 33; 35 34], [1 1; 2 2; … ; 20 20; 21 21;;; 1 1; 2 2; … ; 20 20; 21 21], [2 1; 2 2; … ; 30 10; 42 6;;; 1 2; 2 6; … ; 10 6; 6 2], Euclidean(2, 36, 15; field=ℝ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), [0.0 0.7220930438571685 … -0.00031548321621645433 7.329805839740776e-5; 0.0 0.6737150271731301 … -0.000238545923187427 3.675185262709774e-5;;; 0.0 0.7220930438571725 … -0.00014203200001346817 1.2704812666770377e-5; 0.0 0.6737150271731338 … 4.2467959784908605e-5 -2.4421560732962452e-5;;; 0.0 0.7220930438571717 … -4.251684149520557e-5 -3.411269266790761e-7; 0.0 0.6737150271731331 … 0.00011163361359871251 -2.6256371316818507e-5;;; … ;;; 0.0 0.7220930438571759 … 4.887690201761155e-5 -2.107365614107228e-5; 0.0 0.6737150271731374 … 1.407238410700124e-6 -2.912204190819436e-5;;; 0.0 0.7220930438571794 … 8.800809751004255e-5 -2.3290535413681465e-5; 0.0 0.6737150271731404 … 0.00011782158898518405 1.0868224063628672e-5;;; 0.0 0.7220930438571695 … -4.743270939876633e-6 4.2133701789233576e-5; 0.0 0.6737150271731314 … -0.00017306343412461505 9.128717536598859e-5], QPPolynomial{2, 4, 7, 0, 7, ℝ}([0 0 … 6 7; 0 0 … 1 0; 0 0 … 0 0; 0 1 … 0 0], [0 0 … 6 7; 0 0 … 1 0; 0 0 … 0 0; 0 1 … 0 0], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  321, 322, 323, 324, 325, 326, 327, 328, 329, 330], [5 4 3 2; 12 9 7 6; … ; 329 326 325 324; 330 329 328 327], [1 1 1 1; 2 2 2 2; … ; 209 209 209 209; 210 210 210 210], [1 1 1 1; 1 1 1 2; … ; 6 2 1 1; 7 1 1 1], [15 14 13 12; 32 29 27 26; … ; 329 326 325 324; 330 329 328 327;;; 14 11 10 9; 29 23 21 20; … ; 326 320 319 318; 329 326 325 324;;; 13 10 8 7; 27 21 18 17; … ; 325 319 317 316; 328 325 323 322;;; 12 9 7 6; 26 20 17 16; … ; 324 318 316 315; 327 324 322 321], [1 1 1 1; 2 2 2 2; … ; 125 125 125 125; 126 126 126 126;;; 1 1 1 1; 2 2 2 2; … ; 125 125 125 125; 126 126 126 126;;; 1 1 1 1; 2 2 2 2; … ; 125 125 125 125; 126 126 126 126;;; 1 1 1 1; 2 2 2 2; … ; 125 125 125 125; 126 126 126 126], [2 1 1 1; 2 1 1 2; … ; 30 10 5 5; 42 6 6 6;;; 1 2 1 1; 1 2 1 2; … ; 10 6 2 2; 6 2 1 1;;; 1 1 2 1; 1 1 2 2; … ; 5 2 2 1; 6 1 2 1;;; 1 1 1 2; 2 2 2 6; … ; 5 2 1 2; 6 1 1 2], Euclidean(2, 330, 15; field=ℝ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), [0.0 1.734723475976812e-18 … 9.665411776457724e-19 0.0; 0.0 3.8293175567871474e-32 … -2.2614227046207647e-19 3.7947076036992655e-19;;; 0.0 6.07153216591883e-18 … 3.4022683551138e-19 2.710505431213761e-20; 0.0 -5.551115123125777e-17 … 9.381589043055761e-19 1.0842021724855044e-18;;; 0.0 -1.7347234759768042e-18 … 6.293283677224624e-19 4.336808689942018e-19; 0.0 -1.110223024625156e-16 … -1.9979665549775845e-19 7.589415207398531e-19;;; … ;;; 0.0 -1.734723475976804e-18 … 1.9476139309456128e-18 -1.4907779871675686e-19; 0.0 -5.551115123125775e-17 … 5.559435726844167e-19 -1.5720931501039814e-18;;; 0.0 -1.734723475976798e-18 … -6.026800587308207e-19 -1.6263032587282567e-19; 0.0 -5.4484646120902355e-33 … -7.725057425165442e-19 8.131516293641283e-19;;; 0.0 -8.67361737988408e-19 … 3.903308430219113e-18 2.710505431213761e-19; 0.0 -1.110223024625159e-16 … 2.058340115959879e-18 -1.0570971181733668e-18], QPPolynomial{2, 2, 7, 0, 7, ℝ}([0 0 … 6 7; 0 1 … 1 0], [0 0 … 6 7; 0 1 … 1 0], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  27, 28, 29, 30, 31, 32, 33, 34, 35, 36], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], [6 5; 9 8; … ; 35 34; 36 35;;; 5 4; 8 7; … ; 34 33; 35 34], [1 1; 2 2; … ; 20 20; 21 21;;; 1 1; 2 2; … ; 20 20; 21 21], [2 1; 2 2; … ; 30 10; 42 6;;; 1 2; 2 6; … ; 10 6; 6 2], Euclidean(2, 36, 15; field=ℝ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), [0.0 0.9487048512517642 … 0.00025875978576902613 -0.00015179905445110523; 0.0 0.166440513030325 … -0.00012476403232781726 -2.2742936356266073e-5;;; 0.0 0.948704851251769 … 0.00024149645890246578 -0.00011386088888304673; 0.0 0.16644051303032575 … -6.58657198744062e-5 -2.472527198045435e-5;;; 0.0 0.9487048512517685 … 0.00017960435260394315 -8.006440494626287e-5; 0.0 0.1664405130303256 … -4.091627675545365e-5 -1.3054815484393508e-5;;; … ;;; 0.0 0.9487048512517734 … 0.00023686436220334768 -0.00014617890417471807; 0.0 0.1664405130303266 … -0.0001460951720448794 -5.285485922019452e-6;;; 0.0 0.9487048512517781 … 0.00027208624827358534 -0.00016942302736344702; 0.0 0.1664405130303273 … -0.00016311149617382584 -8.531777736003058e-6;;; 0.0 0.9487048512517651 … 0.0002681253270619172 -0.00017072064539183006; 0.0 0.1664405130303251 … -0.00015816006476494034 -1.364306540319621e-5], QPPolynomial{2, 4, 7, 0, 7, ℝ}([0 0 … 6 7; 0 0 … 1 0; 0 0 … 0 0; 0 1 … 0 0], [0 0 … 6 7; 0 0 … 1 0; 0 0 … 0 0; 0 1 … 0 0], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  321, 322, 323, 324, 325, 326, 327, 328, 329, 330], [5 4 3 2; 12 9 7 6; … ; 329 326 325 324; 330 329 328 327], [1 1 1 1; 2 2 2 2; … ; 209 209 209 209; 210 210 210 210], [1 1 1 1; 1 1 1 2; … ; 6 2 1 1; 7 1 1 1], [15 14 13 12; 32 29 27 26; … ; 329 326 325 324; 330 329 328 327;;; 14 11 10 9; 29 23 21 20; … ; 326 320 319 318; 329 326 325 324;;; 13 10 8 7; 27 21 18 17; … ; 325 319 317 316; 328 325 323 322;;; 12 9 7 6; 26 20 17 16; … ; 324 318 316 315; 327 324 322 321], [1 1 1 1; 2 2 2 2; … ; 125 125 125 125; 126 126 126 126;;; 1 1 1 1; 2 2 2 2; … ; 125 125 125 125; 126 126 126 126;;; 1 1 1 1; 2 2 2 2; … ; 125 125 125 125; 126 126 126 126;;; 1 1 1 1; 2 2 2 2; … ; 125 125 125 125; 126 126 126 126], [2 1 1 1; 2 1 1 2; … ; 30 10 5 5; 42 6 6 6;;; 1 2 1 1; 1 2 1 2; … ; 10 6 2 2; 6 2 1 1;;; 1 1 2 1; 1 1 2 2; … ; 5 2 2 1; 6 1 2 1;;; 1 1 1 2; 2 2 2 6; … ; 5 2 1 2; 6 1 1 2], Euclidean(2, 330, 15; field=ℝ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), [0.0 -4.9873299934333435e-18 … -0.025813233322878437 -0.013516641917457424; 0.0 1.0000000000000009 … 0.01750128328909125 -0.007810768725316193;;; 0.0 1.778091562876229e-17 … 0.02531200794822723 -0.013181471829324219; 0.0 1.0000000000000004 … 0.018483064200764298 0.010967729261208991;;; 0.0 -8.23993651088983e-18 … 0.02169790253884563 0.0060421473245413975; 0.0 1.0 … -0.0324489229907649 0.010595338071659794;;; … ;;; 0.0 -7.806255641895627e-18 … 0.010211488184901122 -0.009580633579270243; 0.0 0.9999999999999999 … 0.023209420493846427 0.002028130929866181;;; 0.0 1.1275702593849262e-17 … 0.02424038341030582 0.0025181837323918063; 0.0 1.0 … -0.019040349529343005 0.006438780488369999;;; 0.0 -8.67361737988405e-19 … -0.020813467437709978 0.003702278439191366; 0.0 1.0000000000000002 … -0.02988799680808855 -0.007831186789279541])

The results are plotted

That, Rhat_r, rho, gamma = MAPFrequencyDamping(MW, XW, MRf, XRf, dispMaxAmp)
foilAmp = range(0, dispMaxAmp, length=1000)
foilFreq = abs.(That.(foilAmp)) ./ Tstep
foilDamp = -log.(abs.(Rhat_r.(foilAmp))) ./ abs.(That.(foilAmp))

lines!(axFreq, foilFreq, foilAmp, label="FOIL O($(ODE_order)) A=$(Amplitude)", linestyle=:dashdot, linewidth=2)
lines!(axDamp, foilDamp, foilAmp, label="FOIL O($(ODE_order)) A=$(Amplitude)", linestyle=:dashdot, linewidth=2)
CairoMakie.Screen{SVG}

ROM identification from data

Set-up

Importing a library to load/store data and specifying parameters

using BSON: @load, @save
# parameters for the data-driven part
maxSimAmp = 1.0     # maximum initial condition measured from the torus
maxAmp = 1.0        # filter out data with amplitude greater than maxAmp
dataRatio = 1.0     # the proportion of data to be retained when filtering (all of it)
R_order = 7         # polynomial order of R
U_order = 7         # polynomial order of U
V_order = 7         # polynomial order of V
S_order = 7         # polynomial order of S
STEPS = 800         # number of optimisation steps to use

# the names of data files
revision = "shawpierre-2p$(log2(scale_epsilon))-SIMAMP$(maxSimAmp)-AMP$(Amplitude)-F$(fourier_order)-R$(R_order)-U$(U_order)-V$(V_order)-S$(S_order)-MODE$(SEL[1])"
datarevision = "shawpierre-SIMAMP$(maxSimAmp)-AMP$(Amplitude)-F$(fourier_order)"
"shawpierre-SIMAMP1.0-AMP0.25-F7"

Creating data

We create 600 trajectories 60 points long each using generate and save the data for future use.

dataX, dataY, thetaX, thetaY, thetaNIX, thetaNIY = generate(NDIM, shawpierre!, ones(NDIM) * maxSimAmp, 600, 50, fourier_order, omega_ode, Tstep, Amplitude, XWre)
size(b) = (4, 100)
Periodicity of solution = 1.2145799233068746e-10 norm of solution = 1.0164815066919473

Finding an Approximate Linear Model

An approximate linear model is found about the invariant torus using findLinearModel.

A, b, MK1, XK1 = findLinearModel(dataX, thetaX, dataY, thetaY, omega)
([0.4086242429013703 0.24570929548960743 0.6172777849477425 0.08002782001583861; 0.24471607626341502 0.4512346778802403 0.08038228603204274 0.6275629403839792; -1.2399926016923382 0.4584560734027012 0.3876967323002161 0.26002105167589357; 0.45237546558995384 -1.1748209705653996 0.25985941439496774 0.4144041444258385;;; 0.42632853647074104 0.24729217107041498 0.6228267590663608 0.08114351418520128; 0.24584419953486936 0.4512958403207193 0.08060461463723662 0.6276261415087447; -1.2067643350490482 0.46397404682914156 0.40471216200329424 0.26276512723689605; 0.4575496781451996 -1.1744670696327812 0.2611435900221003 0.414707021512516;;; 0.4429559752827833 0.24772046005774376 0.6260087328346333 0.08118566015932674; 0.24714415437255685 0.4513023592362667 0.08075844610511702 0.6276352204103651; -1.1855022007739902 0.46508422951565037 0.41047001200908656 0.2623767898797001; 0.46299617399044773 -1.1743897547480555 0.26197646972012745 0.41473598434014275;;; … ;;; 0.42587409364445966 0.24568989962722695 0.6178038414023674 0.08031381295973287; 0.24653574637064796 0.4512654978458718 0.08055398043217464 0.6275986400380157; -1.2379391285586372 0.4540142045613102 0.3740005748944737 0.2587624321735838; 0.4591834347358553 -1.1747201598460826 0.2604500044910681 0.41453932287290424;;; 0.4050036311193514 0.24556549191388033 0.6135794356089468 0.0797660728239932; 0.2450236909626848 0.45127360954216506 0.08035690786810828 0.6275572980090773; -1.2694682043105436 0.45338091723637475 0.3645193827686903 0.25807904193958026; 0.452638210931208 -1.1747174260103674 0.25937898910153384 0.4143635945785931;;; 0.3988174604628312 0.24497354650427483 0.6125276204016864 0.07999456893182505; 0.24431809598597204 0.4512345781443257 0.0802230088646962 0.6275637536226277; -1.2678519821894076 0.45342219840893844 0.3689923217739597 0.2593877455762531; 0.4500648581210388 -1.1749032552261858 0.25886751948219944 0.41440569860250753], [0.05021660556126542 0.041513559301964154 … 0.047620845128613296 0.0500414649041144; 0.06807830947070428 0.052945255388523134 … 0.06311365942241876 0.07157296888884442; 0.11389670839387196 0.08032099438708029 … 0.1154990407381776 0.11995234844404365; 0.14309300277485607 0.09833190859299178 … 0.15742070670577504 0.16342678530522262], QPConstant{4, 7, ℝ}(Euclidean(4, 15; field=ℝ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), [-0.5874229919596033 -0.5049821138811313 … -0.44821036706860035 -0.5672378651305572; -0.608331751231711 -0.5180007802023747 … -0.4760712976330928 -0.5933937161133606; 0.0920659444167287 0.3726752727491479 … -0.4652862441715301 -0.2061190780076265; 0.11122799474332956 0.3982829900397215 … -0.4677658763695057 -0.1949298693157646])

Identifying Invariant Vector Bundles and Transforming Data

We calculate invariant vector bundles from the linear model and project the data into this coordinate system using QPPreProcess

thetaTX, dataTX, thetaTY, dataTY, dataScale, preId, R1, S1, W1 = QPPreProcess(XK1, A, omega, thetaX, dataX, thetaY, dataY, SEL; Tstep=Tstep, maxAmp=maxAmp, data_ratio=dataRatio)
([0.08214481661386881 0.04822944537631952 … -0.060466527784439686 -0.06380744939668545; -0.09968236131744047 -0.05103271544982402 … 0.06910110892233844 0.06496174888271398; … ; 0.06761028781160111 0.04930142307272221 … -0.054218105069802104 -0.07046210573455357; -0.0725513844907744 -0.04769403008009382 … 0.055922750868088865 0.0655064164962598], [0.0051215384718913065 0.0019795561426470686 … 0.1825500069192415 -0.0018701665342527685; -0.0010492323634180407 -0.004049664674296287 … -0.1650787112151768 -0.24007528585448598; -0.0009976108200256595 0.003782924596626657 … 0.08849407628653828 -0.06954682639976897; 0.002060139516033514 0.0006752129563365986 … -0.048190627686765956 -0.055716342611507735], [0.04822944537631952 -0.008481422544779649 … -0.06380744939668545 -0.026713868919854878; -0.05103271544982402 0.008096353120604827 … 0.06496174888271398 0.02447148416250472; … ; 0.04930142307272221 -0.010905038381890173 … -0.07046210573455357 -0.038629793764613375; -0.04769403008009382 0.009333225399794778 … 0.0655064164962598 0.030895134887977793], [0.0019795561426470686 -0.002435523643813768 … -0.0018701665342527685 -0.177684142302258; -0.004049664674296287 -0.0038359919902718363 … -0.24007528585448598 -0.16324102405633623; 0.003782924596626657 0.002232357365092972 … -0.06954682639976897 -0.10541586075415102; 0.0006752129563365986 -0.0021957904287901106 … -0.055716342611507735 0.028631207487507228], 1.52758929291458, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  30591, 30592, 30593, 30594, 30595, 30596, 30597, 30598, 30599, 30600], [0.6653633281317204 0.748230282656151; -0.6984293799211735 0.6803914149169781;;; 0.6651824404591532 0.7351161059944153; -0.7117980311008933 0.6792550524517846;;; 0.6630662083360375 0.7273870201949437; -0.7187502170166337 0.6818751835657852;;; … ;;; 0.6586968544730138 0.7457067302865419; -0.6999279420141792 0.6873349890007365;;; 0.6617222107106276 0.7559220690406279; -0.688937776195794 0.6853396752436565;;; 0.6637403681916549 0.7565312126125183; -0.6901282610964792 0.6834890697816636], [0.1480763554062821 1.663486379100509; -0.5418509000146974 0.181123690058478;;; 0.14995447638635775 1.6511246596144866; -0.5446752583937363 0.18305467685273655;;; 0.1511146324645231 1.644221579020313; -0.5469818872804821 0.1846202993461625;;; … ;;; 0.14920438226942445 1.662897310047963; -0.5412451943409523 0.1840650375878355;;; 0.14870735137132368 1.674516472891695; -0.5385941526161989 0.18157097474301048;;; 0.14699595511963104 1.673773826576939; -0.538718856029738 0.17964862998309553], [-0.002764489327837947 0.6746185201630671 -0.0007195816796218986 -0.6956347621731589; 0.003695694899653377 0.7412893183050949 0.007238238201222368 0.7214318596508607; -0.7161134870874257 0.002640053726505548 0.7374544603399755 -0.013148387876245398; -0.7000292061393409 -0.0027282252316265716 -0.677567234728704 -0.005067817850113792;;; -0.021666242193223596 0.6806639659741923 0.011602584070616304 -0.7231676491696283; 0.011290870896278504 0.7332476246264413 0.016751333452934034 0.6919447422647891; -0.7027222695105149 0.0067882723597026805 0.734536705502992 -0.018337826391390116; -0.7133493267928397 -0.014890803661048434 -0.6798602956939654 -0.014044200327484517;;; -0.02153489413544878 0.6928922526977811 0.01562523347329142 -0.7359292615657115; 0.012313299971328302 0.7226906101313578 0.021965546021035388 0.6783802595647508; -0.6796785966603791 0.010152253355659627 0.7176859881104837 -0.012543578123898176; -0.7344589433426322 -0.018422161713172312 -0.6979242868179965 -0.007781110237348229;;; … ;;; 0.021220666122702805 0.699198705320369 -0.02108979520265602 -0.7373686226587729; -0.004227037777373873 0.7160008605949262 -0.005018608846844598 0.6770648463249777; -0.6698201881282898 -0.0029528660146011755 0.7160266564108254 0.006922139150258803; -0.74489604404665 0.018465635047340716 -0.699840271083632 0.025415490641732957;;; 0.026743504664513956 0.6861796651013425 -0.02811227466882429 -0.722242358206079; -0.0056695495722592385 0.7288093962659346 -0.011284926566145761 0.6921093473340896; -0.6892284161084979 -0.005058322637359919 0.7293800379539379 -0.005062643646452488; -0.7243064950295437 0.02436422567602862 -0.6850963538280179 0.013316244752296508;;; 0.01646982701426667 0.6775140162473786 -0.012950137049536013 -0.7037517050811792; -0.002821569056781108 0.7370311882805061 0.0012378477308987853 0.7123947726187274; -0.7146892248295533 -0.003885678280140689 0.7405669550387988 0.00025813027638928843; -0.702195790521915 0.016945865984045398 -0.6745243094254985 0.014503598229909164])

Here, the inverse transformation W1 was also created.

Setting up the Function Approximators

We create the data structure holding our two invariant foliations using QPCombinedFoliation and set up the data cache that speeds up calculations using makeCache. Our optimisation algorithm requires setting up so-called trust-region radii. For housekeeping purposes we also specify using dataIdV that all data points are used for identifiction.

MCF, XCF = QPCombinedFoliation(NDIM, 2, fourier_order, R_order, U_order, S_order, V_order, R1, S1, MK1, zero(XK1), sparse=false)
# creating a cache
XCFcache = makeCache(MCF, XCF, thetaTX, dataTX, thetaTY, dataTY)
dataIdV = Array{Any,1}(undef, 1)
dataIdV[1] = 1:size(thetaTX, 2)
radii = to_zero(XCF)
([0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0], [0.0])

Performing the Optimisation

Using function QPOptimise, the optimal parameter values of our foliation are found:

QPOptimise(MCF, XCF, thetaTX, dataTX, thetaTY, dataTY;
        maxit=8, gradient_ratio=2^(-7), gradient_stop=2^(-29), steps=STEPS, name=revision,
        cache=XCFcache,
        omega=omega, Tstep=Tstep, dataId=dataIdV[1], radii=radii)

Analysis of the result

The post processing step using QPPostProcess extract a reduced order model together with an invariant manifold. The manifold immersion $\boldsymbol{W}$ is MW, XW.

MSn, XSn, MFW, XFWoWdoWn, MW, XW = QPPostProcess(MCF, XCF, W1, omega)
(InvariantModels.QPFourierPolynomial{2, 2, 7, 1, 7}([0 1 … 6 7; 1 0 … 1 0], [4 3; 5 4; … ; 34 33; 35 34], [1 1; 2 2; … ; 26 26; 27 27], [1 2; 2 1; … ; 6 2; 7 1], Euclidean(2, 35, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; … ;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im;;; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im … 0.0 + 0.0im 0.0 + 0.0im], InvariantModels.QPFourierPolynomial{4, 2, 7, 0, 7}([0 0 … 6 7; 0 1 … 1 0], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], Euclidean(4, 36, 15; field=ℂ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), ComplexF64[3.748154945535523e-6 - 3.4761333331970852e-6im 5.650773112172178e-5 - 0.00010773518132552737im … 0.004165756839416029 - 3.5971947736521516e-5im 0.002036822452356148 - 0.0008596109444622484im; 2.5727289900518823e-6 - 2.1079357145511634e-6im 8.125072691527097e-5 - 0.00018397337006952684im … -0.006080998073588077 + 0.002937996755224843im 0.004281723103034676 + 0.0021055981282649574im; -1.967244535651938e-6 - 2.3910614996249956e-6im 0.00018550294551551456 + 9.47988055155538e-5im … 0.003364114239012189 - 1.2846113663589342e-5im -0.0025821872527712936 - 0.0017565352574889699im; -1.254370626746792e-6 - 1.285376608318247e-6im 9.471520134132838e-5 + 9.757991328829008e-5im … 0.0007273573184245898 + 0.0024261408191444973im 0.0037429285336958475 - 0.005891709316150719im;;; 8.549219314622099e-6 - 6.415832757193304e-7im 9.500171506521035e-6 + 8.585403255072211e-6im … -0.008679379770562821 + 0.0064892502265998im 0.001480108238682105 - 0.00040380162051474697im; 5.423922615092493e-6 - 6.234461287562797e-7im 1.0620807200681042e-5 + 8.262787527217473e-6im … 0.008653532806263128 - 0.005895523417716988im -0.0017891373745823544 + 0.0008973009437978202im; 6.901926192737942e-7 + 5.011031306903287e-6im 6.393433731230431e-6 - 1.4547104308052248e-6im … 0.013967646560412628 + 0.013647181689687336im -0.003032256665917407 + 7.031868575712707e-5im; 5.905681403773868e-7 + 3.560244781991945e-6im 5.9135832994424715e-6 + 1.1084338224480613e-6im … -0.016021978544987783 - 0.01060373871991645im 0.0018524924532686356 + 0.0002869703488352643im;;; -1.5911959515382435e-5 - 3.03576319374013e-6im 1.4146823048675298e-8 - 7.542502988710771e-6im … -0.03276038830303475 - 0.022898163819335588im -0.006538870395308595 - 0.0052969393379729175im; 1.066235941721551e-5 + 3.5092837825593427e-6im -4.1424798863382904e-7 - 4.783577167113424e-6im … -0.03389676075173591 - 0.022968163683066523im -0.0024116506528880067 + 0.009874176338548217im; 4.242295955573176e-6 - 2.9563399475654668e-5im 9.133864021164256e-6 + 1.3841172792086036e-6im … -0.01884005863399038 + 0.03457158307846331im 0.013004062120916768 - 0.0025942861172511085im; -6.535831777701701e-6 + 1.977910621576006e-5im 3.3618993252691184e-6 - 1.5424111144394281e-6im … -0.022830919674341664 + 0.03315329191146816im -0.017833742021896547 - 0.005682854751637289im;;; … ;;; -1.5911959515382435e-5 + 3.03576319374013e-6im 1.697703492675864e-6 + 4.310372642111495e-8im … -0.007755960342991895 - 0.00034441576809611063im 0.003355609748154119 + 0.0012806139150630309im; 1.066235941721551e-5 - 3.5092837825593427e-6im 3.503578181783641e-7 + 2.99918510557406e-7im … 0.005643727643878025 - 0.0017517625058193983im -0.0021818627368493195 - 0.0004802595351819194im; 4.242295955573176e-6 + 2.9563399475654668e-5im 1.3541942770852174e-6 - 2.7798577157580317e-6im … -0.00281430549098232 - 0.009155234273621644im 0.004982055760425196 - 0.004749547583478261im; -6.535831777701701e-6 - 1.977910621576006e-5im 3.004218218473836e-8 + 4.128373983887246e-7im … 0.0021498339733805483 + 0.0055162533431175245im -0.004227223710277882 + 0.0032706317708465016im;;; 8.549219314622099e-6 + 6.415832757193304e-7im 2.6037283490381162e-5 - 2.100307006640644e-5im … -0.0010040169121699955 + 0.0001220095744865438im -0.0033263615022458785 - 0.0016914548795633568im; 5.423922615092493e-6 + 6.234461287562797e-7im -3.812400167647038e-5 + 2.5872649171284368e-6im … 0.0017371503576953937 + 0.0027924421236045433im 0.0013813082617924269 - 0.0028255149069001465im; 6.901926192737942e-7 - 5.011031306903287e-6im -3.5204520425944564e-5 - 4.823204357330762e-5im … -0.0005317479970875858 + 0.0015802638972775155im -0.001115547217309372 + 0.005143980735737015im; 5.905681403773868e-7 - 3.560244781991945e-6im 1.1693192271568192e-5 + 5.554689555692434e-5im … 0.002703979008379667 - 0.0014724593047048334im -0.0031732502614707173 - 0.003792992746868568im;;; 3.748154945535523e-6 + 3.4761333331970852e-6im 7.054174191524772e-5 + 2.7047943712876494e-5im … 0.005766399736424097 + 0.016590587038670602im 0.00012020885656378717 - 0.0012457894871738251im; 2.5727289900518823e-6 + 2.1079357145511634e-6im 8.048719553406593e-5 + 2.6940039398187684e-5im … -0.004665554799155152 - 0.013117392608790116im 0.0004910358437925037 + 0.0019182403218646745im; -1.967244535651938e-6 + 2.3910614996249956e-6im 2.651042847048994e-5 - 6.819345290096938e-5im … 0.02585514108788617 - 0.0165875224047351im 0.00030354237507421427 + 0.002687824360939107im; -1.254370626746792e-6 + 1.285376608318247e-6im 2.864087343841629e-5 - 8.799610968050241e-5im … -0.01924419025923426 + 0.015491127497662295im -0.000547247673639499 - 0.001983148839963727im], QPPolynomial{4, 2, 7, 0, 7, ℝ}([0 0 … 6 7; 0 1 … 1 0], [0 0 … 6 7; 0 1 … 1 0], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  27, 28, 29, 30, 31, 32, 33, 34, 35, 36], [3 2; 5 4; … ; 35 34; 36 35], [1 1; 2 2; … ; 27 27; 28 28], [1 1; 1 2; … ; 6 2; 7 1], [6 5; 9 8; … ; 35 34; 36 35;;; 5 4; 8 7; … ; 34 33; 35 34], [1 1; 2 2; … ; 20 20; 21 21;;; 1 1; 2 2; … ; 20 20; 21 21], [2 1; 2 2; … ; 30 10; 42 6;;; 1 2; 2 6; … ; 10 6; 6 2], Euclidean(4, 36, 15; field=ℝ), ManifoldsBase.ExponentialRetraction(), ParallelTransport()), [0.0014809240361318996 0.6763781671143722 … 0.039637206599849144 0.017394876770976873; 0.0027257635592839924 0.739506111573531 … -0.02408551437167548 -0.015116198445182343; -0.0005311978036175632 0.0031536541534532593 … 0.025653150075140317 0.01817177214542095; -0.0006616417430165843 -0.0012813074573411207 … -0.025788547550563243 -0.019032233607934265;;; 0.0011570478220440427 0.6820750508739165 … -0.0337972383435625 0.014934031325871656; 0.002217196295516878 0.7320124402801984 … 0.024825986244242666 -0.007582286394094781; -0.0012836637861892656 0.009227962361463216 … -0.08829216280660256 -0.0296675572307942; -0.0018543912022699063 -0.014915613100258492 … 0.09172916730256912 0.023494383890657458;;; 0.0006934821636589871 0.6939409221558587 … -0.006075138063008068 0.005134897574939863; 0.001392836923230783 0.7216694109212408 … 0.006706014802676321 -0.011526412906674244; -0.001633363698487912 0.0101862464837176 … 0.00366742288270383 -0.031714554123578015; -0.002865282448176778 -0.016318771071899035 … 0.009102547087697229 0.03312451879875488;;; … ;;; 0.0008491784314074188 0.6991413835235177 … -0.009121304194821671 -0.015579319198881727; 0.0013712335380870737 0.716012792693284 … 0.0033118763399096586 0.02244201578790692; 0.0015390142897013614 -0.00217862287874621 … -0.022461763714390248 -0.024031734887851236; 0.002850227171376962 0.020145177783048817 … 0.022342786840250744 0.029185390109352404;;; 0.0013468256490110164 0.6884139903832518 … 0.04149848963758048 -0.022556495902726438; 0.002201601407633179 0.7266441330893807 … -0.04523404930491276 0.024687039251808342; 0.0009527543304164576 -0.004759576171177643 … -0.0016428382856081417 -0.01168808405252398; 0.002004382008091476 0.025832102116310207 … -0.0050707960345592025 0.012668650659340312;;; 0.001549932584757837 0.6783908977456736 … 0.029591997994725786 0.006577932019771862; 0.0026462224202564015 0.7361312116472596 … -0.027872311940278633 -0.0010556111849949839; 0.0003195160454726562 -0.002070924762742828 … 0.007215339547157344 0.012798083447965578; 0.0006866879766167695 0.01772433033202346 … -0.00508406468309915 -0.013650340015773614])

The instantaneous frequencies and damping ratios are also plotted

That, Rhat_r, rho, gamma = MAPFrequencyDamping(MFW, XFWoWdoWn, MSn, XSn, dispMaxAmp / dataScale)
r_data = range(0, dispMaxAmp / dataScale, length=1000)
dataFreq = abs.(That.(r_data)) / Tstep
dataDamp = -log.(abs.(Rhat_r.(r_data))) ./ abs.(That.(r_data))
dataAmp = r_data .* dataScale
# plotting
lines!(axFreq, dataFreq, dataAmp, label="DATA O($(R_order), $(S_order))", linestyle=:solid, linewidth=2)
lines!(axDamp, dataDamp, dataAmp, label="DATA O($(R_order), $(S_order))", linestyle=:solid, linewidth=2)
CairoMakie.Screen{SVG}

Error Analysis

To make sense of the results, we analyse the fitting error, which is the residual of the invariance equation divided by the distance from the invariant torus

\[ E = \frac{\left| \boldsymbol{R}\left(\boldsymbol{U}\left(\boldsymbol{x}_k,\theta_k\right),\theta_k\right) - \boldsymbol{U}\left(\boldsymbol{y}_k,\theta_k +\omega\right) \right|}{\left| \boldsymbol{x}_k - K(\theta_k) \right|}\]

The amplitude of a given point is calculated by

\[ A = \left| \boldsymbol{W}(\boldsymbol{U}(\boldsymbol{x}_k, \theta_k)) \right|\]

The following call calculates $A$ as OnManifoldAmplitude and various statistics about E as a function of the amplitude.

OnManifoldAmplitude, hsU, errMaxX, errMaxY, errMinX, errMinY, errMeanX, errMeanY, errStdX = ErrorStatistics(MCF, XCF, MW, XW, thetaTX, dataTX, thetaTY, dataTY; dataScale=dataScale, cache=XCFcache)
den = Makie.KernelDensity.kde(sort(vcat(OnManifoldAmplitude, -OnManifoldAmplitude)))
atol = eps(maximum(den.density))
den.density[findall(isapprox.(den.density, 0, atol=atol))] .= atol
dataDensityX = den.density
dataDensityY = den.x
-1.1531227997789346:0.0011266466045715042:1.1531227997789346

The error and data distribution is the plotted

axErr = Axis(fig[1, 2], xscale=log10)
axDense = Axis(fig[1, 1], xscale=log10)
xlims!(axErr, 1e-6, 1e-1)
ylims!(axErr, 0, maxAmp)
xlims!(axDense, 1e-2, 100)
ylims!(axDense, 0, maxAmp)

lines!(axDense, dataDensityX, dataDensityY)
heatmap!(axErr, hsU, colormap=CairoMakie.Reverse(:greys))
lines!(axErr, errMaxX, errMaxY, linestyle=:solid, linewidth=2, color=:red)
lines!(axErr, errMinX, errMinY, linestyle=:solid, linewidth=2, color=:red)
lines!(axErr, errMeanX, errMeanY, linestyle=:dash, linewidth=2, color=:blue)
lines!(axErr, errMeanX .+ errStdX, errMeanY, linestyle=:dot, linewidth=2, color=:green)
MakieCore.Lines{Tuple{Vector{GeometryBasics.Point{2, Float32}}}}

Finally plot legends are displayed

# creating the legend
fig[1, 5] = Legend(fig, axFreq, merge=true, unique=true, labelsize=16, backgroundcolor=(:white, 0), framevisible=false, rowgap=1)
resize_to_layout!(fig)
CairoMakie.Screen{SVG}

Functional Representation

All data structures follow the conventions of ManifoldsBase.jl. For example, the representation of the submersion $\boldsymbol{U}$ is given as two components MU and XU, where MU describes the function and XU contains the parameters of $\boldsymbol{U}$.

Interpolation in Fourier Space

When representing periodic functions, we use Fourier collocation. The grid of the collocation is given by

\[\vartheta_{1}=0,\ldots,\vartheta_{2\ell+1}=\frac{2\ell}{2\ell+1}2\pi,\]

where $\ell$ is the highest resolved Fourier harmonic. The grid has $2\ell+1$ nodes, which is the same as the number of Fourier coefficients that can uniquely represent a periodic function on the grid.

Given a set of points on the grid $\boldsymbol{x}_{j}=\boldsymbol{x}\left(\vartheta_{j}\right)$, we can use an interpolation to reconstruct the periodic function on the interval $[0,2\pi)$ using

\[\boldsymbol{x}\left(\theta\right)=\sum_{j=1}^{2\ell+1}\gamma\left(\theta-\vartheta_{j}\right)\boldsymbol{x}_{j},\]

where

\[\gamma\left(\theta\right)=\frac{1}{2\ell+1}\frac{\sin\frac{2\ell+1}{2}\theta}{\sin\frac{1}{2}\theta}.\]

Let us define the matrix $\boldsymbol{X} = (\boldsymbol{x}_{j}, \ldots, \boldsymbol{x}_{2\ell+1})$ and the interpolation vector for a given value of $\theta$ by

\[\boldsymbol{t}_\theta = (\gamma\left(\theta-\vartheta_{1}\right), \ldots, \gamma\left(\theta-\vartheta_{2\ell+1}\right))^T,\]

then the function evaluation $\boldsymbol{x}\left(\theta\right)$ is a matrix-vector product

\[\boldsymbol{x}\left(\theta\right)= \boldsymbol{X} \boldsymbol{t}_\theta.\]

Input to Data-driven Methods

We assume that the data is produced by an unknown map $\boldsymbol{F}$, such that

\[\boldsymbol{y}_{k} = \boldsymbol{F}\left(\boldsymbol{x}_{k},\theta_{k}\right).\]

When considering trajectories for some set of consecutive indices we have $\boldsymbol{x}_{k+1}=\boldsymbol{y}_{k}$. The data we analyse is organised into matrices

\[\begin{aligned} \boldsymbol{X} &= (\boldsymbol{x}_1,\ldots, \boldsymbol{x}_N) \\ \boldsymbol{Y} &= (\boldsymbol{y}_1,\ldots, \boldsymbol{x}_N) \\ \boldsymbol{\Theta}_x &= (\boldsymbol{t}_{\theta_1},\ldots, \boldsymbol{t}_{\theta_N}) \\ \boldsymbol{\Theta}_y &= (\boldsymbol{t}_{\theta_1+\omega},\ldots, \boldsymbol{t}_{\theta_N + \omega}) \end{aligned}\]

API Reference

Finding a Linear Model from Data

InvariantModels.findLinearModel — Function
findLinearModel(dataX, thetaX, dataY, thetaY, omega; ratio=0.5, steps=3)

The input arguments arguments

  • dataX is matrix $\boldsymbol{X}$. Each column corresponts to a data point.
  • dataY is matrix $\boldsymbol{Y}$.
  • thetaX is matrix $\boldsymbol{\Theta}_x$
  • thetaY is matrix $\boldsymbol{\Theta}_y$
  • ratio the ratio of data to be kept at the final step
  • steps the number of steps to take to find

Finds a linear model in the form of

\[\begin{aligned} \boldsymbol{x}_{k+1} &= \boldsymbol{A}(\theta_k) \boldsymbol{x}_k + \boldsymbol{b}(\theta_k) \\ \theta_{k+1} &= \theta_k + \omega \end{aligned}\]

and calculates the invariant torus represented by $\boldsymbol{K}$ from the equation

\[\boldsymbol{K}(\theta+\omega) = \boldsymbol{A}(\theta) \boldsymbol{K}(\theta) + \boldsymbol{b}(\theta)\]

The output is 

A, b, MK, XK
  • A is a 3-index array representing $\boldsymbol{A}(\theta)$
  • b is a 2-index array representing $\boldsymbol{b}(\theta)$
  • MK, XK represent to torus parametrisation $K$
source
InvariantModels.getgrid — Function
getgrid(n::Integer)

Returns a uniform grid on the interval $[0, 2 \pi)$ with $2n + 1$ number of points. The end point $2\pi$ is not part of the grid.

The number $n$ corresponds to the number of Fourier harmonics that can be represented on the grid.

source
InvariantModels.QPConstant — Type
M = QPConstant(dim_out, fourier_order, field::AbstractNumbers=ℝ)

Creates a vector valued periodic function with fourier_order Fourier harmonics. The parameters are

  • dim_out number of dimensions,
  • fourier_order number of Fourier harmonics,
  • field whether it is real valued field=ℝ or complex valued field=ℂ.
source
InvariantModels.QPPolynomial — Type
M = QPPolynomial(dim_out, dim_in, fourier_order, min_polyorder, max_polyorder, field::AbstractNumbers=ℝ; perp_idx=1:dim_in)

Creates a vector valued polynomial that also a periodic function of another variable.

\[ \boldsymbol{y} = P(\boldsymbol{x}, \theta)\]

The parameters are

  • dim_out number of output dimensions, i.e., the dimensionality of $\boldsymbol{y}$
  • dim_in number independent variables, i.e., the dimensionality of $\boldsymbol{x}$
  • fourier_order number of Fourier harmonics for the independent variable $\theta$
  • min_polyorder the lowest monomial order contained in polynomial $P$
  • max_polyorder the highest monomial order contained in polynomial $P$
  • field whether it is real valued field=ℝ or complex valued field=ℂ
  • perp_idx the indices of idependent variables in $\boldsymbol{x}$ that must appear at least linearly in each monomial. For example if we have 4 variables and perp_idx = [1 2], the monomial $x_3 x_4$ is not part of the polynomial, because they do not contains any of $x_1$ or $x_2$. However, $x_1 x_3$ is part of the polynomial, because $1$ is among the indices given by perp_idx.
source
InvariantModels.fromFunction! — Function
fromFunction!(M::QPPolynomial, W, fun)

Taylor expands Julia function fun to a polynomial QPPolynomial, whose structure is prescribed by M. Function $f$ = fun must have two arguments

\[ f : \mathbb{R}^n \times [0,2\pi) \to \mathbb{R}^m,\]

where $n$ is dim_in and $m$ is the dim_out parameter of polynomial M. The result is copied into W.

source
InvariantModels.Eval — Function
res = Eval(M::QPPolynomial, X, theta::Matrix, data::Matrix; cache=makeCache(M, X, theta, data))

Same as Eval!, except that the result is returned in res.

source
InvariantModels.Eval! — Function
Eval!(res, M::QPPolynomial, X, theta::Matrix, data::Matrix; cache=makeCache(M, X, theta, data))

Evaluate the polynomial M, X at multiple data points. The data points are given by theta and data. Each column of theta and data corresponds to a single data point. The result is copied into res.

The matrix theta contains the interpolation weights of Fourier collocation and not the actual values of $\theta$. Generally, these interpolation weights are pre-calculated, which avoids repeated calculation. The number of rows of $\theta$ is the same as $2 n + 1$, where $n$ is the number of Fourier modes resolved by polynomial M.

source
InvariantModels.mapFromODE — Function
mapFromODE(M::QPPolynomial, W, fun, par, omega, dt, step)

Creates a Poincare map from the vector field given by fun and places it into the polynomial M, X. The parameters are

  • M ,W the polynomial to contain the result.
  • fun the vector field in the form of fun(dx, x, p, theta), where x is the input state var, dx is the time-derivative $\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t}$, p is the parameter, which equals to par and theta is the phase of forcing on the interval $[0, 2\pi)$.
  • par is the parameter (vector) of fun.
  • omega is the forcing frequency. The independent variable t of fun is calculated as theta = $\omega t$, where $t$ is the independent variable (time).
  • dt is the time step of the integrator.
  • step is time period that the differential equation is solved for.
source
InvariantModels.toFourier — Function
MF, XF = toFourier(M::QPPolynomial, X)

Converts the polynomial M, X to another polynomial that contains Fourier coefficients as opposed to collocated values of the polynomial.

source

Semi-analytic Reduction to Invariant Manifolds

InvariantModels.QPMAPTorusManifold — Function
QPMAPTorusManifold(MP::QPPolynomial, XP, omega, SEL; 
                        threshold = 0.1, 
                        resonance = false, 
                        rom_order = max_polyorder)

The inputs are the following

  • MP, XP the Poincare map that describes the dynamics.
  • omega the phase shift at each time-step
  • SEL which spectral bundles to use for the calculation of the invariant manifold
  • threshold if a near resonance is greater than this value, it is ignored
  • resonance: when false the conjugate map is reduced to an autonomous system, even if there are near resonances. Otherwise all near resonances are accounted for in the conjugate map.
  • rom_order what is the highest order for which near resonances are taken into account. It defaults to the order of the Poincare map MP, XP.

This function does three calculations

  1. Calculates the invariant torus using Newton's method. The startnig iteration is at the origin.
  2. Calculates the linear system about the invariant torues and decomposes it into invariant vector bundles. The bundles are ordered by the corresponding spectrum. The bundles with the greatest spectral radius appear first.
  3. Given the invariant bundles about the torus, calculate the invariant manifold tangent to the selected vector bundles. The selection is done by SEL.

Return values are:

MK, XK, MSn, XSn, MWn, XWn, MSd, XSd, XWre
  • MK, XK represent the invariant torus
  • MSn, XSn the normal form of the conjugate dynamics on the invariant manifold
  • MWn, XWn the manifold immersion in the coordinate system of the invariant vector bundle, selected by SEL
  • MSd, XSd the full system MP, XP transformed to the coordinate system of the vector bundle decomposition
  • XWre the transformation from the vector bundles to the original coordinate system
source
InvariantModels.QPODETorusManifold — Function
QPODETorusManifold(MP::QPPolynomial, XP, omega, SEL; 
                        threshold = 0.1, 
                        resonance = false, 
                        rom_order = max_polyorder)

This is the equivalent version of QPMAPTorusManifold for vector fields.

The inputs are the following

  • MP, XP the vector field that describes the dynamics.
  • omega is the forcing frequency
  • SEL which spectral bundles to use for the calculation of the invariant manifold
  • threshold if a near resonance is greater than this value, it is ignored
  • resonance: when false the conjugate map is reduced to an autonomous system, even if there are near resonances. Otherwise all near resonances are accounted for in the conjugate map.
  • rom_order what is the highest order for which near resonances are taken into account. It defaults to the order of the Poincare map MP, XP.

This function does three calculations

  1. Calculates the invariant torus using Newton's method. The startnig iteration is at the origin.
  2. Calculates the linear system about the invariant torues and decomposes it into invariant vector bundles. The bundles are ordered by the corresponding spectrum. The bundles with the greatest spectral radius appear first.
  3. Given the invariant bundles about the torus, calculate the invariant manifold tangent to the selected vector bundles. The selection is done by SEL.

Return values are:

MK, XK, MSn, XSn, MWn, XWn, MSd, XSd, XWre, Lambda
  • MK, XK represent the invariant torus
  • MSn, XSn the normal form of the conjugate dynamics on the invariant manifold
  • MWn, XWn the manifold immersion in the coordinate system of the invariant vector bundle, selected by SEL
  • MSd, XSd the full system MP, XP transformed to the coordinate system of the vector bundle decomposition
  • XWre the transformation from the vector bundles to the original coordinate system
  • Lambda is the linear part of the system about the invariant torus in the coordinate system of the vector bundles. This is for debug purposes only, it should be a diagonal matrix for each collocation point.
source

Manifolds via Invariant Foliations

InvariantModels.QPGraphStyleFoliations — Function
QPGraphStyleFoliations(MP::QPPolynomial, XP, omega, SEL; dataScale=1, resonance=false, threshold=0.1, rom_order=max_polyorder)

The purpose of this function is to calculate two complementary invariant foliations about an invariant torus and then extract the invariant manifold that is the zero level set of the second foliation. The invariant manifold has the same dynamics as the first foliation.

The inputs are the following

  • MP, XP the Poincare map that describes the dynamics.
  • omega the phase shift at each time-step
  • SEL which spectral bundles to use for the calculation of the invariant manifold
  • dataScale is the scaling factor that scales the coordinate system of the invariant manifold. This scaling is used when we would like to compare the invariant manifold identified from data that has been scaled. Using the same scaling factor makes this comparison possible.
  • threshold if a near resonance is greater than this value, it is ignored
  • resonance: when false the conjugate map is reduced to an autonomous system, even if there are near resonances. Otherwise all near resonances are accounted for in the conjugate map.
  • rom_order what is the highest order for which near resonances are taken into account. It defaults to the order of the Poincare map MP, XP.

Return values are:

MRf, XRf, MW, XW, MRt, XRt, MUV, XUVflat, MSt, XSt, MVU, XVUflat
  • MRf, XRf the conjugate dynamics of the first invariant foliation (SEL) in normal form coordinates and with Fourier coefficients
  • MW, XW the immersion of the invariant manifold in the original coordinate system

The outputs that could be used to set initial conditions to a foliation

  • MRt, XRt the conjugate dynamics of the first invariant foliation (SEL) when the foliation is parametrised graph-style
  • MUV, XUVflat
  • MSt, XSt the conjugate dynamics of the second invariant foliation when the foliation is parametrised graph-style
  • MVU, XVUflat
source

Frequencies and Damping Ratios

InvariantModels.MAPFrequencyDamping — Function
MAPFrequencyDamping(MW::QPFourierPolynomial, XW, MR::QPFourierPolynomial, XR, amp_max; output = Diagonal(I,dim_out))

Calculates corrected frequencies and damping ratios for a centre type equilibrium. The inputs are

  • MW, XW is the manifold immersion $\boldsymbol{W}: \mathbb{R}^2 \times [0,2\pi) \to \mathbb{R}^n$
  • MR, XR is the conjugate dynamics $\boldsymbol{W}: \mathbb{R}^2 \times [0,2\pi) \to \mathbb{R}^2$
  • amp_max is the maximum amplitude to be considered
  • output is a matrix $\boldsymbol{M}$ that is used to calculate the instantaneous amplitude by pre multiplying the immersion: $\boldsymbol{M} \boldsymbol{W}$.

The inputs are such that they satisfy the invariance equation

\[\boldsymbol{W}\left(\boldsymbol{R}\left(\boldsymbol{z},\theta\right),\theta+\omega\right)=\boldsymbol{F}\left(\boldsymbol{W}\left(\boldsymbol{z},\theta\right),\theta\right),\]

where $\boldsymbol{F}$ is the map of the original dynamics. Given that the invariant manifold is in a normal form, we have the following structure to the conjugate map

\[\boldsymbol{R}\left(\boldsymbol{z}\right)=\begin{pmatrix}s\left(z,\overline{z}\right)\\ \overline{s}\left(z,\overline{z}\right) \end{pmatrix}\]

Using a the transformation $z=\rho\left(r\right)\exp\left(i\beta+i\alpha\circ\rho\left(r\right)\right)$ we can write the manifold immersion as

\[\widehat{\boldsymbol{W}}\left(r,\beta,\theta\right)=\boldsymbol{W}\left(\rho\left(r\right)\exp\left(i\beta+i\alpha\circ\rho\left(r\right)\right),\rho\left(r\right)\exp\left(-i\beta-i\alpha\circ\rho\left(r\right)\right),\theta\right)\]

so that the invariance equation holds

\[\widehat{\boldsymbol{W}}\left(R\left(r\right),\beta+T\left(r\right),\theta+\omega\right)=\boldsymbol{F}\left(\widehat{\boldsymbol{W}}\left(r,\beta,\theta\right),\theta\right),\]

where $R$ and $T$ are calculated using the same transformation.

The unknown functions $\rho$ and $\alpha$ are calculated such that

  1. the amplitude of a tori as a function of $r$ is the same as $r$ and
  2. there is zero phase shift between any two tori with different amplitudes with respect to parameter $\beta$.

Return values are

T, R_r, rho, alpha
  • T is the function $T : [0,\infty) \to \mathbb{R}$
  • R_r is the function $r \mapsto R(r) / r$
  • rho is the same as $\rho$
  • alpha is the same as $\alpha$

The instantaneous frequency is calculated as

\[\omega(r) = \frac{T(r)}{\Delta t},\]

the instantaneous damping is calculated as

\[\zeta(r) = -\frac{\log r^{-1} R(r)}{T(r)},\]

where $\Delta t$ is the sampling period.

source
InvariantModels.ODEFrequencyDamping — Function
ODEFrequencyDamping(MW::QPFourierPolynomial, XW, MR::QPFourierPolynomial, XR, amp_max; output = Diagonal(I,dim_out))

The inputs and outputs are the same as MAPFrequencyDamping, except that the interpretation of the input/output is slightly different. The input satisfies the invariance equation

\[D_1\boldsymbol{W}\left(\boldsymbol{z},\theta\right) \boldsymbol{R}\left(\boldsymbol{z},\theta\right) + D_1\boldsymbol{W}\left(\boldsymbol{z},\theta\right) \omega = \boldsymbol{F}\left(\boldsymbol{W}\left(\boldsymbol{z},\theta\right),\theta\right),\]

where $\boldsymbol{F}$ is a vector field.

From the output, the instantaneous frequency is

\[\omega(r) = T(r),\]

the instantaneous damping is calculated as

\[\zeta(r) = -\frac{r^{-1} R(r)}{T(r)}.\]

source

Generating Data

InvariantModels.generate — Function
generate(NDIM, rhs!, maxIC, nruns, npoints, fourier_order, omega, Tstep, A, XWre)

Creates a data set from a differential equation, represented by rhs!. An invariant torus is found by observing a convergent simulation to the torus. The initial conditions are sampled from an ellipsoid about the torus. The directions of the ellipsoid are given by the linear transformation XWre. For each column of XWre the maximum amplitude is specified by the corresponding element of the vector maxIC.

The input parameters are

  • NDIM the dimensionality of the system
  • rhs! the vector field in the form of fun(dx, x, p, theta), where x is the input state var, dx is the time-derivative $\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t}$, p is the parameter, which equals to par and theta is the phase of forcing on the interval $[0, 2\pi)$.
  • maxIC is a vector of size NDIM, containing the maximum simulation amplitudes in each direction of XWre
  • nruns the number of simulation runs
  • npoints the length of the trajectory for each simulation run
  • fourier_order the number of Fourier harmonics to consider when calculating $\boldsymbol{\Theta}_x$, $\boldsymbol{\Theta}_y$
  • omega the forcing frequency, as in $\dot{\theta} = \omega$
  • Tstep the time step between two samples of the trajectory, denoted by $\Delta t$
  • A the parameter p of the vector field
  • XWre is the transformation that defines the axes of the ellipsoid about the invariant torus

The function returns

dataX, dataY, thetaX, thetaY, thetaNIX, thetaNIY
  • dataX is matrix $\boldsymbol{X}$. Each column corresponts to a data point.
  • dataY is matrix $\boldsymbol{Y}$.
  • thetaX is matrix $\boldsymbol{\Theta}_x$
  • thetaY is matrix $\boldsymbol{\Theta}_y$
  • thetaNIX is the vector $(\theta_1, \ldots, \theta_N)$
  • thetaNIY is the vector $(\theta_1 + \omega \Delta t, \ldots, \theta_N + \omega \Delta t)$

space

source
InvariantModels.generateMap — Function
generateMap(MF, XF, MK, XK, maxIC, nruns, npoints, fourier_order, omega, Tstep)

Similar to generate, except that it simulates a discrete-time system.

The inputs differ from generate. Instead of the vector field, the following must be specified

  • MF, XF the discrete-time map of the system
  • MK, XK the invariant torus arbout which the system is simulated.

There is no way the specify an ellipsoid for initial conditions, they are taken from a sphere of radius maxIC.

source

Identfying Invariant Foliations from Data

InvariantModels.QPPreProcess — Function
QPPreProcess(K, A, omega, thetaX, dataX, thetaY, dataY, sel; Tstep=1.0, maxAmp=Inf, data_ratio=1.0, preId=[])

This function is called after a linear model has been identified using findLinearModel. The following calculations are carried out

  1. decompose the linear dynamics into invariant vector bundles.
  2. project the data onto two invariant vector bundles, one that corresponds to the spectral circles selected by sel and the the remaining bundle. In this coordinate system the approximate linear model represented by A becomes block-diagonal.

The input arguments arguments

  • K approximate invariant torus
  • A approximate linear map about the torus
  • dataX is matrix $\boldsymbol{X}$. Each column corresponts to a data point.
  • dataY is matrix $\boldsymbol{Y}$.
  • thetaX is matrix $\boldsymbol{\Theta}_x$ These are the interpolations vectors for Fourier collocation corresponding to $\theta_k$.
  • thetaY is matrix $\boldsymbol{\Theta}_y$
  • sel indicates the selected vector bundles.
  • Tstep is the sampling period $\Delta t$
  • maxAmp is the maximum amplitude of the data along the selected vector bundle
  • data_ratio $\in [0,1]$ the proportion of the data to be kept and the rest is removed so that it is uniformly distributed along in the neighbourhood of the selected vector bundles.
  • preId if non-empty, this replaces the data filtering process.

The output is

thetaTX, dataTX, thetaTY, dataTY, dataScale, id, R1, S1, W1
  • thetaTX, dataTX, thetaTY, dataTY are the filtered and scaled data. The approximate invariant torus is subtracted and the data is transformed into the coordiante system of the approximate vector bundles calculated from A.
  • dataScale is the scaling factor that relates the input and output dataX[:,id] = dataScale * dataTX.
  • id the indices of the input that are brought to the output
  • R1 is the invariant part of the the linear system A that is invariant under the selected sel vector bundles
  • S1 is the invariant part of the the linear system A that is invariant under the not selected (setdiff(1:data_dim, sel)) vector bundles
  • W1 the inverse transformation that brings that data back into the physical coordinate system.
source
InvariantModels.QPCombinedFoliation — Type
QPCombinedFoliation(dim_data, dim_latent, fourier_order, R_order, U_order, S_order, V_order, R1, S1, MK, XK; sparse=false)

Creates the data structures to store a pair of foliations that might be interconnected through scaling. The invariance equations represented by a Combined Foliation are

\[\begin{aligned} \boldsymbol{R}\left(\boldsymbol{U}\left(\boldsymbol{x},\boldsymbol{\theta}\right),\boldsymbol{\theta}\right) &= \boldsymbol{U}\left(\boldsymbol{F}\left(\boldsymbol{x},\boldsymbol{\theta}\right),\boldsymbol{\theta}+\boldsymbol{\omega}\right)\\ \boldsymbol{S}\left(\boldsymbol{V}\left(\boldsymbol{x},\boldsymbol{\theta}\right),\boldsymbol{\theta}\right) &= \boldsymbol{V}\left(\boldsymbol{F}\left(\boldsymbol{x},\boldsymbol{\theta}\right),\boldsymbol{\theta}+\boldsymbol{\omega}\right) \end{aligned}\]

The input parameters are

  • dim_data system dimension
  • dim_latent the dimension of the reduced order model
  • fourier_ordernumber of Fourier harmonics used in discretisation
  • R_order the polynomial order of function $\boldsymbol{R}$
  • U_order the polynomial order of function $\boldsymbol{U}$
  • S_order the polynomial order of function $\boldsymbol{S}$
  • V_order the polynomial order of function $\boldsymbol{V}$
  • R1 the initial linear part of function $\boldsymbol{R}$
  • S1 the initial linear part of function $\boldsymbol{S}$
  • MK, XK the invariant torus
  • sparse use compressed polynomials (default = false)

The output is the data structure

MCF, XCF
source
InvariantModels.makeCache — Function
makeCache(M::QPCombinedFoliation, X, thetaX, dataTX, thetaY, dataTY)

Creates a cache object for a combined foliation. The cache stores intermediate results that are not re-calculated between various steps of the optimisation. This improves performance significantly.

Input:

  • M, X the combined foliation generated by QPCombinedFoliation
  • dataTX is matrix $\boldsymbol{X}$. Each column corresponds to a data point.
  • dataTY is matrix $\boldsymbol{Y}$.
  • thetaX is matrix $\boldsymbol{\Theta}_x$ These are the interpolations vectors for Fourier collocation corresponding to $\theta_k$.
  • thetaY is matrix $\boldsymbol{\Theta}_y$

The output is a cache object cache.

source
InvariantModels.updateCache! — Function
updateCache!(cache::FCache, M, X, thetaX, dataTX, thetaY, dataTY)

The same argument as makeCache, expect that it updates cache.

source
updateCache!(cache::FCache, M, X, thetaX, dataTX, thetaY, dataTY, sel)

The same input argument as makeCache, expect that it updates cache only for the component of the foliation given by sel.

source
InvariantModels.QPOptimise — Function
QPOptimise(MCF::QPCombinedFoliation, XCF, thetaTX, dataTX, thetaTY, dataTY;
    maxit=24, 
    gradient_ratio=2^(-7), gradient_stop=2^(-29), 
    steps=400, 
    name="default",
    cache=makeCache(MCF, XCF, thetaTX, dataTX, thetaTY, dataTY),
    omega=1.0, Tstep=1.0, dataScale=1.0, dataId=1:size(thetaTX, 2), radii=to_zero(XCF), skipRS=false)

Use optimisation to fit the parameters of the combined foliation to data.

Input parameters are

  • MCF, XCF the combined foliation
  • thetaTX, dataTX, thetaTY, dataTY the data
  • maxit maximum number of iterations within one trust-regions subproblem solution
  • gradient_ratio the subproblem solution stops once the norm of the gradient is reduced by this factor
  • gradient_stop the subproblem solution stops when the norm of the gradient is below this value
  • steps the number of subproblem solutions as we iterate over the components of the combined foliation
  • name the name of the problem, used for logging purposes
  • cache is created by makeCache
  • omega the forcing shift-angle
  • Tstep the sampling period of the data $\Delta t$
  • dataScale the scaling factor used to scale the data
  • dataId the indices of the data points from the original data set
  • radii initial trust-regions radii for each component of the combined foliation
source
InvariantModels.QPPostProcess — Function
QPPostProcess(MCF::QPCombinedFoliation, XCF, omega, resonance=false, threshold=0.1)

Extracts the invariant manifold from the combined foliation together with its dynamics. The result is in a normal form with mixed power series Fourier series form.

Paramaters:

  • MCF, XCF the combined foliation.
  • W1 the inverse trabsformation produced by QPPreProcess
  • omega the forcing shift-angle.
  • resonance true if to account for non-autonomous internal resonances. When false, the result will be an autonomous normal form.
  • threshold if the distance between the two side of the non-resonance criterion is less than threshold, it counts as an internal (or external) resonance

Returns:

MSn, XSn, MFW, XFWoWdoWn, MW, XW
  • MSn, XSn the conjugate dynamics on the manifold in normal form
  • MFW, XFWoWdoWn manifold immersion in the physical coordinate system and adapted to correspond to the normal form
  • MW, XW manifold immersion in the coordinate system of the combined foliation (not physical coordinates)
source
InvariantModels.EvalU — Function
EvalU(M::QPCombinedFoliation, X, thetaX, dataTX; cache)

Evaluates the submersion $\boldsymbol{U}$ at data points thetaX, dataTX.

source
InvariantModels.EvalV — Function
EvalV(M::QPCombinedFoliation, X, thetaX, dataTX; cache)

Evaluates the submersion $\boldsymbol{V}$ at data points thetaX, dataTX.

source
InvariantModels.Loss — Function
Loss(M::QPCombinedFoliation, X, thetaX, dataTX, thetaY, dataTY; cache)

Returns the loss function of the combined foliation M, X at data points thetaX, dataTX, thetaY, dataTY.

source
InvariantModels.ResidualU — Function
ResidualU(M::QPCombinedFoliation, X, thetaX, dataTX, thetaY, dataTY; cache)

Returns the residual of the invariance equation for $\boldsymbol{U}$ divided by the norm of $\boldsymbol{V}$.

source
InvariantModels.ResidualLoss — Function
ResidualLoss(M::QPCombinedFoliation, X, thetaX, dataTX, thetaY, dataTY; cache)

Returns the loss function only for the component with $\boldsymbol{U}$ of the combined foliation M, X at data points thetaX, dataTX, thetaY, dataTY.

source
InvariantModels.QPKernel — Function
QPKernel(MF::QPCombinedFoliation, XF)

Returns the root of the combined foliation in the form of MK, XK, which is the same as the invariant torus.

source
InvariantModels.copyMost! — Function
copyMost!(MCFD::QPCombinedFoliation, XCFD, MCFS::QPCombinedFoliation, XCFS)

Make a best-effort attempt to copy the combined foliation MCFS, XCFS into MCFD, XCFD.

source
InvariantModels.SetTorus! — Function
SetTorus!(MCF::QPCombinedFoliation, MK::QPConstant, XK)

Set the location of the torus MK, XK for the combined foliation MCF.

source
InvariantModels.bincuts — Function
bincuts(long_amp, trans_amps, max_amp, max_points, nbins; exponent=1)

This function is used to filter out data points that are far from an invariant manifold.

  • long_amps vector of amplitudes along the invariant manifold
  • trans_amps vector of distances from the invariant manifold
  • max_amp discard points that have greater amplitude than max_amp
  • max_points the maximum number of points to return
  • nbins the number of bins along which the distribution is made uniform
  • exponent used to calculate bin sizes. The bins fall between the nodes of range(0, max_amp, length=nbins + 1) .^ exponent.

Returns the indices of long_amps that create a uniform distribution with respect to the bins and only contain the smallest trans_amps values from those bins.

source
InvariantModels.to_zero — Function
to_zero(x::ArrayPartition)

Creates an ArrayPartition the same structure as x expect that each component is a single element array with zero value.

source