SIMULATION OF COMPLEX DYNAMICS USING POD ‘ON THE FLY’ AND RESIDUAL ESTIMATES

. Proper orthogonal decomposition (POD) is a very eﬀective means to identify dynamical information contained in sets of snapshots that cover numerically computed trajectories of dissipative systems of partial diﬀerential equations. Such information is organized in a hierarchy of POD modes and the system can be Galerkin-projected onto the associated linear subspace. Quite frequently, the outcome is a low dimensional model of the problem. Flexibility and eﬃciency of the approximation can be enhanced if POD is applied ‘on the ﬂy’, adaptively combining a standard numerical solver (which provides the necessary snapshots) with the reduced system in interspersed intervals, as both time and a bifurcation parameter are varied. Residual estimates are introduced to make this adaptation accurate and robust, preventing possible mode truncation instabilities in the presence of complex dynamics. All ideas are illustrated in some bifurcation scenarios including quasi-periodic and chaotic attractors, which highlights a good computational eﬃciency.


1.
Introduction. Analysis of transitions and bifurcations has stimulated a large amount of research works in the last decades, thanks to the wide range of scientific interests and fields the topic embraces [10]. Also, the increasing need to account for nonlinearities in industrial environments is making this area highly relevant in engineering applications [22].
Depending on the particular class of problems, there exist several families of methods to deal with bifurcation phenomena. For instance, specific bifurcations can be analyzed via projection onto the center manifold to derive the normal form equations [19] or numerical continuation [1] can be performed to more flexibly follow branches of steady and/or periodic solutions. More complex (quasi-periodic and chaotic) dynamics, instead, require a Poincaré map, whose construction relies on time dependent simulations, disregarding transients. Unfortunately, for realistic descriptions of many problems modeled by partial differential equations or systems, current computational tools may require non-affordable CPU time and memory to construct Poincaré maps involving complex attractors. Indeed, a numerical solver must be run for each value of the bifurcation parameter over a sufficiently large time span, to reach the asymptotic state. Thus, new methodologies are strongly demanded to support with computational efficiency the huge body of knowledge developed on nonlinear dynamics and bifurcations.
If the system is dissipative and the spatial domain is bounded, the large-time behavior is contained in a finite dimensional inertial manifold [12], which opens the possibility of low dimensional descriptions. Hence, reduced order models (ROMs) [21] may be feasible and flexible tools for reducing the number of involved degrees of freedom. In particular, ROMs based on proper orthogonal decomposition (POD) consist in low dimensional systems obtained upon projection of the governing equations onto 'customized' linear subspaces; these have been recently developed and applied to a variety of physical problems [2,3,8,23,11,9,20,27,29]. The POD modes are constructed from appropriate sets of snapshots, often calculated at the outset by means of a standard numerical solver. This offline preprocess may be a computational bottleneck but can be minimized by an adaptive, alternate integration (in interspersed time intervals) of the full and a reduced systems [23,29]. A major drawback of POD-based ROMs concerns the high-order mode truncation instabilities, which may promote deviations of the approximation from the actual dynamics in a somewhat unpredictable way [25]. So far, various explanations and solutions have been provided by the ROM community on the issue [7,9,26].
Bifurcation problems have also been tackled by means of ROMs [16,18,15]. In general, snapshots are calculated in a neighborhood of the bifurcation point to be analyzed and are of the same type as the relevant states in the involved attractors. In other words, the preprocess is performed 'ad hoc', depending on the structure of the bifurcation diagram, which strongly penalizes flexibility. As far as the authors are concerned, simulation via ROMs of quasi-periodic and chaotic attractors has received less attention. A quite flexible approach has been proposed in [30,31] for the general dissipative system of semilinear parabolic equations where q is the state vector, µ is the bifurcation parameter, L is a linear operator such that L −1 is well-defined and compact, and f is a nonlinear function such that the operator q → L −1 f (q, µ) is compact. These assumptions hold, in particular, for reaction-diffusionconvection systems, micro-fluidic systems, and pattern-forming thermal convection systems; they have been discussed in [30] and heuristically justify the feasibility of low dimensional (linear) descriptions. The expected continuous dependence of POD modes on time and the bifurcation parameter (also observed in related works; see [30] and references therein) suggests that POD subspaces can contain the dynamics of the system for values of t and µ not used for their computation, as long as some more POD modes than necessary are retained. In addition, a convenient set of snapshots can be calculated from transient orbits, which are generically expected to cover a large portion of the phase space. Thus, information extracted from such behaviors may reasonably yield rich POD bases. These ideas have been used in [31] as building blocks of a general, adaptive ROM to fast construct bifurcation diagrams of (1). The proposed method relies on a Poincaré map based on the solution of a reduced system obtained from some POD modes. These are constructed at the first value of µ and updated 'on the fly' by means of new snapshots (calculated out of short temporal transients of the full system) when some transition occurs and new dynamical features need to be taken into account in the ROM. Two appropriate error estimates reveal whether a local update of the POD basis is necessary: the first one guarantees the desired accuracy of the approximation, while the second one is designed to detect possible mode truncation instabilities.
The latter estimate, based on imposing consistency with the results provided by an instrumental ROM of a larger dimension, turns out to be effective but someway lacks flexibility and undesirably increases the effort of calculations. In this context, the goal of this paper is to improve the way truncation instabilities are monitored, in terms of both robustness and computational cost. To this end, the second instrumental ROM will be skipped and the mode truncation instabilities will be monitored using an estimate of a normalized residual of the discretized full system, which enhances the performance of the adaptive ROM.
The remaining of the paper is organized as follows. In §2 POD and Galerkin projection are briefly recalled in the context of bifurcation problems, and an adaptive ROM to fast compute bifurcation diagrams is introduced [31]. Then, §3 illustrates a residual-based strategy to control truncation instabilities that improves the previous work by the authors. Finally, some results on bifurcation scenarios for the complex Ginzburg-Landau equation are shown in §4, followed by concluding remarks in §5.

2.
Adaptive POD to compute time dependent bifurcations. Among the model reduction methods, POD is a widely studied technique already exploited in a variety of applications [17,21]. Thanks to redundancies existing in either numerical or experimental data associated with high dimensional systems (due to the underlying physical laws), this methodology is able to extract the most coherent features in terms of an optimal linear basis, thus providing low dimensional descriptions.
In this framework, a possible approach to evolution problems [28] starts up by numerically computing a particular trajectory of the system that is sampled in order to get some snapshots, namely a set of spatial portraits q 1 , . . . , q N of the state vector at N representative time instants. Then, the orthonormal POD modes are constructed as where α k i are the components of the eigenvectors of the (Hermitian, positive semi-definite) covariance matrix R ik = q i , q k , while the singular values σ 1 ≥ σ 2 ≥ . . . ≥ σ N ≥ 0 are the square roots of the corresponding eigenvalues. It is worth observing that, in order to avoid round-off errors in the above product of snapshots and enhance the accuracy of calculations, POD can be performed via singular value decomposition [31]. Now, POD-based reduction techniques aim at predicting the dynamical behavior of the system at future times by supposing that the state vector belongs to the linear subspace generated by the POD modes, namely q(x, t) n j=1 A j (t)Q j (x), where n is suitably selected according to the singular values spectrum. Thus, such expansion is substituted into (1), which is finally Galerkin-projected onto the POD basis to yield a reduced system of n ordinary differential equations for the time evolution of the unknown mode amplitudes A j (t). Computational efficiency of POD and Galerkin projection can be increased (as done throughout this paper) by using a nonstandard inner product based on few discretization mesh points [2,3,23,29] (similar ideas can also be found in [5]).
The (stable parts of the) bifurcation diagrams of (1) can be efficiently approximated by means of the adaptive ROM proposed in [31]. At the first value of µ in the desired parameter span, an initial condition that breaks all the symmetries of the system is set and a non-small transient of the corresponding trajectory is numerically computed and sampled to get a sufficiently large set of snapshots. In this way, a POD basis of nonsymmetric modes is generated, thus preventing a spurious restriction of the ROM approximation to symmetric branches. Then, a Poincaré map for the resulting reduced system is defined from the outward intersections of the orbits with the Poincaré hypersurface where q n ROM is the state vector reconstructed by the n most-energetic (i.e., with the largest singular values) POD modes. Note that the outward/inward intersections are associated with local maxima/minima of the function t → q n ROM 2 , which means that H contains all steady states, at least two points of each periodic orbit, and at least two points of any time oscillation of q n ROM for more complex attractors. Finally, the bifurcation diagram plots a quantity related to the computed states at such maxima (or minima) after disregarding all transients. At subsequent values of µ, the initial condition for the reduced system comes from the final attractor at the previous value of µ (slightly perturbed to break all symmetries); furthermore, according to remarks in §1, the same reduced system is maintained and integrated as long as the approximation is satisfactory. Upon forward continuation along the bifurcation diagram, new dynamics will likely appear but the POD modes computed at the outset might no longer be able to provide a good description. Accuracy of q n ROM within the error bound ε can be checked by imposing where the two expressions are equivalent thanks to the POD modes orthonormality. Thus, the integration of the reduced system for each value of µ is combined with the online calculation of E n1 n , which is a standard error estimate in spectral methods [13]. Note that its calculation requires the reduced system to retain n 1 > n POD modes. If control (4) fails, the POD subspace is updated by mixing the available modes with some new ones constructed out of new snapshots for the current value of µ. These snapshots are solutions of the full system, which is simulated over a small transient taking as initial condition the last computed state (slightly perturbed). In general, only few ones are enough, as the expected continuous dependence of the POD basis on t and µ makes the old and updated POD subspaces somewhat close to each other. In the updating process, all modes are suitably weighted as to be adapted to the local dynamics [31].

A residual-based control of truncation instabilities.
In the context of the method introduced so far, mode truncation instabilities (which may be enhanced by unexpected transitions as µ is varied) are monitored in the attractor by simultaneously integrating a higher dimensional 'instrumental' reduced system and testing consistency with the solution provided by the basic one [31]. Specifically, for each value of µ, integration of the basic ROM assisted by control (4) is performed over a suitable transient; then, both reduced systems are run at large time to additionally prevent truncation instabilities over the attractor, whose good approximation is the main concern of the strategy. Actually, the use of two ROMs can be awkward if the reduced system is hard to construct and entails an increase of the CPU effort.
Therefore, a different approach is proposed, which has already been used by us in a related context [24] and arises from the observation that truncation instabilities worsen the extent to which the full time-marching scheme is satisfied by the solution of the reduced system. Thus, in addition to (4), an estimate E n1 res of the normalized residual of the discretized system (1) is imposed to be smaller than a prescribed tolerance, namely where ε is the desired error bound and k is a tunable factor. A failure of condition (5) would then suggest a deviation of the ROM approximation from the actual dynamics. Note that only one reduced system is now integrated and mode truncation instabilities are monitored over the whole time span. On the other hand, the form of E n1 res depends on the used discretization but its calculation is quite computationally fast, as based on the values of the residual at a limited number of mesh points, in the spirit of the present adaptive ROM (see [31] and §2). In general, in order to get a good estimate, such grid points must be as many as few times the number of retained POD modes.

4.
Results. The improved strategy to detect truncation instabilities (see §3) is applied to the complex Ginzburg-Landau equation (CGLE) which is a widely studied nonlinear equation [4] that exhibits intrinsically complex (chaotic) dynamics, due to the modulational instability, if αβ < −1 and µ exceeds a threshold value. The CGLE (6) is invariant under transformations x → 1 − x, u → u e ic , with c ∈ R, and has attractors of the form u(x, t) = e iδt u 0 , with δ ∈ R. Thus, the associated bifurcation diagrams show a variable complexity depending on whether |u| = |u 0 | is (i) steady and spatially uniform, (ii) steady and spatially nonuniform, or (iii) time periodic (with a period incommensurable with 2π/δ). In order to get both the necessary snapshots and the reference ('exact') solutions, the full equation (6) is integrated by a standard Crank-Nicolson plus Adams-Bashforth scheme, based on time step ∆t = 10 −4 and spatial centered finite differences in a uniform grid of 1000 points. Furthermore, u(x, 0) = 0.1 (1 + i) cos(πx) + 3i cos(2πx) is taken as (nonsymmetric) initial condition and the accuracy required to the ROM is given by ε = 10 −2 . In the first example, the parameter values α = 5, β = −17 are set and the bifurcation diagram in the range 0 < µ ≤ 50 (which is discretized by ∆µ = 0.01) is illustrated in Figure 1. It shows the values of |u 0.25 | ≡ |u(0.25, t)| (left) and |u 0.75 | ≡ |u(0.75, t)| (right) at the outward intersections of the orbits with the Poincaré hypersurface (3) in the time interval (14,15], after neglecting all transients, as obtained by the full equation (top) and the adaptive ROM described in §2 (bottom). The approximation is fairly good and quite small (stable) branches are also correctly captured (see Figures 2 and 3). Some apparent discrepancies between top and bottom plots can be observed (e.g., in the vicinity of µ = 46), which however are just artifacts of symmetry-breaking and might appear due to the different errors of the two used methods. Indeed, these may follow two reflection-symmetric branches, hence a cross (left-right) comparison between top and bottom plots reveals a correct description in the problematic regions. The (n, n 1 ) = (7,9) POD modes computed at the beginning of the continuation process are updated 18 times in the range 0 < µ ≤ 50,  namely the POD subspace is either rotated or enlarged to account for new transitions in the dynamics. In particular, at µ = 20.18, 24.78, 48.95, 49.52, and 49.95 updating is performed after the residual-assisted control (5) reveals the occurrence of a truncation instability. And this is done with a little computational effort, since the evaluation of the residual requires terms already computed for the reduced equation integration. On the other hand, it is interesting to observe that instabilities are detected in regions (near µ ∼ 24 and µ ∼ 49) where the dynamics are highly complex. In these regions, several updating events occur when trajectories are not confined to the current low dimensional POD subspace, which lacks important directions of the attractor. After conveniently updating the POD manifold, at µ = 50, the numbers of POD modes are (n, n 1 ) = (24, 36).   In the second example, the parameter values α = −30, β = 30 yield an almost conservative case, which is a quite demanding test for a POD-based approach. Approximation of the bifurcation diagram in the range 0 < µ ≤ 50 (constructed as before) is satisfactory and can be appreciated in Figure 4. After few bifurcations, the dynamics become (over 36 < µ < 47) quasi-periodic and then chaotic. Figure 5 shows a comparison between some representative attractors computed by the full equation and the adaptive ROM. The latter starts up with (n, n 1 ) = (4, 6) POD modes, which are updated 4 times to adapt to the various transitions appearing in the considered parameter span. It is interesting to observe that the computed residual exceeds the imposed bound at µ = 12.9. A careful analysis reveals a progressive decrease of the characteristic timescale of the trajectories along the branches this value of µ belongs to. Then, when the dynamics turn fairly oscillatory, an enrichment of the POD subspace becomes necessary. The complex window over 36 < µ < 47 requires a higher dimensional description based on (n, n 1 ) = (10, 13) POD modes. A pointwise comparison between the involved (quasi-periodic and chaotic) asymptotic states is obviously meaningless, due to their intrinsic nature. But a closer look shows that the Poincaré map of the reduced system seems to provide an adequate distribution of values in this range, as evidenced in Figure 4 by the darker portions of the bifurcation diagram. In other words, the adaptive ROM seemingly preserves certain (spatio-temporal) topological features of the most complex attractors also. This can be appreciated in Figure 6, which illustrates the approximation of the attractor at µ = 38. Nevertheless, an additional quantitative (statistical) study should be done in order to assert such conclusions, which is currently under research.
Concerning the computational cost, the adaptive ROM turns out to be fairly efficient. Indeed, it requires less than 14 CPU hours in a 3.46 GHz processor to compute the bifurcation diagram in Figure 4, while the same calculations based on the full equation are performed in 3 CPU days. In other words, the acceleration factor is higher than 5. It is worth noting that the performance of the reduced modeling technique is remarkably improved by the residual-based control of instabilities, in comparison with a previous version of the adaptive ROM [31]. Indeed, according to what emphasized in §3, the speed-up is typically 25% higher than before. 5. Final remarks. This paper describes an adaptive ROM based on POD and Galerkin projection to efficiently compute bifurcation diagrams of some dissipative systems. The authors presented a previous version of the method in [31], which is now improved by a residual-assisted control of truncation instabilities that may play an important role in the low dimensional simulation of complex time dependent dynamics.
Residual estimates have been already used in the ROMs literature (see, e.g., [14,7]). The main difference of the current approach with other works stands in the use of the residual as (i) an indicator of the stability of the reduced system and (ii) a flexible, computationally inexpensive, new ingredient in the framework of the introduced adaptive ROM that enhances its performance. Indeed, the estimate in (5) is calculated by means of only a few mesh points, which are generally enough to predict possible deviations of the approximation from the actual dynamics.
The final method has been applied to the CGLE (6), showing a robust and accurate description of quite complex transitions. On the other hand, the adaptive ROM is expected to work fine also for other dissipative systems that mimic laminar flow instabilities. More complex systems, however, still remain beyond its means as generally require subtle combinations of various techniques [3,6].
In any event, the ideas underlying this paper are hopefully a step further in current efforts to increase the computational efficiency in the analysis of nonlinear systems, which are paramount in many scientific and industrial fields.