Invariant spectral foliations with applications to model order reduction and synthesis

The paper introduces a technique that decomposes the dynamics a nonlinear system about an equilibrium into low order components, which then can be used to reconstruct the full dynamics. This is a nonlinear analogue of linear modal analysis. The dynamics is decomposed using Invariant Spectral Foliations (ISF), which is defined as the smoothest invariant foliation about an equilibrium and hence unique under general conditions. The conjugate dynamics of an ISF can be used as a reduced order model. An ISF can be obtained from vibration data without carrying out model identification first. It is also shown how arbitrary forcing can be taken into account by constructing a parametrisation, so that forcing only minimally breaks the invariance of the ISF.


Introduction
In this paper we highlight how invariant foliations [12,27] of dynamical systems can be used to derive reduced order models (ROM) either from data or physical models. We consider dynamics about equilibria only. We assume a deterministic process, such that future states of the system are fully determined by initial conditions. An invariant foliation is a decomposition of the state space into a family of manifolds, called leaves, such that the dynamics brings each leaf into another (see figure 1). If a leaf is brought into itself, then it is also an invariant manifold. A foliation is generally characterised by its co-dimension, which equals the number of parameters needed to describe the family of leaves so that it covers the state space. The dynamics that maps one leaf of an invariant foliation into another leaf has the same same dimensionality as the co-dimension of the foliation. We call this mapping the conjugate dynamics, which is lower dimensional than the dynamics of the underlying system and therefore suitable to be used as an ROM. Such ROM treats all initial conditions within one leaf equivalent to each other and characterises the dynamics of the whole system. In contrast, the conjugate dynamics (ROM) on an invariant manifold captures the dynamics only on a low-dimensional subset of the state space. The conjugate dynamics on an invariant manifold, however describes the exact evolution of initial conditions taken from the invariant manifold, while the conjugate dynamics on an invariant foliation is imprecise about the evolution, it can only tell which leaves a trajectory goes through. This ambiguity about the state has some advantages: for all initial conditions there is a leaf and a valid reduced dynamics. In contrast, when using invariant manifolds, the initial condition must come from the invariant manifold in order to have a valid prediction.
Multiple foliations can act as a coordinate system about the equilibrium. When individual leaves from different foliations intersect in one point, the dynamics can be fully reconstructed from the foliations. Therefore invariant foliations are fully paralleled with linear modal analysis of mechanical systems [9]: it allows both the decomposition of the system and the reconstruction of the full dynamics. To reconstruct the dynamics one needs to find intersection points of leaves from different foliations which is more complicated than adding vibration modes of linear system. However, such composability is not at all possible with invariant manifolds or any other nonlinear normal mode (NNM) definition [14,21,10,20]. Therefore, an invariant foliation seems to be the closest nonlinear alternative of linear modal analysis. The concept of composition is illustrated in figure 1.
Invariant foliations can be directly fitted to time-series data, because the foliation acts as a projection, much like linear modes. This allows for another parallel to be drawn with modal testing [9], which identifies linear vibration modes from data. For invariant manifolds, such fitting is not available. In [24], a two-step process was used to find invariant manifolds in vibration data. First an underlying high-dimensional model was identified and then the invariant manifold was extracted. Direct fitting of the manifold invariance equation to data would require a projection that maps from the state space to the parameter space of the manifold, which is not a-priori available and cannot be easily calculated. In fact, the calculation of such a projection is equivalent to finding an invariant foliation. Moreover, a leaf of a foliation that is mapped into itself or equivalently in our case contains the equilibrium, is an invariant manifold. Therefore finding two complementary invariant foliations, one transversal to an invariant manifold, another containing the invariant manifold as a leaf can substitute for calculating the invariant manifold and ROM.
The condition for uniqueness of invariant foliations are different from invariant manifolds. Only invariant manifolds about equilibria that are sufficiently smooth are unique. Unique invariant manifolds about equilibria, periodic or quasi-periodic orbits are called spectral submanifolds (SSM) [10]. The Figure 1 Two foliations act as a coordinate system. An initial condition (red dots) is mapped forward by F , however each leaf of a foliation is brought forward by the lower dimensional maps S 1 and S 2 . Due to invariance of the foliation, the full trajectory can be reconstructed from the two maps S 1 and S 2 and the leaves of the foliations. theory behind SSMs was mainly developed in [8], generalised to infinite dimensions in [6] and applied to mechanical systems in [10]. For an SSM to exist, non-resonance conditions need to be satisfied and the dynamics must be smoother than a so-called spectral coefficient, which is calculated from the eigenvalues of the Jacobian about an equilibrium. For an SSM to be interesting, it must contain the slowest dynamics, so that it captures long-term behaviour, rather than just transients (see R2 in [11]). It turns out that the spectral coefficient of such an SSM is also the highest and therefore the SSM requires the highest order of smoothness to be unique. While the concept of smoothness is theoretically well-understood, it is almost impossible to quantify numerically or determine from data. This is one of the reasons why it is challenging to calculate SSMs numerically (in contrast to series expansion [17]) in a reproducible manner. Invariant foliations, as explained below, also need to satisfy non-resonance conditions to exist and sufficiently smooth to be unique. We call a unique invariant foliation tangential to an invariant linear subspace about an equilibrium an invariant spectral foliation (ISF). In contrast to SSMs, ISFs that capture the long term dynamics require the lowest order of smoothness among all ISFs. This eases the difficulty of numerical calculation and allows to put a practically acceptable error bound on numerical approximations (for an ISF with spectral coefficient of one, the errors grows linearly with the distance from the equilibrium).
The existing literature on invariant foliations is rich and difficult to summarise without distracting too much from the purpose of the paper (see e.g. [12,2,19]). However, the setting used here is also different from most of the literature in that we are not dealing with stable or unstable fibres and hyperbolicity is not an important aspect either. The closest results in the literature are the remarks of de la Llave in section 7.3 of [8] and section 2 of [6], that generalise the parametrisation method to foliations.
The paper also introduces a technique to account for forcing in reduced order models derived using SSMs or ISFs. Previous studies [25,5,18] carried out a similar analysis, but they assume that the type of forcing (periodic, quasiperiodic) is known before model order reduction is carried out. The major advantage of knowing the type of forcing ahead of model order reduction is that forcing can be part of the SSM calculation as in [18]. When the type of forcing is unknown, such calculation cannot be carried out. Here we do not assume the type of forcing. The main idea is that the reduced order model should capture as much of the forcing as possible through a re-parametrisation of the ISF or SSM. We achieve this by constraining the parametrisation of the SSM such that the invariant manifold does not move in the tangential direction of itself or that the leaves of an ISF are only perturbed in their tangential directions by the external force. In appendix C we find out that the requirements for the SSMs and ISFs lead to equivalent reduced order models.
The plan of the paper is as follows. We start with introducing invariant foliations and describe their properties. We then state and prove theorems for the existence and uniqueness of ISFs both for discrete-time systems and vector fields. Subsequently, we describe a simple method that allows finding ISFs from time-series data, which is then tested on a simple example alongside with two other approaches. Finally, we describe our technique to account for arbitrary external forcing in lower order models derived using SSMs and ISFs.

Invariant foliations
Consider a dynamical system that is defined by the C r map F : R n → R n . A trajectory of the dynamical system is then defined by recursively applying F to the initial condition x 0 , such that successive points along a trajectory are generated by We assume that the origin is a fixed point, that is F (0) = 0 and the Jacobian at the origin, A = DF (0), is semisimple. The eigenvalues of A are denoted by µ i , i = 1, . . . , n and we have a full set of left and right eigenvectors, v i and v i , that satisfy v i A = µ i v i and Av i = µ i v i , respectively. Let us denote the linear subspace spanned by the first ν eigenvectors as E = span {v 1 , . . . , v ν } and the dual subspace E = span {v 1 , . . . , v ν }.
We are interested in how codimension-ν sets about the origin are brought into each other by F . The manifold of sets is parametrised by an ν-dimensional parameter z ∈ R ν and a single set at point z is denoted by L z . We assume that each L z is a differentiable manifold and L z and Lz are disjointed if z =z. In technical terms this is called a codimension-ν foliation of R n [15] and each L z is a leaf. The foliation is a collection of leaves, that is F = {L z : z ∈ R ν }. A foliation F is invariant under F if there is a map S : R ν → R ν , which brings the leaves into each other in the same way as the high-dimensional dynamics, that is A foliation can be represented by a function U : R n → R ν , called submersion, such that a leaf is the pre-image of the parameter z under the submersion U , that is, Using definition (3), we find that the inclusion (2) translates into an algebraic equation for the submersion U , which is called the invariance equation. Similar to invariant manifolds, we require a tangency condition to a linear subspace. To consider the dynamics corresponding to the linear subspace E , we require that which means that DU (0) is a set of ν linearly independent row vectors from the dual space of R n (row vectors) spanning the whole space E . Figure 2 shows the geometry of an invariant foliation. Each leaf is mapped into another, in particular, the green solid line representing L z is mapped into the red solid line by F . A leaf that corresponds to a fixed point of S is an invariant manifold as it is mapped into itself. In particular, we have S (0) = 0, hence L 0 is an invariant manifold.
The solution of the invariance equation (4) with the tangency condition (5) is not unique for a number of reasons. Firstly, assuming that there exist a pair of functions U and S satisfying (4) and (5) a large class of diffeomorphism Φ : R ν → R ν can be used, such thatŨ = Φ • U andS = Φ • U • Φ −1 are also solutions of (4) and (5). However if two pairs of solutions of (4) and (5) are conjugate through a diffeomorphism Φ, they represent the same invariant foliation F. The kind of non-uniqueness that is problematic when multiple solutions of (4) and (5) are not conjugate and do not represent the same foliation. To fix possible non-uniqueness, we impose extra smoothness conditions on the submersion U in addition to being differentiable and tangent to E . We then call the smoothest and unique foliation the invariant spectral foliation (ISF) corresponding to the linear subspace E.

Vector fields
In many applications the dynamics is defined by a vector field. Here we recall that there is a one-to-one relationship between invariant foliations of maps and vector fields [1]. Consider the vector fieldẋ = G (x), which has a fundamental solution Φ t (x), such that Here Φ t is a one-parameter group, because Φ t (Φ s (x)) = Φ t+s (x) and Φ 0 (x) = x. We can now define the map F (x) = Φ t (x), which brings the invariance equation (4) into The conjugate dynamics S must also be a one-parameter group with S t+s (x) = S t (S s (x)) and S 0 (x) = x in order to satisfy the invariance equation, that is, The infinitesimal generator of the group S is denoted by R, such that d dt S t (x) = R (S t (x)). On the other hand U must be independent of time, if it is to define an invariant foliation. We now take the derivative of the invariance equation (6) with respect to time and find Setting t = 0 in equation (7), we get the invariance equation for vector fields in the form of The next example, which aims to illustrate non-uniqueness of foliations, also shows that occasionally, it is easier to find an invariant foliation using (8) than using (4).

Example: smoothness and uniqueness of foliations
Let us consider the discrete-time map x n+1 y n+1 = e −λ x n e −µ y n , λ > 0, µ > 0 for which we can find an equivalent vector field in the form of such that x n = x (n) and y n = y (n). The solutions of system (9) lie on the curves y (x) = ce xµ/λ , c ∈ R, x ≥ 0, as we only consider the right half-plane. The invariance equation (8), when (9) is substituted, becomes where r describes the dynamics among the leaves of the invariant foliation.
Here, we have used non-bold, lower-case letters to represent U = u and R = r, because they assume scalar values. Without restricting generality, we prescribe the parametrisation of the foliation by setting u (x, 0) = x, which implies that r (x) = −λx. We note that any other parametrisation for whicĥ u (x, 0) is a strictly monotonous (invertible) and smooth function of x can be brought into the special parametrisation that we have just chosen, that is u (x, y) =û û −1 (x, 0) , y . Using this parametrisation, the invariance equation then simplifies to The solution of (10) is sought in the form of u (x, y) = xw (x, y), where w has to satisfy the somewhat simpler equation Using the method of characteristics and assuming the boundary condition w (1, y) = f (y) gives the general solution where f is an unknown, continuously differentiable function with f (0) = 1 due to the constraint on the parametrisation. We now assume that f is m times differentiable, such that f (x) = 1 + m k=1 a k x k +O x m+1 and that λ/µ > m. In this case the k-th order term of f leads to an order 1+k (1 − µ/λ) > 1+k−k/m term in u, that are continuously differentiable if and only if k ≤ m. This implies that if m < λ/µ ≤ m + 1, function f must assume the form Figure 3 Uniqueness of foliations. Dashed blue lines are trajectories of (9), red continuous lines are the contours of u (x, y) and represent the leaves of the foliation. (a) f = 1 + x/5, λ = 2, µ = 3, the resulting u does not define a differentiable foliation; (b) f = 1 + x/5, λ = 3, µ = 2, the foliation is once differentiable but not unique; (c) f = 1, λ = 2, µ = 3 leads to the unique and differentiable foliation.
for u to be once differentiable. This means that the foliation is non-unique and has m free parameters. Repeating the same argument but stipulating that the foliation must be m-times continuously differentiable we find that f = 1, which has no parameters and therefore the invariant foliation becomes unique. Indeed, after differentiating (11) m-times, a k-th order term in f results in an hence none of the terms apart from the constant one will lead to an m-times differentiable u, hence the only solution is f = 1 and the unique submersion is u (x, y) = x.
We also note that in this example, the ISF is as smooth as the vector field, that is analytic. The x variable represents the slow dynamics if λ < µ. In this case λ/µ < 1, which means that a differentiable foliation is already unique. The result of this section is graphically illustrated in figure 3.

Existence and uniqueness of invariant foliations
In this section we generalise the findings from the example in section 2.2 and provide a sufficient condition for the existence of a unique invariant foliation, i.e., an ISF. We start with a definition.

Definition 1
The number E = min k=1...ν log |µ k | max k=1...n log |µ k | is called the ISF spectral coefficient of the linear subspace E about the origin.
Theorem 1 Assume that there exists an integer σ ≥ 2, such that E < σ ≤ r. Also assume that n k=1 µ m k k = µ j , j = 1, . . . , ν for all m k ≥ 0, 1 ≤ k ≤ n with at least one m l = 0, ν + 1 ≤ l ≤ n and with n k=0 m k ≤ σ−1. Then in a sufficiently small neighbourhood of the origin there exists an invariant foliation F tangent to the invariant linear subspace E of the C r map F . The foliation F is unique among the σ times differentiable foliations and it is also C r smooth.
Proof The proof is carried out in appendix A.

Remark 1
The value of the ISF spectral coefficient, when all eigenvalues µ k are within the complex unit circle must be 1 ≤ E . In this case theorem 1 implies that the invariant foliations must be at least twice differentiable to be unique. This is a stricter requirement than what we have found for example (9) in section 2.2, where simple differentiability already implied uniqueness.
The proof follows the same lines that Cabré et al. [6] employ. First, a low-order series expansion is carried out avoiding possible resonances. For higher order terms, where no resonance is possible, Banach's contraction mapping principle is applied to find a unique correction. Since the series expansion allows a number of free parameters, we also show that the choice of these parameters does not influence the geometry of the foliation, only its parametrisation. This results in a unique foliation. Differentiability directly follows from formally differentiating the invariance equation.
Theorem 1 also applies to C r vector fieldsẋ = G (x). Again, we assume that the origin is the equilibrium, that is G (0) = 0 and that the Jacobian B = DG (0) is semisimple. The eigenvalues of B are denoted by λ i , i = 1, . . . , n and we have a full set of left and right eigenvectors, v i and v i , that satisfy v i B = λ i v i and Bv i = λ i v i , respectively. The invariant linear subspaces are defined as before: Using the spectral mapping theorem for the equivalence A = exp Bτ , τ > 0, we find that the ISF spectral coefficient for a vector field is Due to the equivalence between discrete-time dynamics and vector fields, the following corollary is a direct consequence of theorem 1.
Corollary 1 Assume that there exists an integer σ ≥ 2, such that E < σ ≤ r. Also assume that n k=1 m k λ k = λ j , j = 1, . . . , ν for all m k ≥ 0, 1 ≤ k ≤ n with at least one m l = 0, ν + 1 ≤ l ≤ n and with n k=0 m k ≤ σ − 1. Then in a sufficiently small neighbourhood of the origin there exists an invariant foliation F tangent to the invariant linear subspace E of the C r vector field G. The foliation F is unique among the σ times differentiable foliations and it is also C r smooth.

Fitting an invariant foliation to data
In this section we outline how to find the immersion U and conjugate map S from a time-series without identifying map F . We assume that the time-series is already split into two linear invariant subspaces E 1 and E 2 . To find these linear subspaces one can use standard linear regression that fits a linear model to the data in the neighbourhood of the equilibrium [4]. Therefore we assume that the underlying map has the form where We do not assume the values of matrices A 1 and A 2 , however their initial guess from the linear estimate can speed up convergence. The time-series is denoted by (x k , y k ) , k = 1, . . . , N , which is produced by the map F , such that x k+1 , y k+1 = F (x k , y k ). The invariance equation must therefore be satisfied for each point along the time-series, that is where r k is a small residual, which is due to the noise in the data. The obvious strategy to finding an invariant foliation is to minimise the residual r k , which we do using the least-squares method. We use the norm (43) from the proof of theorem 1 in appendix A.2, which guarantees a unique solution in the proof of theorem 1. Therefore our version of least-squares optimisation problem can be written as where we require that U (x, 0) = x. This might not be the best parametrisation (which is problem dependent), for example it cannot account for cases when the foliation is folded over the x coordinates. However, the simplicity of this parametrisation leads to a straightforward implementation.

Polynomial representation for optimisation
Here we use a polynomial representation to carry out the optimisation given by equation (15). We represent the unknown functions as polynomials of finite order α, such that where the powers are interpreted as with m ∈ N ν or n ∈ N n−ν , respectively. We further introduce the notation |m| = m, where the summation applies to all elements of m, depending on the dimensionality of m. The multi-index notation implies that coefficients of linear terms have indices given by unit vectors e k = 0 Matrices are consequently denoted as multi-indexed vectors, that is, the element of a matrix in the j-th row and k-th column is written as S e k j or just simply the k-th column vector of a matrix is written as S e k . We also prescribe the parametrisation of the foliation by stipulating that U (x, 0) = x, which becomes U e k ,0 j = δ jk and U m,0 = 0 ∀ |m| = 1 using the polynomial representation, where δ jk is the Kronecker delta. The polynomial representation of the objective function in the optimisation problem (15) can be written as Note that the composition S • U is not evaluated to the highest order in (17), instead the resulting polynomial is truncated at order α so that U and S • U can be fully matched. This results in a consistent representation. We have also observed more accurate results in our example in section 4.3 (data not shown) as opposed to not limiting the polynomial order of the composition S•U . From experience with other model identification studies, we believe that accuracy can be improved if not just two consecutive points, but multiple points along a trajectory are taken into account. This leads to a so-called 'multiple shooting' technique [3], which will be part of a further investigation. In our implementation we use the Optim.jl [16] package of the Julia programming language and choose the BFGS method to find an optimal solution. This only requires the calculation of the gradient of f , which is simple to calculate.

Reconstructing the dynamics
Two or more carefully selected ISFs can act as a nonlinear coordinate system of the state space and therefore can be used to reconstruct the dynamics of F .
Let E 1 and E 2 be invariant linear subspaces, both satisfying the conditions of theorem 1, such that E 1 ∩ E 2 = {0} and E 1 ⊕ E 2 = R n , the corresponding submersions of the ISFs are U 1 and U 2 , respectively. Let us now assume trajectories (x k , y k ), z 1,k and z 2,k , k ∈ N satisfying x k+1 , y k+1 = F (x k , y k ), z 1,k+1 = S 1 (z 1,k ), and z 2,k+1 = S 2 (z 2,k ) with matching initial conditions, that is, Because of the invariance of ISFs, these trajectories satisfy the equivalence equations Due to our assumptions about E 1 and E 2 , we can invert U in a neighbourhood of the origin and therefore there exist two functions h 1 and h 2 , such that Under h 1 and h 2 , the equivalence is expressed as Equation (19) can be used to reconstruct the full dynamics of the system from the lower order conjugate dynamics of the ISFs. The equivalence is the same between the trajectories of the vector fields G, R 1 and R 2 , except that the subscript k is replaced by time t.
The functions h 1 and h 2 can be obtained by a fixed point iteration in a small neighbourhood of the origin. Due to the constraints U 1 (x, 0) = x, U 2 (0, y) = y, the submersions U 1 and U 2 can be written as Therefore, the inverse of U can be obtained by the convergent iteration For polynomials of a given order the iteration converges in finite number of steps.
In some circumstances, for example when the SSM is sought as a leaf of a foliation or visualising the foliation, it is helpful to represent the leaves of a foliation as graphs over some coordinate axes. Here we assume U (x, 0) = x and look for a function g that satisfies where z is the parameter of the foliation and each leaf is then given by k 2 k 1 m 1 1 Figure 4 Schematic of a single mass, two degrees of freedom (DOF) oscillator. Damping is not represented in the diagram as it is assumed to be linear.
Note that g (0, z) = z due to the parametrisation of U and that L 0 is an SSM. Given the polynomial representation and our choice of the parametrisation, we Therefore we can apply the fixed point iteration to find g. Given that U N is a polynomial of order α, and we seek g as another polynomial of order α, starting with g 1 (y, z) = z, the iteration (22) finishes in α steps.

Numerical example: single mass, two DOF oscillator
Consider a point mass m in figure 4, that is mounted to two springs of stiffnesses k 1 and k 2 , respectively. The point mass is allowed to move in plane.
The mechanism therefore has two degrees of freedom. This system was used to test various techniques that calculate nonlinear normal modes and SSMs [25,26,5], hence we do not repeat the derivation of the model. For simplicity, we assume linear damping, although dampers applied parallel with the springs can also be included at the cost of longer equations of motion. The equations of motion in a first order form reaḋ 2 and x 1 is the horizontal displacement, x 2 is the vertical displacement, x 3 is the horizontal speed and x 4 is the vertical speed of the point mass. We use the parameters m = 1 kg, k 1 = 1 N/m, k 2 = 3 N/m, c 1 = 1/100 Ns/m and Ns/m. The forcing P is set to zero in this section, but it will assume a time-dependent value in section 5.3. The natural frequencies are ω 1 = k 1 /m = 1 1/s and ω 2 = k 2 /m = √ 3 1/s and the corresponding ISF spectral coefficients are This also implies that it takes longer for vibrations with the first natural frequency to decay and therefore the ISF and SSM corresponding to E 1 contain the slowest dynamics. We carry out three kinds of ISF calculations, one directly for the vector field as given in appendix B, another for the discrete time model identified from vibration data as in appendix A.1 and yet another when the ISF is identified from the vibration data directly as outlined in section 4.1. To create fitting data for the last two ISF calculations, we solve equation (23) numerically. We create 10 trajectories of length t = 128 s, and use a T = 0.8 s time-step to sample the solution, which leads to 160 samples per trajectory. The initial conditions are chosen from a normally distributed sample with variance radius of 0.2 (m, m/s) about the origin. The ISF decomposition is always carried out for the two linear subspaces E 1 and E 2 . The reduced order models are then solved either by finding the trajectories of the two-dimensional vector fields R 1 and R 2 or iterating the initial condition using the maps S 1 and S 2 .
In figure 5(a) we have used the submersions, as in equation (18), to map the full solution into the coordinates of the ISF components and compared that with the solutions of the conjugate dynamics given by R 1 , R 2 and S 1 , The error measure is given by where the subscript fw refers to forward prediction. In figure 5(b) we have used equation (19) to fully reconstruct the dynamics from the two ISFs. For the same trajectories as in figure 5(a), the error is now calculated as where the subscript bw refers to backward prediction. When vector fields are used the indices k are replaced with the independent time variable in formulae (24) and (25). It can be seen that decomposing the vector field directly yields similar accuracy to the ISFs of the discrete time map. The directly fitted ISFs however behave differently. They display roughly the same error regardless of the order of polynomials and they are most accurate for the first few timesteps, after which the error significantly grows. This is due to the identification method, which minimises the error for the first time-step, but disregards further evolution. As mentioned, this can be fixed by using a multiple shooting technique [3]. ISFs also contain the SSM as a leaf going through the origin. Figure, 6 shows this SSM as the leaf L 0 of the foliation corresponding to the second natural frequency (green surface), as given by equation (21), which is the SSM for the first natural frequency. We have also tracked a leaf from the foliation corresponding to the first natural frequency (black lines) for four iterations of the conjugate dynamics and compared this to how initial conditions on this leaf are mapped forward by the solution of the vector field (23) (red dots). This was carried out for ISFs identified from the vector fields (23) directly and also ISFs identified from simulation data by optimisation. It can be seen that the SSM is visually the same for the two methods, however the leaves look different. Nevertheless the invariance of the leaves for these 4 steps does hold, because the forward trajectories denoted by the red dots are mapped onto consecutive leaves.

Applying forcing to a reduced order model
So far we were concerned with ISFs that describe autonomous motion, however in many applications, which includes mechanical systems, external forcing is an important part of the problem. In this section we show that forcing can be included into the framework of SSMs and ISFs in an equivalent manner, although the insights gained from the two approaches are different. To set up the discussion, we first define the class of systems we are concerned with, summarise some preliminary results and state our result. In this section we restrict our attention to vector fields.  Figure 6 State space projections:

Preliminaries and statement of the result
Let us assume the differential equatioṅ where G : R n+1 → R n is a C r vector field and u represents a scalar-valued time-dependent external forcing. In this section we assume that there exists a u-dependent family of ISFs with submersion U : R n+1 → R ν and conjugate vector field S : R ν+1 → R ν that satisfy the invariance equation and the tangency condition Similarly, we assume that there exists a family of SSMs, with immersion W : R ν+1 → R n and conjugate dynamics R : R n+1 → R ν that satisfy and the tangency condition The SSM is then the manifold given by We note that the notation has somewhat changed from the previous sections as we needed to introduce SSMs. However, in this section we only consider vector fields and therefore symbols previously used for maps, such as R and S are reused. To illuminate the connection between an ISF and a SSM corresponding to the linear subspace E we also define Υ (y, u) = U (W (y, u) , u) , which is invertible in the first variable in a neighbourhood of the origin, because of the tangency conditions (28) and (30). The inverse of Υ in the first variable is denoted by Υ −1 that satisfies the equation Υ Υ −1 (y, u) , u = y. We will use this notation to represent the inverse of functions of type R ν+1 → R ν in the first variable in this section and in appendix C. The function Υ connects the possibly different parameterisations of the ISF and SSM.

Lemma 1
The conjugate vector fields S and R of an ISF and a SSM corresponding to the invariant linear subspaces E and E , respectively represent the same dynamics, that is, Proof Substitute the immersion W into equation (27), then use equation (29) in the resulting equation to replace G • W with D 1 W R. Then notice that For completeness, we include the basic theorem for existence and uniqueness of SSMs. Here we are only concerned with a single SSM, which is denoted by M. Theorem 2 Assume that there exists an integer σ ≥ 2, such that ℵ E < σ ≤ r. Also assume that ν k=1 m k λ k = λ j , j = ν + 1, . . . , n for all m k ≥ 0 such that ν k=0 m k ≤ σ − 1. Then in a sufficiently small neighbourhood of the origin there exists an invariant manifold M tangent to the invariant linear subspace E of the C r vector field G. The manifold M is unique among the σ times differentiable manifolds and it is also C r smooth.
Proof The theorem is a subset of theorem 1.1 in [6]. Practical implications are discussed in [10,24].
We now state the result of this section. Let us define V (y, u) = D 1 U (W (y, u) , u) , which is a rectangular matrix valued function. The kernel of matrix V is the invariant normal bundle of the SSM or in different words the rows of V span the adjoint tangent space of the SSM. This can be seen at the origin, where by definition (28) we have spanV (0, u) = E . Due to invariance, the adjoint tangency continues to be the case over the whole SSM. Note that V can be calculated directly and more efficiently than U as described in section C.5.
We now want to ensure that a trajectory that is on the SSM, will not be perturbed in the tangential direction of the SSM, because that perturbation can be fully captured by the conjugate dynamics R. Therefore any perturbation due to forcing must fall into the invariant normal bundle of the SSM and therefore be in the kernel of V . This is the same setting as in [23, section 4.3], which is extended here from the perspective of ISFs. The perturbation of the invariance equation (29) by an external forcing is D 2 W (y, u)u, therefore we require that V (y, u) D 2 W (y, u) = 0.
Our goal is to find a parametrisation of the SSM, such that (33) holds. This can be achieved using the following proposition. An exhaustive explanation of why (33) is a good choice from various perspectives is included in appendix C.
Proposition 1 Let Ψ : R ν+1 → R ν be the unique solution of the initial value problem In what follows, equation (36) will be used as a reduced order model capturing the external forcing. We will illustrate for two examples that near the resolved natural frequencies and for harmonic forcing, the solution of (36) is nearly indistinguishable from the full system. The first example also illustrates what happens when forcing drives the system near the radius of convergence of the series expanded ROM.

A modified Shaw-Pierre example
Here we investigate the Shaw-Pierre example [21] modified by Ponsioen, Pedergnana and Haller [18, section 5.1] using our method, where the same example was used to study isolated forced responses (isolas). The differential equations areẋ The SSM of the first natural frequency ω 1 (or invariant subspace E 1 ) is unique among the five times differentiable invariant manifolds as indicated by the SSM spectral coefficient ℵ E1 . The ISF is already unique among the twice differentiable invariant foliations. The comparison between the forced response of the full model and the reduced order model can be seen in figure 7. We find that the reduced order model is able to reproduce the vibration even around an isola. To find the forced response, numerical continuation was used [22], which is able to compute unstable periodic orbits that occur at higher forcing amplitudes. Note that the instability of the forced response about the second vibration mode with natural frequency ω 2 cannot be predicted by the ROM. This is because the instability occurs in the normal direction of the SSM, which is ignored by the conjugate vector field R. The departure of the response of the order-3 ROM seems to be unrelated to this instability.

One-mass oscillator revisited
Here we look at the forced response of equations (23) of the schematic in figure  4. The forcing occurs along the horizontal direction and therefore the second  vibration mode corresponding to the vertical motion does not contribute significantly to the overall response. Here, we only investigate the first vibration mode.
Using the polynomial expansion techniques that we have derived in this section, their forced response of the order 3 and 5 ROMs were compared to the series expanded equation (23) of the same order. The result can be seen in figure 8. The results show good agreement between the full model and the ROM for low forcing amplitudes, but display discrepancy for higher forcing amplitudes. In this case the discrepancy is due to the small radius of convergence of the series expansion.

Conclusions
The paper has introduced invariant spectral foliations (ISF) as a tool to derive reduced order models (ROM) of dynamic systems about an equilibrium. We have carried this out both for forced and autonomous systems. The major advantage of this approach over other methods is that the full dynamics can be reconstructed from a set of ISFs. We have used polynomial expansion about equilibria to calculate the ISFs and ROMs and demonstrated some of the properties of ISFs. Polynomial expansion can be restrictive because it may not converge for the part of the phase space that we would like to resolve, which we have seen in examples. Therefore it would be beneficial to try other kinds of representations, for example the universal function approximators used by machine learning algorithms, such as deep neural networks [13]. We have demonstrated that ISFs can also be directly learned from vibration data, however the accuracy of the ROM deteriorated for long-time predictions. We believe that this can be remedied by using a multiple shooting type optimisation technique, so that longer trajectories are taken into account when learning the ISF. One particularly desirable feature of an ISF is that its representation only needs to be once differentiable when it is calculated for the slowest dynamics (as in section 2.2). This is in contrast with SSMs, where the SSM containing the slowest dynamics must be many times differentiable (as given by the SSM spectral coefficient) to be unique. The required order of smoothness gets higher if there is a time-scale separation with an increasingly fast dynamics in the system. As smoothness is difficult to quantify numerically calculating unique SSMs can be a challenge. Therefore ISFs are better suited for systems with time-scale separation when capturing the slow dynamics is the goal. This aspect will be elaborated elsewhere.

A.1 Polynomial expansion
In this section we find an approximate solution to the invariance equation (4) in the form of a power series. Here we employ a complexification [1] of the vector space R n which is isomorphic to C n and therefore the complexified map F may not be defined on the whole of C n , because F is not an entire function. We now apply a linear transformation T to the map F , such that its Jacobian about the equilibrium becomes a diagonal matrix and the map assumes the form where, If there is a complex conjugate pair of eigenvalues µ k = µ k+1 the corresponding components of the state variable must also be complex conjugate x k = x k+1 . Similarly, if µ k is real, x k must also be real.
To arrive at an approximate solution of the invariance equation (4), we represent the unknowns as power series in the form of where the powers are interpreted as with m ∈ N n or m ∈ N ν , respectively. The notation used here is explained in section 4.1.

The equation for the linear terms in the invariance equation becomes
which does not have a unique solution, however we choose U e k j = δ jk , S e k j = µ k δ jk . This is a normalising constraint that will be taken into account when we prove the uniqueness of the foliation F . The equations for the order |m| terms in the invariance equation are written as where H m j are the terms which are composed of lower order terms of U and S and known |m|-th order terms of F . The equations are written for two different kinds of exponents. Equation (39) is for exponents that exist for both U and S and therefore part of the conjugate dynamics, equation (40) is for exponents that only identify terms in U and they correspond to dynamics that occurs inside the leaves. Equations (39) and (40) are solved recursively starting with |m| = 2 and then in increasing order for |m| > 2.
Equation (39) can be solved under any circumstances, but there are multiple solutions. The terms U m j and S m j can be chosen relative to each other. If there is a resonance or near resonance, that is ν k=1 µ m k k ≈ µ j we can choose the solution otherwise we can also choose or some other combination of U m j , S m j . The choice made here is another normalising condition and as we see it will not affect the uniqueness of the foliation. Equation (40) has a unique solution if n k=1 µ otherwise no solution exists, unless H m j vanishes, which is unlikely. It turns out that if |m| > E , no further resonances are possible. We then choose the smallest σ ∈ N + , such that E < σ and denote the truncated series expanded solution of the invariance equation by U ≤ and S ≤ .

A.2 Contraction mapping
Here we show that once an order σ − 1 asymptotic solution of the invariance equation is found, then there is a unique C σ correction of the asymptotic solution U ≤ and S ≤ , so that the invariance equation (4) is exactly satisfied. For the following argument we re-scale the map F , such that F γ (x) = γ −1 F γ (γx), where γ > 0. We have now solved the invariance equation up to order σ − 1 and therefore the invariance equation (4) has an order σ residual when the approximation is substituted.
Let us now fix S = S ≤ and decompose the exact solution into U = U ≤ + U > , where U > is a C σ function. To obtain an iterative solution, we apply the inverse S −1 to the invariance equation (4) and find that Equation (41) can be re-cast as a fixed point iteration using a nonlinear operator. This nonlinear operator is defined as We define operator T on the space X σ = C σ (B n , R ν ), where B n = {x ∈ R n : |x| ≤ 1} is the closed unit ball. The norm on X σ is chosen to be which makes X σ a Banach space. Next, we show that operator T is a contraction. In particular, we show that there exists a constant L < 1, such that and that T (U ) σ ≤ 1 for all U σ ≤ 1. The latter criteria can be demonstrated through the Lipschitz condition (44) using the estimate which means that we need to have T (0) σ < 1 − L for T to be a contraction. We start by estimating T (0) σ . Due to the polynomial approximation, we have which implies that there exists M > 0 such that Using the scaled nonlinear map F γ and the scaled approximate solutions, the estimate now scales with γ as This result implies that T (0) ≤ γ σ−1 M. Since σ ≥ 2, the norm T (0) can be made arbitrarily small and therefore any Lipschitz constant 0 ≤ L < 1 is sufficient to demonstrate a unique fixed point of T . Next we need to show that there is a Lipschitz constant 0 ≤ L < 1 in equation (44). We use the fundamental theorem of calculus to make a calculation similar to For our operator T , we write that Due to the scaling, S becomes linear as γ → 0, hence we can find 1 (γ) such that lim γ→0 1 (γ) = 0 and This implies the estimate In the next step, we are estimating the effect of the inner function F γ of U > 1,2 by way of the σ-norm, that is Similar to the previous estimates, there exists 2 (γ) such that lim γ→0 2 (γ) = 0 and sup |x|≤1 |F γ (x)| = max k=1...n |µ k | + 2 (γ) .
Putting together the previous estimates we find that (45) Since 1 → 0 and 2 → 0 as γ → 0, we only need to show that max k=1...n so that the Lipschitz constant in equation (45) is less than one. After rearranging the criterion (46) for σ, we find that if operator T is a contraction, which is certainly true because we assumed σ > E .

A.3 Uniqueness
In section A.1 we have made some normalising assumption, which restricted the parametrisation of the foliation. Here we show that those assumptions can be removed and that any sufficiently smooth submersion U satisfying the invariance equation (4) can be transformed so that it satisfies the normalising assumptions and yet represents the same invariant foliation, which then implies uniqueness. Let us write variable x as a tuple x = (x 1 , x 2 ) such that x 1 ∈ C ν and x 2 ∈ C n−ν with the restrictions due to complexification. When finding a solution of (4) as per the argument in sections A.1 and A.2, we can use normalising conditions such that the solution of (39) satisfies U m j = 0 when |m| ≥ 2. This then implies that U (x 1 , 0) = x 1 + O (|x 1 | σ ). Let us now denote the unique solution of (4) under these normalising conditions by U and S and another σ-times differentiable solution of (4) and (5) byÛ andŜ. We now define 0) is an invertible and σ-times differentiable function in a sufficiently small neighbourhood of the origin, because its Jacobian at the origin is invertible due to the tangency condition (5). Consequently, we define Φ (z) = U Ψ −1 (z) , 0 , which is again an invertible σ-times differentiable transformation. Let us now definẽ which satisfy the invariance equation (4), the tangency condition (5) and our normalising conditions. However we have shown that there is a unique solution to (4) and (5) under the normalising conditions, henceŨ = U andS = S. This means that solutions of (4) and (5) can be re-parametrised into each other, therefore the invariant foliation F represented by the submersions U orÛ is unique in a sufficiently small neighbourhood of the origin.
This conclude the proof of theorem 1.

B Series expansion of ISFs for vector fields
In section A.1 we have provided an algorithm to find a power-series expansion of an ISF for a discrete map. Here we modify the algorithm for vector fields. We assume a first order differential equationẋ = G (x), whose vector field is transformed into the frame of the eigenvectors of its Jacobian about the origin, such that where, To arrive at an approximate solution of the invariance equation (4), we represent the unknowns as power series in the form of where the powers are interpreted as with m ∈ N n or m ∈ N ν , respectively. Using the notation in section 4.1, the equation for the linear terms in the invariance equation becomes where H m j are the terms which are composed of lower order terms of U and S and known |m|-th order terms of F . The equations are written for two different kinds of exponents. Equation (48) is for exponents that exist for both U and S and therefore part of the conjugate dynamics, equation (49) is for exponents that only identify terms in U and they correspond to dynamics that occurs insides the leaves. Equations are solved recursively starting with |m| = 2 and then in increasing order for |m| > 2.
Equation (48) can be solved under any circumstances, but it is clear that there are multiple solutions. The terms U m j and S m j can be chosen relative to each other. If there is a resonance or near resonance, that is ν k=1 m k λ k ≈ λ j we can choose the solution otherwise we can also choose or some other combination of U m j , S m j . Equation (49) has a unique solution if n k=1 m k λ k = λ j , which is otherwise no solution exists, unless H m j vanishes.
C Derivation of proposition 1 from the ISF and SSM perspectives

C.1 Invariant foliation perspective
Here we discuss how the invariance of an ISF can be approximately preserved if external forcing is applied. Let us assume a time-dependent u and equate the time-derivative of y = U (x, u) with the conjugate dynamicsẏ = S (y, u) to arrive at the equation Equation (50) becomes the invariance equation (27) if D 2 U (x, u) = 0, which is not possible for all x, because U would then be independent of u. However, it is only necessary to minimise D 2 U (x, u) about the SSM, so that invariance is minimally affected by the forcing. The simplest strategy is to require that Due to continuity of the derivative D 2 U , the loss of invariance will be small in a small neighbourhood of the SSM if (51) holds.
In order to ensure that (51) holds, we find a new parametrisation of the submersion U . We introduce the diffeomorphism Φ : R ν+1 → R ν with conditions Φ (0, u) = 0 and Φ (y, 0) = y, such that the new submersion and conjugate vector field that satisfy the invariance equation (27) becomẽ respectively. Taking the derivative ofŨ with respect to u and substituting W (y, u) yields We notice that U (W (y, u) , u) = Υ (y, u), which is invertible, hence the differential equation for Φ becomes Equation (54) can be solved by the method of characteristics. The reduced order model that fully captures forcing on the SSM iṡ z =S (z, u) .
However it might be convenient to not carry out the transformation, but equivalently introduce an additional term to the conjugate dynamics of the ISF. Let us now differentiate z = Φ (y, u) with respect to time, which leads tȯ Substituting z = Φ (y, u) into (56) and multiplying it by [D 1 Φ (y, u)] −1 from the left yields an alternative reduced order model in the original frame aṡ Equation (57) has some advantages over (55), because it does not require the inverse of Φ, only the inverse of its Jacobian. The form of S may also be simpler thanS, as it can be reduced to the set of terms that correspond to (near) internal resonances.

C.2 Invariant manifold perspective
In this section we investigate the deviation of solutions of the externally forced (26) from the SSM up to linear order and thereby prove proposition 1. We then restrict the parametrisation of the SSM such that forcing does not perturb the autonomous solution in the tangential direction of the SSM up to a linear order. This is the same setting as in [23]. The technical difficulty is that one requires a projection that projects the difference between the perturbed and the autonomous solution onto the tangential direction of the SSM. Since invariant tangential and normal directions about the SSM are not orthogonal, the calculation of the projection is not trivial. We assume that the solution of system (26) is represented by the dynamics on the SSM and a perturbation z, such that We now derive the equation that perturbation z must satisfy up to linear order. Taking the time derivative of (58), we geṫ Then taking into account the invariance equation (29) and that G (W (y, u) + z, u) = G (W (y, u) , u) + D 1 G (W (y, u) , u) z + O |z| 2 , the linearised equation for z becomeṡ However, we want to avoid that z has components in the tangential direction of the manifold, because that can be captured by the conjugate vector field R, if the immersion W is appropriately parametrised.
To capture the tangential component of the perturbation z, we now define a family of projections that project from R n to the tangent space of the SSM. We require that the projection is invariant under the linear dynamics about the SSM and therefore it can also be thought of as the normal vector bundle of the SSM [27]. Assume a linear projection V : R ν+1 → L (R n , R ν ), where L (R n , R ν ) stands for the space of linear operators from R n to R ν (ν ×n matrices). In order to ascertain that this projection fully captures the tangential direction of the SSM, we require that V (0, u) D 1 W (0, u) = I. (For the most general result, identity is not necessary, an invertible V (0, u) D 1 W (0, u) would suffice). The projection is invariant if there exists a linear (time-dependent) vector field in the form ofξ = C (y, u) ξ that describes the evolution of ξ = V (y, u) z if z is the solution ofż = D 1 G (W (y, u) , u) z. Taking the time derivative of ξ we find thaṫ Substituting ξ = V (y, u) z,ẏ = R (y, u) and simplifying by z, the differential equation for the projection V becomes We require that D 2 W does not have a tangential component as that can be captured by the conjugate vector field R, which occurs when Again, we introduce a re-parametrisation by a diffeomorphism Ψ , so that the new quantities becomeW (y, u) = W (Ψ (y, u) , u) , which is an ordinary differential equation for Ψ (y, u) in the independent variable u with initial condition Ψ (y, 0) = y. The ROM that captures all the forcing in the tangential direction of the SSM becomesẏ =R (y, u). As before in section C.1, instead of usingR, an extra term added to the conjugate dynamics R can capture the forcing. The alternative version of the ROM can be written aṡ which requires the inverse of Ψ , but preserves the structure of (near) internal resonances captured by R.

Remark 3
We note that projection V is directly related to the invariance equation (27) of the ISF, that is V (y, u) = D 1 U (W (y, u) , u) and C (y, u) = D 1 S (Υ (y, u) , u). This result is obtained by differentiating (27) with respect to x and substituting x = W (y, u).

C.3 Equivalence of the ISF and SSM perspectives
Here we show that ROMs derived from the ISF and the SSM perspectives represent the same dynamics. First we start with a lemma.
Lemma 2 Let U be a submersion of an ISF and W an immersion of an SSM, tangential to the linear subspaces E and E, respectively. Let us define Υ (y, u) = U (W (y, u) , u). If the transformations Φ and Ψ satisfy any two out of the three equations (53), (62) and then the third equation also holds.
Proof To simplify the notation, composition of two functions is now denoted by •. First, Υ = U • W is substituted into (53), which gives Next, we substitute V = D 1 U • W into (62) and find and then using (66) results in Finally, substituting Ψ −1 into (68) gives (65), which is what we needed to prove. This calculation can be carried out in reverse and therefore any of the three equations can be reached from the other two.
ThereforeS andR represent the same dynamics even when u is time-dependent.
Proof We have the two transformed reduced order modelsS andR that are defined by and the relation between S and R as given by lemma 2 as Substituting Υ into (70) and using (71) yields Then substituting Ψ and using (69) yields Now denote Ξ (y) = Φ (Υ (Ψ (y, u) , u) , u), which is independent of u because equation (64) holds if the defining equations of Φ and Ψ , that is (53), (62) hold by lemma 2. Then notice that (72) can be written asS which is a conjugacy independent of parameter u, thereforeS andR represent the same dynamics even if u is time-dependent.

C.4 Notation for parameter-dependent calculations
In what follows we take a somewhat simplified approach so that previous methods on calculating autonomous SSMs and ISFs can be reused without modification even with the additional parameter u. Parameter u can be made part of the state space by attachinġ u = 0 to the unforced equations. We therefore define ξ = (x, u) T , the extended vector field and the new parameter of the ISF and SSM as η ∈ R ν+1 , η = (y, u) T . Such an extended equation is not covered by corollary 1 of the ISF or theorem 2 of the SSM, because neither E nor ℵ E are defined for Ge. However, series expansion still works, except that all terms will be u dependent and none of the powers of u can be excluded from the conjugate vector fields Se, Re, because they all lead to internal resonances. For practical purposes, such a series expansion can be helpful.

C.5 Calculating the invariant vector bundle
While it is possible to obtain the invariant vector bundle projection V from the submersion of the ISF and the immersion of the SSM, tangential to E and E, respectively, it is computationally more efficient to calculate the projection directly from equation (59). Using the concise notation of section C.4, and dropping the e subscript from various functions, the invariance equation (59) for the bundle projection becomes DV (η) R (η) + V (η) DG (W (η)) = C (η) V (η) .
Firstly, we assume a series expansion of the known quantities R (η) and DG (W (η)). The unknown quantities are then represented as where the powers are interpreted as η m = η m 1 1 · · · η m ν+1 n , with m ∈ N ν+1 . The quantities V m and C m can be represented as three-dimensional arrays and we denote the elements by V m pq , C m pq , respectively. When substituting η = 0 into (74) we get V 0 B = C 0 V 0 whose solution is C 0 = B 1 and V 0 = (I, 0) .
Note that B and B 1 are defined in (47). Terms of order |m| must satisfy where H m pq are composed of lower order terms of V and C and known |m|-th order terms of D 1 G and R. The equations (75) and (76) are written for two different set of indices. Equations (75) and (76) are solved recursively starting with |m| = 1 and then in increasing order for |m| > 1.
Equation (75) can be solved under any circumstances, but it is clear that there are multiple solutions. The terms V m pq and C m pq can be chosen relative to each other. If there is a resonance or near resonance, that is ν k=1 m k λ k ≈ λp − λq, we can choose the solution otherwise no solution exists, unless H m j vanishes, which is unlikely. The non-resonance conditions ν k=1 m k λ k = λp − λq p = 1, . . . , ν, q = ν + 1, . . . , n are related to the conditions (13). Let us definem such thatm k = m k for k = 1, . . . , ν and m k = δ kq for k = ν + 1, . . . , n, in which case we have n k=1m k λ k = λp.
We now note that (77) is the same condition as (13) except that a smaller set of exponents m are considered. Therefore if the conditions (13) hold, then (77) must also hold. This is not surprising, because according to remark 3, V can be calculated from the submersion of an ISF.

C.6 Calculating the re-parametrisation
While there are two equivalent ways to obtain a reduced order model that includes forcing, we choose to re-parametrise the SSM, because it is computationally more efficient. The SSM needs to be calculated even when the ISF approach is used, but only the projection V and transformation Ψ is necessary. All of these functions depend on low number of variables in contrast to the submersion of an ISF, which has the same number of parameters as the dimension of the state space. We now want to solve equation (79) for Ψ . Therefore, we account for various terms in equation (79) in terms of the notation introduced in section C.4. The solution of (74), according to remark 3, can be written as which is a square matrix. Formula (78) includes the term D 2 U (W (y, u) , u) that does not appear in equation (79) and we need to eliminate it. Therefore, we definẽ V e (y, u) = I 0 0 0 V e (y, u) + 0 0 0 1 and find thatṼ e (Ψ e (y, u)) DW e (Ψ e (y, u)) D 2 Ψ (y, u) = 0 1 is the same equation as (62), which can be used to find Ψ . Therefore we simply need to solve the ordinary differential equation D 2 Ψ e (y, u) = Ṽ e (Ψ e (y, u)) DW e (Ψ e (y, u)) −1 0 1 with independent variable u and parameter y. Using a polynomial representation, the most efficient solution of (79) can be carried out by Picard iteration [7] Ψ k+1 e (y, u) = y 0 + u 0 Ṽ e Ψ k e (y, τ ) DW e Ψ k e (y, τ ) which incorporates the initial condition Ψ e (y, 0) = (y, 0) T . The iteration (80) converges in finite number of steps if we limit the order of Ψ e to a finite value. The resulting ROM is then given by equation (57).
Funding: The author is supported by a University Research Fellowship awarded by the University of Bristol.