Predicting the bounds of large chaotic systems using low-dimensional manifolds

Predicting extrema of chaotic systems in high-dimensional phase space remains a challenge. Methods, which give extrema that are valid in the long term, have thus far been restricted to models of only a few variables. Here, a method is presented which treats extrema of chaotic systems as belonging to discretised manifolds of low dimension (low-D) embedded in high-dimensional (high-D) phase space. As a central feature, the method exploits that strange attractor dimension is generally much smaller than parent system phase space dimension. This is important, since the computational cost associated with discretised manifolds depends exponentially on their dimension. Thus, systems that would otherwise be associated with tremendous computational challenges, can be tackled on a laptop. As a test, bounding manifolds are calculated for high-D modifications of the canonical Duffing system. Parameters can be set such that the bounding manifold displays harmonic behaviour even if the underlying system is chaotic. Thus, solving for one post-transient forcing cycle of the bounding manifold predicts the extrema of the underlying chaotic problem indefinitely.


Introduction
Coupled systems of non-linear ordinary differential equations (ODEs) are ubiquitous in scientific literature and appear as the end result of efforts to mathematically model phenomena from practically every branch of science. The desire for quantitative accuracy often leads to such models containing large numbers of variables, for example, as a consequence of a highly resolved finite element mesh, an in-silico model of human tissue-be it heart, tumour or brain -containing many simulated cells, or an electrical circuit with many nodes. Systems of ordinary differential equations with non-linear terms have the potential to display chaotic behaviour, in which case their solutions evolve on "strange" attractors of fractal dimension.
The ability to precisely predict the extrema of large chaotic systems has obvious utility across many fields. While bifurcation analysis serves to characterise the qualitative behaviour of a system, extrema prediction adds quantitative knowledge. Perhaps the clearest benefit is where state variables have critical values, e.g., in predator-prey systems where reaching a a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 possible future states of a system relatively narrowly in phase-space and that they can be applied to high-D chaos in high-D phase space. The obvious drawback is their local and consequently short-term nature. The effort towards determining quantitatively accurate long-term global bounds of chaotic systems has mainly been focused on obtaining closed-form solutions to specific systems on an ad-hoc basis. As the information required to store a manifold is exponentially dependent on dimension, existing methods of bounding strange attractors break down in high-D phase space. Efforts to bound chaotic systems have, so far, mainly focused on systems of three variables [40,41] with attractor dimension close to phase space dimension. A more generalised method for the bounding of low-D chaotic systems was proposed in [42,43], which took the approach of defining localising sets with bounding manifolds defined as polynomials.
The present method is applicable in situations when the phase-space is high-D, but the attractor is low-D. In this situation it is computationally feasible to define a long-term bound for the system. Although study of high-D chaos and the factors impacting on attractor dimension has taken place since the 1980s [39,[44][45][46], as discussed above, many natural processes, even if high-D in terms of phase space, generate low-D attractors. It has been and still is commonplace to reduce large systems to the minimum dimension necessary for qualitative accuracy [5], which in turn tends to be defined as preservation of the topology of the attractor [9]. However, qualitative methods sacrifice quantitative accuracy. And quantitative accuracy may be crucial. Even periodic systems, are often simulated by very large quantitatively accurate high-D models. The study of dynamic phenomena with linear finite element models is an example of this.
The method presented here finds its niche in situations where model reduction is unacceptable and the extrema of system variables are important. It is deliberately designed to be intuitive and easy to apply, i.e. "wrap around" an existing system. The aim is that, once applied and proven on attractors of higher (but still low) dimension, it can be used by non-experts in chaos and non-linear dynamics in a plug and play fashion.

Results
Two manifolds are defined in the following: The auxiliary manifold (AM) and the bounding manifold (BM). The AM is an M-dimensional manifold, where M is an integer. It is embedded in the phase space of the system to be bounded and, in turn, the BM is embedded on the AM. The hope is that the AM and the BM display stationary harmonic behaviour in time, such that they only need to be calculated for a couple of post-transient forcing cycles.

The bounding-and auxiliary manifolds
The manifolds are non-autonomous extensions of the concept of invariant and inertial manifolds and draw inspiration from existing literature [47,48]. However, the following does not adhere strictly to the quite specific definitions and nomenclature of inertial manifold theory. Firstly, inertial manifolds are most often employed as a means to directly solve the original equations on the manifold. They allow the reduction of dimensionality from either infinity (for partial differential equations) or high-finite dimension (ODEs with many degrees of freedom). Here, no attempt is made at solving the underlying system of ODEs on the manifolds. The AM serves the purpose of reducing the dimensionality of the BM-not that of the underlying problem. Secondly, inertial manifold theory rests on the premise that solutions in full phase-space approach a solution on the manifold. The requirement here is simply that solutions approach the manifold, which is less restrictive. So, although the AM draws from inertial manifold theory, it is named as it is to emphasize the differences and avoid confusion.
If all unstable directions are contained on the AM, solutions that do not start out on the manifold, will decay exponentially onto it and, post-transiently, an attractor will lie on the AM. Let c denote the bounding manifold. Then, the premise of the method is that where A is an attractor, B is the auxiliary manifold and C is an M-manifold with boundary. Thus, c is a closed (M − 1)-manifold. Note that "bounding manifold" is formally a (potential) misnomer as extrema of the system may lie on interior points (and not on the BM itself). Further, this necessitates the post hoc calculation of interior points, which is discussed below.
The AM is embedded in the full phase space x 2 R N of the system of ODEs: As it shall become apparent in the following, it is a requirement that f is smooth. Furthermore, it is a requirement that M > d, where d is the (fractal) dimension of the attractor. The most obvious choice is to select M as the nearest integer greater than d. As mentioned above, the discipline of bounding chaotic systems has previously focused on N − 1 dimensional boundaries to systems with N dimensional phase spaces [40,42,43]. Such approaches will fail rapidly as N increases to that of even modest sized problems, since the computational cost of operating on a discretised manifold depends exponentially on its dimension. Here, it is exploited that the strange attractor is normally a low-D entity in a high-D phase space, i.e, that d < M < <N. So, the BM, c(t), need not be more than M − 1 dimensional. If c(t) is harmonic in time, the solution can be stopped after a couple of post-transient forcing cycles. Fig 1 shows this scenario for M = 2. For the method to work, the AM must be invariant under the flow, meaning that solutions that start on it, stay there. The manifold is allowed to vary with time. So, labelling the manifold "invariant" may cause confusion. Nevertheless, this is the naming convention widely adopted in literature and time dependent inertial manifolds are well known [49]. The key requirement is that the manifold displays regular motion. Clearly, if the manifold itself is chaotic, hence the BM also, the method has achieved very little.
To introduce key concepts, we consider the case of a 2-D AM in a R N phase space and then generalise the ideas to M-D manifolds in R N phase spaces. Fig 2 shows a schematic of the BM and AM. It also introduces the basis vectors of the AM (w 1 ,w 2 ), the inward normal u 1 and the BM unit tangent vector u 2 . Furthermore, the figure shows an off-manifold point which is projected along a stable direction onto the manifold. If certain criteria are met, it can be proven that the off-manifold point converges exponentially with a solution on the manifold [9,48,50,51]. However, what matters is that off-manifold points approach the manifold asymptotically, end up somewhere on it and stay there, which is uncontroversial. Here, the notion of points being "on" the manifold is taken to include points infinitesimally close to it, which is the case with asymptotically approaching points. Obviously, once the distance drops below machine precision, there is no practical distinction.
The motion of the BM is the sum of the tangent on-manifold motion and the orthogonal off-manifold motion. The latter is readily determined through Gram-Schmidt processing of the flow. The expression for BM motion then becomes: where f(x, t) = f(c(m), t) is defined in Eq (2), _ cðm; tÞ is the rate of the BM, d i are the on-manifold velocities and u i is a local orthonormal basis spanning the tangent space of the AM. For AMs of dimension 3, 4 and so on, Eq (3) is easily expanded with d 3 , u 3 , d 4 , u 4 , etc. That is, generally The flow and the vectors are evaluated locally along the length of the BM and depend on time, i.e. Initially they lie as a cloud in phase space, but they rapidly contract along stable directions until they converge to the AM. As shall become apparent, it is sufficient to calculate the tangent space of the AM to orientate the BM c(m, t). This is done with a spectral method.
The fundamental idea behind the method is that the locally unstable directions, the directions of stretch, are contained in the on-manifold terms ("k"). If that is the case, the remaining offmanifold motion ("?") will consist only of exponential decay onto the manifold.
The nomenclature is such that u 1 is the inward normal (tangent to the AM and orthogonal to the BM) with all other vectors spanning the space tangent to the BM (tangent to both the AM and the BM). The inward motion of the boundary is taken to follow the heuristic expression: where β is the linear impedance and κ = κ(c(m, t)) is the signed dimensionless curvature signal passed through a low-pass filter with time constant τ: The equation of motion Eq (6) automatically bounds the underlying ODE and filters away high frequency oscillations. It will straighten the BM locally until the curvature term vanishes, at which point the BM is at equilibrium. Since the BM is a closed curve, the global net effect of the curvature term is to contract the BM. The responsiveness of the BM is determined by the values of β and τ. Although the expression is non-linear in d 1 , it has the merit of being local, so does not require the (iterative) solution of a system of non-linear equations. A spatially coupled formulation, e.g., based on constrained optimization and a Lagrange multiplier field, would come with this drawback. The present study considers 1-D BMs, so the measure of curvature is given. For higher-dimensional BMs (2-D surfaces and beyond), a scalar quantity such as the mean curvature might be used. The motion of the BM is discretised using Hermite finite elements. These have the sufficient order to allow the calculation of curvature. Axial motion. The axial motion of the BM seeks to space nodes evenly along the BM. This is no more than an adaptive meshing exercise. For the sufficiently well discretised BM, the axial motion has no appreciable impact on the shape of c. The below expressions are valid for a 1-D BM, so a general (M − 1)-D version is needed. Consider the rate of change of elastic energy with the addition of a dissipative term Now, applying Eq (3) and requiring stationarity with respect to the d 2 field, one gets the expression where prime denotes differentation with respect to m. As with the d 1 field described above, the d 2 field is expanded in terms of Hermite shape functions and the finite element method is used to solve Eq (9) for nodal values of d 2 . The details are given in the methods section.

Interior points
Unlike in some other methods from the literature, the interior of the manifold is not necessary for the determination of manifold orientation. Despite this, it must still be considered, as an extremum of a state variable might lie on the interior and not on the BM. One must balance the added computational cost of sampling the interior against the need for a high number of samples. The equation of motion for the interior points is a compromise between simplicity and the desire for evenly spaced points and could be defined in a variety of ways. What is unequivocal though, is that interior points must ride on the tangent space associated with the AM. The orientation at interior points is determined in the same way as on the BM. A practical and fairly straight forward way of ensuring evenly spaced points throughout a simulation, where the BM may change shape significantly, is to connect the interior points to their neighbours and the BM with linear truss elements of zero initial length. Then, point positions are updated dynamically as nodes of truss elements. For the purposes outlined here, this is nothing more than a simple, and probably sub-optimal, adaptive meshing technique. The advantage of truss elements is that they are readily defined in higher dimensions as their length can be easily calculated in N-dimensional Cartesian coordinates. Thus, they are practical for point spacing in phase spaces above 3-D. The equation of motion for interior points is identical in form to Eq (3), except that now d 1 and d 2 are determined by force summations of the truss elements. That is, the interior points, i.e., truss element nodes, are restricted to move on the AM onto which all truss element forces, accelerations and velocities are projected. The truss element normalised stiffness and damping, k t /m t and d t /m t , respectively, are given in Table 1. The motion of truss element nodes on the BM c is prescribed to follow the BM in a one-way coupling fashion, so the motion of the BM is not influenced by the interior nodes. This means that the motion of the BM can be computed first, and interior point sampling performed in postprocessing.
The equation of motion for an interior point, where d i are determined as and h 2 R N is the resultant mass-normalised force from trusses pulling on an interior point.

Example: The duffing oscillator
The Duffing oscillator is among the simplest systems known to exhibit chaos, yet still finds use for a wide variety of applications including, e.g., nano-mechanical systems [52]. Consider a system consisting of a Duffing oscillator in series with a number of linear spring-mass-damper oscillators. The linear oscillators add to the phase space of the system, and quantitatively affect Table 1. Model parameters for the simulations. The parameter τ is the time constant for curvature filtering, β is the proportional impedance on curvature, β 2 is the impedance on axial motion, k t , d t and m t are the stiffness, damping and mass of interior point truss elements, α is the non-linear Duffing stiffness, k 1 , d 1 and m 1 are the linear stiffness, damping and mass of the Duffing oscillator, f 0 is the forcing term amplitude, ω is the forcing term frequency, k 2 , d 2 and m 2 are the stiffness, damping and mass of the linear oscillators, Δt is the time step for the time integration scheme (BM and interior point position update), Δt λ is the time step for eigenvalue calculation and orientation update, K is the number of nodes (and elements) used to discretise c, and K G is the number of Gauss points for numerical integration in space along c.
All the solution to the equations, but do not add to the dimension of the attractor or influence its topology. If the duffing oscillator is assigned the first two degrees of freedom, the system of equations can be written as where f 0 cos(ωt) is a harmonic forcing term, k 1 , d 1 , α and m 1 are the parameters of the Duffing oscillator and all the linear oscillators are assigned the identical parameters k 2 , d 2 and m 2 . The governing equations are on the form of Eq (2) so can be plugged into the method as described above. In order to reduce the mesh refinement and time step requirements associated with curvature and stretch along fast, linear terms, the equations are transformed as ., x N ] T = [y 1 , y 2 , y 3 /10, . . ., y N /10] T and transformed back in post-processing. This transformation is beneficial and possible because it is known a-priory that the non-linearity of the system appears in the first two equations-thus "compressing" all other variables and effectively increasing mesh refinement in the non-linear directions. For some systems, all variables might appear in non-linear terms, or it may not be known which ones that do. Then, scaling becomes futile, potentially necessitating finer meshes and shorter time steps. Parameters for all simulations are shown in Table 1. In general terms, a fast system requires a fast BM. So a system containing high frequencies and rapid changes in J will require low impedances on the BM, lower time constants for the first order filter and, of course, smaller time steps. In Fig 3 is shown the convergence of the AM with twenty solutions for parameter set II in Table 1. The phase space is 42-D, and only the (y 1 , y 2 , y 42 ) projection is shown, so more information is contained in the simulation than can be displayed on the figure. Note that the solutions converge to the AM and stay there. The figure also shows the truss elements connecting interior points to each other and to the BM. A video of a similar process for parameter set III is shown in S1 File. Fig 4 shows snapshots of projections of the BM onto a selection of coordinate surfaces for parameter set III in Table 1. The projections vary considerably in shape and size, indicating differences in phase and magnitude, respectively, of the harmonic components of the state variable motion. A video of a similar process for parameter set III is shown in S2 File.  Fig 3). The motion of the AM and the BM is calculated for five forcing cycles and the subsequent 95 cycles are copies of the fifth, which is post-transient. Twenty time series are calculated for 100 forcing cycles in order to validate the BM and it is seen that they all remain within the predicted values. The BM is noted to be reasonably tight and to be following the low frequency dynamics of the system. As one would expect, not all variables in the system influence the BM equally. The geometric nature of the method means that for the BM to contract, curvature must exist in the coordinate direction of the relevant variable. With the proposed equation of motion, the BM will tend towards being convex everywhere. The approach works well with the here presented Duffing-type system, but there is no proof that it will in general, so future modifications to the equation of motion Eq (6) producing a locally concave BM might be beneficial.

Example: The parametric duffing oscillator
To illustrate the method's ability to deal with situations where the Jacobi matrix depends explicitly on time, a parametrically excited version of the duffing equation is considered. The equations for this system are given in Eq (13). Notice that the non-linear term now varies explicitly with time.

Local and instantaneous stretch and compression
Let us explore the behaviour of two infinitesimally nearby trajectories in a smooth, non-linear and non-autonomous flow. Due to smoothness, the state can be represented by its Taylor expansion. Consider a perturbation in time by Δt and in space by a v, where a is a scalar and v 2 R N is a unit vector: where Δt = t − t 0 and O(3) represents terms of order 3 and above. Fig 10 shows a state, x 1 , perturbed forward in time and another, x 2 , perturbed in space and time. Representing the perturbed states x 1 and x 2 as their Taylor expansions and substituting _ x ¼ f we obtain where f = f(x 0 , t 0 ). Now, from Eq (15) and recognising that @x @a ¼ v, the rate of stretch is readily obtained as lim a!0 It is easily seen that Where J is the Jacobi matrix. If v is an eigenvector with a real eigenvalue, the local stretch is parallel to it. If it is complex, the stretch will lie in the space spanned by v and its complex conjugate (since the flow is real). Thus, imagnining for a moment a full spectral decomposition, it is easily seen that, locally, the only directions of positive stretch are those associated with eigenvalues with a positive real part. Nearby solutions will converge exponentially in all other directions-i.e. contract. This holds for either autonomous or non-autonomous systems as no assumptions have been made in regards to the time dependence of _ xðx; tÞ apart from the requirement of smoothness. So, as is evident, the spectrum plays a pivotal role in manifold orientation, even for non-autonomous systems.
Although the above is not a novel result, its implications justify its repetition here. Chaos is a global phenomenon and eigenvalue decompositions cannot be utilised globally to describe stretching and folding. However, locally and instantaneously, stretching and folding are described by the spectrum of the (ever changing) Jacobi matrix and the recognition of this underpins the present method. Spatial discretisation. The choice of spatial discretisation method is not central to the method, and the approach taken here may well be substituted with another method depending on application. Here, all fields on the BM are discretised using finite element theory with Hermite shape functions: where ϕ i (m) are the shape functions, K is the number of degrees of freedom used to discretise the BM and, e.g., are nodal values of the position of the BM and N is the dimension of the system. Time derivatives and spatial derivatives for the BM are easily obtained from these expressions. The use of Hermite (higher order) shape functions allows for the definition of curvature, which appears in Eq (6). All simulations use the same 40-node (K = 80) mesh of 1-D elements with isoparametric Hermite shape functions. Numerical integration over m is performed with a full Gauss integration scheme with 7 Gauss points per element (K G = 280). Nodal updates are determined by solving for nodal rates and integrating numerically over t with an explicit time integration scheme. The forward Euler method is used for ease of implementation. Introducing the approximations in Eq (18), multiplying Eq (3) by ϕ j and integrating over m on both sides, we obtain I Which can be written in matrix form and solved as where Now, the time integration scheme can be utilised to update c. The system matrices from the finite element formulation do not change between time steps, which can be exploited for efficiency.
For a higher order attractor, the expression for the right hand side is modified according to the attractor dimension M as The inward velocity d 1 can be evaluated at each point on c directly using Eq (6). The tangential velocity d 2 is found by the finite element solution to Eq (9): where Auxiliary manifold orientation -concepts. There exists an extensive body of literature on how to define and practically determine invariant manifolds. A comprehensive overview is given by Gorbana and co-workers in [53]. This is an active area of research and new highly advanced methods are continually being developed [54]. Here, an approach is taken, which is somewhat related to what is labelled the "method of invariant grids" in [53]. It differs substantially in that the tangent field of the invariant manifold is not found by numerical differentiation in phase space. This approach, though intuitive, requires refined meshes and is prone to numerical instability. Instead, a hybrid approach is taken, where AM orientation is determined by a spectral method, but selection of eigenvectors from this spectrum, and thus orientation, makes use of the grid. The reasoning behind this choice is pragmatic and based on experience: Even if very large time steps are taken between two orientation updates, the tangential motion of the AM in between these updates will be small when compared to typical mesh refinement. Furthermore, the use of a spectral method renders the orientation of the inward normal u 1 independent of interior points. The obvious drawback is the need for solution of the eigenvalue problem at each Gauss point and at every time step. A redeeming feature is that this can be done efficiently by use of sparse iterative solvers that solve for only the M necessary eigenvalues and utilise eigenvalues from previous time steps as initial guesses. The current software implementation does not exploit these features and computes the entire spectrum from full matrices using direct methods.
Locally, the tangent space of the AM is spanned by the tangent vectors w i 2 R N where i = 1, 2, . . ., M. Chaotic motion is often described as successive processes of stretching and folding. As discussed, locally, these processes correspond to exponential growth and decay, respectively. The selection of the tangent vectors w i is such that every direction in which stretch occurs is contained within the AM. In this way, every off-manifold motion will be an exponential decay onto the manifold.
Consider the Jacobi matrix of Eq (2) In any non-linear system, the Jacobi matrix depends on the state, i.e., J = J(x). In non-linear non-autonomous systems it may also depend explicitly on time, i.e. J = J(x, t). The method is equally applicable in both of these situations. The AM is oriented such that all directions of stretch are contained within it at all times. Locally and instantaneously, these directions are those spanned by eigenvectors associated with eigenvalues with a positive real part. That is, the manifold tangent vectors w i are selected such that they span the same space as solutions to the local system associated with directions of stretch. The number of eigenvalues with positive real parts determines the required dimension of the AM. For two positive real parts at least a surface is required, and for three a 3-D hypersurface, and so on. Knowing the dimension of the attractor a priori is advisable, since then a good (if not perfect) indication of the dimension of the AM is then known. A good account of the estimation of attractor dimension, for both experimental data and systems with known ODEs, is given by Wolf and co-workers [55].
When a direction of stretch with associated eigenvalue λ is encountered somewhere in phase space, there are two situations where v 1 and v 2 are the eigenvectors associated with eigenvalues λ and " l.
2. l 2 R: The direction of stretch is parallel with the eigenvector associated with λ.
In case 1, one may define the local tangent vectors in any way that span the same space as v 1 e lt þ v 2 e " lt , e.g., by selecting two distinct phases: 0 (giving v 1 + v 2 ) and +/−i (giving i v 1 − i v 2 ). If the flow is real v 1 and v 2 appear as complex conjugate, so the space is spanned by Re(v 1 ) and Im(v 1 ). These vectors are then Gram-Schmidt orthonormalised to yield the orthonormal basis (w 1 ,w 2 ) of the AM. With the requirement of smoothness, the orientation of the AM is defined even in regions of phase space with no directions of stretch. This implies the requirement that J(x) is continuous and differentiable. Even if a direction of stretch is associated with a positive real eigenvalue, there may be another region of phase space where this eigenvalue becomes complex (with positive or negative real part) and thus requires a second direction of expansion. The only way to know this is to trace the eigenvalue along the AM and track its evolution. A schematic illustrating this process is shown in Fig 11. The phenomenon of directions of stretch varying in number across phase space is well known and is what is labelled "unstable dimension variability" in, e.g., [39].
Auxiliary manifold orientation-initialisation. Selecting the appropriate eigenvectors during initialisation takes a two step approach: First, the eigenvalues of all Gauss points on the BM are determined until a node is reached where at least one eigenvalue has a positive real part. Subsequently, the rest of the Gauss points are scanned incrementally along the BM selecting the ones corresponding with minimum increments in modulus. If the discretisation of the BM and AM is sufficiently fine, the minimum increment selection criterion will correspond to continuity. Thus, eigenvalues for neighbouring points on the AM and the BM can be directly organised and numbered consistently. Although it is somewhat computationally wasteful, it is conceptually more straight forward to perform the sorting process before the selection of eigenvalues. This is illustrated as pseudo-code in algorithm 1. Now, the number of entries in l j equalling 1, should exceed the a-priori predicted dimension of the attractor. If it does not, then the initialisation of the BM is inadequate and another must be tried.
The above pseudocode is written to communicate the concept of the method and is not optimised. For example, there is no need to sort eigenvalues that do not correspond to directions of stretch.
Auxiliary manifold orientation-updates. Updating l ðpÞ j at each time step simply follows the principle of continuity in time, so updates in eigenvalues are those minimizing changes in modulus from the previous time step. At each time step, the inward normal u 1 is projected onto the orthonormal basis associated with the corresponding eigenvectors. Thus, at each time step, Gauss point and interior point the Jacobi matrix is updated, the eigenvalues of the Jacobi matrix are calculated, and the eigenvectors associated with minimum change in eigenvalue modulus are used to update the orientation of the manifold. If Δt λ is sufficiently small, this corresponds to continuity in time of manifold orientation. The process of time integration and manifold updates is outlined as pseudo-code in algorithm 3. Note that the time step between orientation updates is larger than the time step for position updates. For clarity, the algorithm is written for a 2-D AM (M = 2) with l ðpÞ 1 and l ðpÞ 2 corresponding to the eigenvalues of orientation at gauss point m p . Generalisation to higher dimension is straight forward. J = J(c(m), t i ) Solve J v j = λ j v j for λ j and v j , j = 1, 2, ..N Select λ j from the spectrum of eigenvalues such that jl ðpÞ 1 À l j j is minimised l ðpÞ 1 l j Select λ k from the spectrum of eigenvalues such that jl ðpÞ At the first time step, updates in orientation are associated with large changes in u 1 if the initial guess is poor. Conversely, if the initial guess for the AM is close to the real, long term, AM, the projection will not alter u 1 much. In this paper, a flat, circular surface is used as an initial guess and initial values of w 1 and w 2 computed on this surface. This approach works well when off-manifold directions are almost in phase with on-manifold location, but may fail if off-manifold dynamics are out of phase, potentially causing the AM to become degenerate and numerically unstable. Thus, convergence depends on a reasonable initial guess. This which might potentially be generated from Poincaré sections in future incarnations of the method. Fig 3. shows a manifold changing its shape in the initial transient phase of a simulation. Fig 12 depicts the converged manifold and the manifold eigenvalues for parameter set III in Table 1. The manifold is projected to the (y 1 , y 2 , y 81 ) space. Where complex, as expected, the two eigenvalues of orientation are complex conjugate. There is a central region of stretch where λ 1 is positive. Conceptually, this corresponds to the area around the unstable saddle equilibrium at (0, 0) in the unmodified 2-D Duffing system.

Discussion
A practical method for calculating extrema of chaotic systems has been presented here. The equations of the method are written such that they can be wrapped around an underlying system of which the detailed structure is less important, as long as its flow is continuous and differentiable. The simulation results are post-transient extrema of systems with fairly fast, close to in-phase, off-manifold dynamics. Such systems pose the least challenging in terms of defining (guessing) an initial shape, that will rapidly converge onto the real AM. Addressing this challenge, for example by use of fits to Poincaré sections, would increase the chances of convergence.
The presented method has drawbacks which should be addressed in future work: Firstly, it will not distinguish between basins of attraction as the BM is defined a priori as a single closed curve. For this, an approach based on the level set method would be superior but perhaps at the cost of added computational work. Secondly, the presence of noise would also make the BM itself noisy. The method in its current incarnation offers no strategies to deal with this added challenge, which will inevitably come with the application to experimental data. Thirdly, the adaptation to experimental data would entail a phase space reconstruction step. The method, in its current incarnation, requires knowledge of the rates of all phase space variables, so will need to be modified if phase space reconstruction techniques based on time-delay coordinates [56,57] are to be employed. However, it may well be that a hybrid scheme where some rates are inferred from knowledge of the rest would suffice.
The computational benefits associated with defining the BM as low-D as possible are potentially enormous. Consider the cost of operating on a conventional N − 1 BM rather than an M − 1 ditto. The memory needed to store the BM, let alone operate on it, scales as K D , where D  (12). The parameters (set III) are given in Table 1.
https://doi.org/10.1371/journal.pone.0179507.g012 is BM dimension and K is the resolution of the BM, e.g., as in this case, the number of nodes of a finite element formulation. Plugging in the numbers for even modest sized systems quickly reveals why the D = N − 1 approach is infeasible. For example, storing an N − 1 BM, for an N = 82 phase space system with identical resolution to the one used in the presented analyses requires 1.13 Á 10 155 bytes if stored as double precision.

Conclusion
Based on the presented results, it is concluded that it is practically possible to predict extrema for large differentiable chaotic systems and that these extrema have the potential to display harmonic behaviour. This enables the efficient calculation of extrema of chaotic systems in the long term once the post-transient regime is reached.