Uncertainty quantification of time-average quantities of chaotic systems using sensitivity-enhanced polynomial chaos expansion

We consider the effect of multiple stochastic parameters on the time-average quantities of chaotic systems. We employ the recently proposed \cite{Kantarakias_Papadakis_2023} sensitivity-enhanced generalized polynomial chaos expansion, se-gPC, to compute efficiently this effect. se-gPC is an extension of gPC expansion, enriched with the sensitivity of the time-averaged quantities with respect to the stochastic variables. To compute these sensitivities, the adjoint of the shadowing operator is derived in the frequency domain. Coupling the adjoint operator with gPC provides an efficient uncertainty quantification (UQ) algorithm which, in its simplest form, has computational cost that is independent of the number of random variables. The method is applied to the Kuramoto-Sivashinsky equation and is found to produce results that match very well with Monte-Carlo simulations. The efficiency of the proposed method significantly outperforms sparse-grid approaches, like Smolyak Quadrature. These properties make the method suitable for application to other dynamical systems with many stochastic parameters.


I. INTRODUCTION
The performance of real world systems is significantly affected by the uncertainty of the parameters that define these systems.Large research effort has focused on the quantification of the effect of such stochastic variations to a Quantity of Interest (QoI), usually a time-averaged quantity.This field of research is commonly known as Uncertainty Quantification (UQ) and efficient UQ methods have been developed for static and dynamic problems [2][3][4][5].
In the standard generalized polynomial chaos (gPC) method, originally proposed in [3,6], an orthonormal polynomial base that spans the stochastic space is used for the spectral representation of uncertain quantities.The spectral coefficients are computed with Galerkin projection, allowing for the efficient evaluation of the statistics of the QoI.However, the cost of gPC scales as ∼ m p , where m is the number of stochastic parameters and p the polynomial order of the expansion.This exponential growth is known as the 'curse of dimensionality' [4].Various approaches have been proposed to mitigate the rapid growth of the computational cost, such as sparse grid approaches like Smolyak grids [7], or adaptive methods that build a sparse polynomial chaos expansion (PCE) basis using least-angle regression [8].gPC methods have been successful in predicting the statistics of the QoI in many applications, such as fluid dynamics, mechanics, space, medicine, see for example [3,5,[9][10][11].Applications of gPC to chaotic systems have also appeared in the literature [12][13][14][15][16].
Another method for the computation of the spectral coefficients is based on the least-squares approach [4,17].To reduce the computational coast, the method is usually coupled with efficient multidimensional sampling techniques; for an overview of the different sampling algorithms, see [18].Recently, new sampling approaches that incorporate real-world data into the computation of the gPC coefficients have been introduced [19][20][21].However the cost still grows exponentially with the number of uncertain parameters, making such methods difficult to apply in systems with a large m.
In this paper we use the sensitivity-enhanced gPC, or se-gPC, a least squares approach for the computation of the spectral coefficients that is augmented with the sensitivity of the QoI with respect to the uncertain parameters, see [1] for a detailed description of the method and a review of previous works in this area.When all the sensitivities are computed efficiently in a single step with the adjoint method, the computational cost of se-gPC is reduced by a factor m, i.e. it scales as ∼ m p−1 , and the method becomes increasingly useful as the number of stochastic inputs grows.Efficient sampling algorithms can by employed to reduce the number of required evaluations, see [1] for comparison between two such algorithms.For the special case of first order spectral representation, i.e. for p = 1, the spectral coefficients can be estimated with a single direct and a single adjoint evaluation, regardless of the number of stochastic inputs.
The se-gPC is efficient because the sensitivities of the QoI with respect to all stochastic inputs at one sampling point can be estimated using a single adjoint evaluation.Adjoint methods have long been successfully used to estimate derivatives (sensitivities) for stationary or nonchaotic systems in aerodynamics [22], structural optimization [23], chemical kinetic systems [24] among many others.However, when the underlying system is chaotic, a small variation in an input parameter causes large deviation in the trajectory of the system in phase space (with respect to the reference trajectory); this is popularly known as 'butterfly effect'.Under these circumstances, standard sensitivity analysis tools (such as adjoint) fail to produce physically meaningful results [25,26].Math-ematically, the deviation between the two trajectories is due to the presence of one or more positive Lyapunov exponents.To address this problem, the least-squares shadowing method (LSS) and its variants were proposed in a series of papers [27][28][29].This method is relies on the shadowing lemma [30,31] and provides a systematic and rigorous approach for the computation of sensitivities of time-average quantities of chaotic systems.
Assume a reference trajectory of a dynamical system evaluated for a parameter value s, u(t; s).Shadowing methods aim to compute another trajectory at s + δs, u(τ ; s + δs) that 'shadows', or stays close to, the reference trajectory for the time frame T of the analysis.The sensitivity problem is reformulated as a minimization problem between the reference and the shadowing trajectories.This formulation regularizes the problem, the two trajectories remain close to each other, and can be used to compute accurate sensitivities for a long T .Note that the two trajectories start from different initial conditions, but for ergodic system this does not affect the time-average QoI.Also, the time variable in the shadowing trajectory is not be the same as t, thus a different symbol, τ , is used.
Various approaches have been proposed to improve the computational efficiency of the original LSS method.The Multiple Shooting Shadowing (MSS) algorithm [32] was introduced to reduce the memory requirements of the standard LSS.When coupled with a matrix-free pre-conditioner to improve the convergence rate [33], MSS has lower computational cost and memory requirements than standard LSS.In [34], the non-intrusive LSS (NILSS) approach was derived and applied successfully to the 2D flow over a backward facing step and later to 3D flows inside a channel and around a cylinder [35,36].The computational cost of NILSS methods scales linearly with the number of positive Lyapunov exponents (PLEs).Theoretical predictions indicate that the highest Lyapunov exponent scales with the inverse of the Kolmogorov time scale [37,38].Yet another approach is the formulation of the shadowing problem in the frequency domain; to this end, the shadowing harmonic operator was introduced recently [39].The cost of the method is case-dependent, but for the Kuramoto-Sivashinsky system, sensitivities are computed at a cost roughly equal to that of the baseline solution.
In this paper, we derive the adjoint of the shadowing harmonic operator.The adjoint formulation allows us to compute the sensitivity of a time-average QoI of a chaotic dynamical system with respect to multiple parameters at a cost independent of the number of parameters.These sensitivities are then used to compute the UQ spectral coefficients in the context of se-gPC.As mentioned earlier, for spectral order p = 1, the cost of se-gPC is independent of m.This is the first application of UQ in chaotic systems with computational cost independent of the number of stochastic parameters.The method is applied and tested to the stochastically forced Kuramoto-Sivashinski equation and the results are compared against reference data obtained from Monte Carlo simulations.
The rest of this paper is structured as follows: in section II an overview of the standard gPC and se-gPC formulations is given.The adjoint shadowing harmonic operator for a general dynamical system is derived in section III.The se-gPC is applied to the Kuramoto-Sivashinsky equation with a single and multiple uncertain forcing parameters in sections IV and V respectively.Finally, in section VI the main findings of the paper are summarized.

II. SENSITIVITY-ENHANCED UNCERTAINTY QUANTIFICATION
Consider a dynamical system governed by a set of ordinary differential equations, where u(t; s) ∈ R Nu is the vector of state variables and s ∈ R Ns is a set of control parameters that define the dynamics of the system (for example Reynolds number in the case of incompressible fluid flows).We assume that the vector field f : R Nu × R Ns → R Nu varies smoothly with u and s.In most practical applications, we are interested in a time-averaged quantity J(s) : R Ns → R, which is usually referred to as the Quantity of Interest (QoI), for example lift or drag coefficient of an aerofoil.We assume that the control parameters s are functions of m independent stochastic variables ξ i that form the vector ξ = [ξ 1 , . . ., ξ m ].Each random variable ξ i is characterized by a probability density function (PDF), w i (ξ i ) in the domain E i .We seek to estimate the effect of ξ to the statistics of J(s).This effect can be quantified via the generalized polynomial chaos (gPC) expansion.In gPC a complete probability space P = (Ω, Σ, dP) is defined, where Ω refers to the set of random events and the probability measure dP is characterized by the σ-algebra Σ.The vector ξ follows the PDF W = m i=1 w i (ξ i ), defined in the domain E = m i=1 E i .This stochastic space is spanned by a polynomial basis Ψ = {Ψ 0 , Ψ 1 , . . .} which is orthogonal to W with respect to the inner product, The polynomial basis is normalized so that ⟨Ψ j , Ψ j ⟩ = 1.When m > 1, Ψ is defined by the tensor product of the unitary polynomials ψ (i) , as in Ψ := ⊗ m i=1 ψ (i) = {Ψ 0 , Ψ 1 , . . .}.The base is truncated to a finite number of polynomials by limiting the order of Ψ j to p.In that case, the QoI J is written in spectral form as where the number of basis functions is given by, and ϵ(ξ) is the truncation error (due to finite P ).In UQ with gPC, the goal is to compute the spectral coefficients c (i) in a computationally efficient manner.The moments of J(ξ) can be easily computed algebraically from c (i) .
In this paper, the coefficients are computed via a Weighted Least Squares (WLS) approach.To this end, q realizations of ξ are defined, with the i-th realization written as . The QoI J is computed for q realizations and stored in the vector Q = [J(ξ (1) ), . . ., J(ξ (q) )] ⊤ .Defining the vector c = [c (0) , . . ., c (P ) ] ⊤ ∈ R P +1 , equation ( 4) can be written in matrix form as, where ψ is the measurement matrix with elements ψ ij = Ψ j (ξ (i) ), i.e. the i-th row contains the values of the orthogonal polynomial basis Ψ j evaluated at the i-th sample point ξ (i) , and ϵ is the vector of truncation errors.The spectral coefficients are computed by solving the following weighted least squares minimization problem, is a diagonal positive-definite matrix, to be defined later.The solution of ( 7) results in the normal set of equations, For system (8) to be well conditioned q ≫ P + 1. Evaluating the QoI J(ξ (i) ) at the q sample points is computationally expensive, and dominates the cost of the method.
As mentioned in the introduction, when the number of uncertain parameters m is large, the number of spectral coefficients P grows exponentially (equation ( 5) indicates P + 1 ∼ m p ), leading to large computational cost, known as the 'curse of dimensionality'.The problem can be mitigated by enriching system (6) with gradient information.This method, called the sensitivity enhanced gPC, or se-gPC, is presented in [1].Differentiating eq. ( 4) with respect to the k-th random variable at the j-th sample point we get, ∂J ∂ξ For each random variable k, the block of q equations (9) can be written in matrix form as, where matrix ∂ψ ∂ξ k contains the gradients of the basis functions, . . .∂Ψ P (ξ (1) ) and vector ∂Q ∂ξ k stores the gradient of J with respect to the k-th random parameter at the q sample points, By stacking together Q and ∂Q ∂ξ k , we define the following block column vector Similarly, by stacking together the measurement matrix ψ and its sensitivity ∂ψ ∂ξ , we define the following block matrix ϕ ∈ R (1+m)q×(P +1) We can therefore write the following (1+m)×q equations for the spectral coefficients c, As before, to compute the coefficients c, we solve the following minimization problem where W ′ is a block diagonal weighting matrix, consisting of 1 + m blocks W .The solution ĉ is obtained via the normal equations, Notice the similarity between equations ( 8) and (17).The weights W 1 2 are computed with asymptotic sampling, a version of coherence sampling, see [40].In this paper for simplicity (and without loss of generality) only Gaussian inputs are considered, and the weights can be computed analytically as, For more details on the weights calculation and for extension to other input distributions, see [40].
To avoid large values of q, and thus keep the computational cost low, it is important to sample the QoI effectively.Different algorithms to sample the stochastic space are presented in [18].In this paper, we apply QR decomposition.This is a greedy algorithm that maximises the determinant of a matrix; in this sense it is a D-optimal design method, see [18,41,42].The process is as follows: A large pool of q random sample points is generated (from the prescribed probability density functions) and the measurement matrix ψ is formed.The question is how to select a subset of at least P + 1 points from this large pool.To this end, we multiply equation ( 6) with the row selection matrix P ∈ R (P +1)×q , thus we have, At each row of P all elements are 0, except the element at the column that corresponds to the selected sampling point, which takes the value of 1.In D-Optimal experiment design, P is found as a solution to the following maximization problem, where det() denotes the determinant of a matrix.The solution to this problem is given via the pivoted QR decomposition, The index matrix P is chosen so that the diagonal elements r ii of R are ranked in descending magnitude Since each sample point offers 1+m equations, at least 1+m samples with the highest r ii scores are retained.Thus sensitivity enhancement reduces the computational cost compared to standard gPC by a factor m at the cost an adjoint evaluation at each sample point.We could have applied the same approach to system (15) directly.However, applying QR decomposition to matrix W 1 2 ψ ⊤ and computing at these points the QoI and its gradient is preferable compared to QR decomposition of the weighted augmented matrix W ′ 1 2 ϕ ⊤ , see [1] for a comparison between the two approaches.

III. ADJOINT SENSITIVITY ANALYSIS OF CHAOTIC SYSTEMS
In this section we present a method for the computation of sensitivities of time-average quantities of chaotic systems to multiple parameters.To this end, we derive the adjoint version of the shadowing harmonic operator introduced in [39].
The goal of the LSS is to find a shadowing trajectory at s + ds that stays in close proximity (i.e.shadows) the reference trajectory at s.If the underlying system is uniformly hyperbolic, the shadowing trajectory is guaranteed to exist.To achieve this goal, LSS solves the following minimization problem, see [32], where v(t) = du(τ (t);s) ds is the sensitivity of the solution u(t; s) to a change δs of s, η(t) = d ds dτ dt is the timedilation, while (22c) denotes the orthogonality between the vectors f (u; s) and v(u; s) at each point along the trajectory.The gradient dJ ds is then given by the following expression, The solution of ( 22) has shown to produce accurate sensitivities [27][28][29]35].
In [39] the shadowing operator was formulated in the frequency domain.The key idea is to replace the minimization (22a) with the periodicity condition, yielding the following set of equations, The new contribution of the present paper with respect to [39] is that an adjoint approach is taken, since sensitivities with respect to a large number of parameters are required.To this end, a Lagrangian function is defined, where R v ∈ R Nu and R η ∈ R are the residuals of (24a) and (24b), while λ ∈ R Nu and µ ∈ R are the adjoint state variables.This is expanded as, and using integration by parts, We seek to make the Lagrangian independent of v(t) and η(t).This is achieved by solving the following field adjoint equations, The gradients can be computed from, Notice that (26a) is similar to (22b); the difference is that the transpose of the Jacobean is used and the time dilation term ηf is replaced by the adjoint term µf .The adjoint normality constraint (26b) has as forcing the residual J − J.Note also that the periodicity condition (24c) for v extends also to λ, see (26c).
The above equations can be formulated in the frequency domain by expanding λ(t) and µ(t) in Fourier series as, where λk , μk are the Fourier coefficients and ω 0 = 2π T is the fundamental frequency.The index k ∈ [−ℓ, ℓ] denotes the harmonics with frequencies ω k = kω 0 .Introducing expansions (28) into (26a) and (26b), the system that yields the Fourier coefficients is written in compact form as where We define the block diagonal matrix, which is a block Toeplitz matrix, because each diagonal has the same block.Using T (T m ), system (29) can be written in matrix form as, where D = diag[D −q , . . ., D 0 , . . ., D q ] is a block diagonal matrix with and Λ = Λ −q , . . ., Λ 0 , . . ., Λ q
Matrix H = D − T (T m ) is also known as a Hill matrix.Defining the adjoint shadowing harmonic operator as, the solution of system (32) can be written symbolically as The adjoint operator A maps the forcing R to the vector Λ that contains the unknown Fourier coefficients of the adjoint variables.More details can be found in [39,43,44].We do not directly compute the adjoint shadowing operator A. Instead we apply LU decomposition to the Hill matrix H = D − T (T m ) and find Λ by forward and back substitution.The gradients can then be computed from, which is the equivalent of ( 27) in Fourier space.The cost of solving system (32) is independent of the number of parameters and thus the adjoint formulation can provide the sensitivities of time-average quantities of chaotic systems is a single step.This information is used to augment the se-gPC system as explained in the previous section.

IV. APPLICATION TO THE KURAMOTO-SIVASHINSKY SYSTEM
The aforementioned methodology is now applied to the forced Kuramoto-Sivashinsky (KS) equation, where x ∈ [0, L] and boundary conditions u(0, t) = u(L, t) = 0 and u x (0, t) = u x (L, t) = 0. Two QoIs are defined, where J (1) (x, t) = u and J (2) (x, t) = u 2 , that represent time-and space-average values of the state and its energy.
The KS equation was discretised in space using secondorder accurate central finite differences and integrated with a variable time step Runge-Kutta method.We take L = 128 that results in chaotic behaviour.The state u(x, t) was stored every dt = 0.1 time units and used as input to the adjoint system.The first 1000 time units were discarded to ensure that the system has reached a chaotic attractor.The next T = 100 units were considered as the time horizon for the UQ analysis.
We first consider the bell-shaped forcing distribution ϕ(x) shown in figure 1, where Φ denotes the forcing amplitude.The smooth profiles from 0 → Φ (close to the left boundary) and from Φ → 0 (close to the right boundary) are obtained using the error function, erf , see [45] for more details.This forcing is infinitely differentiable and satisfies both Dirichlet and Neumann boundary conditions.
The variations of J (1) and J (2) with respect to amplitude Φ are shown in figures 2a and 3a respectively.The sensitivities of these QoIs obtained using finite differences and the adjoint of the shadowing harmonic operator are shown in figures 2b and 3b respectively.To form the harmonic operator, we consider frequencies in the range f ∈ [−0.3, 0.3] that captures the active frequency band of the unforced KS system, see [39].Notice that the two approaches are in very good agreement for both QoIs.In this case, where the sensitivity with respect to a single parameter is considered, the adjoint shadowing operator does not provide any computational advantage over the standard shadowing harmonic operator, hence this test case is only used as a benchmark to evaluate the accuracy and computational implementation of the method.
Contour plots of the state and adjoint variables in the (t, x) plane for the unforced system, i.e. for ϕ(x) = 0, FIG.2: Variation of (a) J (1) with the forcing amplitude Φ and (b) Sensitivity dJ (1)  dΦ computed using the adjoint operator and finite differences.Results are averaged over 200 random initial conditions.FIG.3: Variation of (a) J (2) with the forcing amplitude Φ and (b) Sensitivity dJ (2)  dΦ computed using the adjoint operator and finite differences.Results are averaged over 200 random initial conditions.
are shown in figure 4. Notice that the adjoint variable λ(x, t) does not have the same spatio-temporal streaky structure as the state variable u(x, t).This has been observed before in sensitivity analysis with shadowing method [35].Similarly to the state u(x, t), the spatiotemporal structure of the adjoint state is characterised by high sensitivity to initial conditions.We now evaluate the effect of stochastic variation of Φ to J (1) and J (2) using the se-gPC method and the sensitivities produced with the adjoint shadowing operator to augment the least squares system.We assume that Φ follows Gaussian distribution with Φ ∼ N (0, σ), and σ = 0.01.The standard deviation σ is small, but the response of the chaotic system to even small values of Φ is large, as shown in figures 2a and 3a.For example, within the range Φ ∈ [−3σ, +3σ], the value of J (2) doubles and the sensitivity varies between −100 to +100; this suggests that the effect of Φ is strong.In the range of Φ values considered the system is chaotic.Larger forcing amplitudes were also tested, however they result in forced oscillations with non-chaotic behaviour.
A Monte-Carlo simulation with 5000 samples was performed and used as a benchmark to evaluate the accuracy of the se-gPC method.For a single stochastic variable (m = 1) and polynomial order p = 1, there are P + 1 = 2 spectral coefficients, while for p = 2 there are P + 1 = 3 coefficients.In the se-gPC we used q = 4 samples for p = 1, that were augmented with another 4 se-gPC p = 1 se-gPC p = 2 Monte-Carlo QoI mean std mean std mean std J (1) 0.0196 0.6286 0.0196 0.6395 0.0195 0.6370 J (2) 2.3478 0.5251 2.3505 0.5174 2.3552 0.5192 TABLE I: Comparison between Monte Carlo simulations with 5000 samples and se-gPC for the KS system.The stochastic input is Φ ∼ N (0, 0.01).
equations for the sensitivity of the QoI with respect to Φ.For p = 2, q = 6 samples were used, augmented with 6 additional equations for the sensitivity.We used more equations than the number of coefficients to account for the (small) variation of the sensitivities to initial conditions.The results for the mean and standard deviation of the QoIs are summarized in table I, where se-gPC is compared with Monte Carlo (MC).There is very good agreement between the two methods.Errors in the standard deviation are less than 1.5% for p = 1 and less than 0.5% for p = 2 for both QoIs.
To better understand the effect of stochastic forcing ϕ(x), in figure 5 we plot the expectation of the time-average state, E J (1) (x) , where J (1) We compare against the profile of the unforced system (ϕ(x) = 0), where E J (1) (x) = J (1) (x).The results are averaged over 1000 random initial conditions.It's interesting to notice that the stochastic forcing smooths out the spatial oscillations of the time-average state of the unforced system.We also compute the standard deviation of J (1) (x) across x, and we superimpose the extent of one standard deviation above and below the expectation (the area between the two boundaries is marked grey).The large spread indicates that the actual time-average J 1 (x) oscillates wildly and can take values much larger that the expectation.Therefore the stochastic forcing drastically affects the output of the system.This explains why the standard deviation of J (1) is much larger that the mean (expectation) in table I.

V. MULTIDIMENSIONAL UNCERTAINTY FORCING
We now consider the stochastic forcing shown in figure 6, which is a continuous and differentiable profile that contains 10 peaks and troughs.The local amplitudes Φ i are the m = 10 independent stochastic variables considered.Such cases are usually found in control problems, where spatially complex forcing allows for more accurate control of the desired output quantities.The mean values of the localised forcing amplitudes are Φ 1 = 0.001, Φ 2 = −0.001,Φ 3 = 0.005, Φ 4 = 0.002, Φ 5 = 0.007, Φ 6 = −0.003,Φ 7 = −0.001,Φ 8 = −0.002,Φ 9 = 0.0005, Φ 10 = −0.002.These values were selected arbitrarily, since the main objective of the section is to demonstrate the efficiency of the proposed computational approach to conduct UQ.The standard deviations are taken equal to 20% of the mean value, i.e. σ Φi = |Φ i |/5.The spectra of λ(x, t) at x = L 4 , L 2 and 3L 4 are shown in fig.8 for the unforced and forced systems.The spectra were computed with a time window of T = 100 and were smoothed with a 5-th order Savitzky-Golay convolution filter with 5 averaging windows.For the unforced system, the spectra are very similar at the three locations.However, for the forced system the PSD values are larger and vary significantly between the locations due to the spatially varying forcing profile.The lower frequencies are also damped.The accuracy of the adjoint shadowing harmonic operator is assessed in figure 9.The values of the m = 10 sensitivities computed by the adjoint operator are compared against reference finite difference results.To evaluate the reference results we varied each Φ i separately and averaged over 100 initial conditions (we performed in total 1000 simulations with random initial conditions for all Φ i 's).The results for the adjoint shadowing operator were averaged over 100 random initial conditions.It is clear from the figure that the adjoint approach computes accurate sensitivities that are in very good agreement with finite differences for both J (1) and J (2) .

B. Uncertainty quantification with se-gPC
We now proceed to conduct UQ with se-gPC; the independent stochastic variables are the m = 10 amplitudes Φ i , as already mentioned.We first check if the system maintains its chaotic behaviour with the stochastic forcing.To this end, the system was integrated for 2000 random Φ i inputs over a time horizon of T = 200 and random initial conditions for each input.The Lyapunov exponents were computed for each realisation using the methodology presented in [46]; the maximum exponent λ max is shown in figure 10.For all realisations the forced system maintained its chaotic behaviour, with a maximum Lyapunov exponent that fluctuates around  FIG. 9: Comparison of the sensitivities of J (1) and J (2)  wrt to the amplitudes Φ i computed using finite differences and the adjoint shadowing operator.Results are averaged over 100 initial conditions for the adjoint shadowing operator.
the value of the unforced system (solid black line).For all systems it was observed that λ max > 0.04.In figure 11, the convergence rates of the mean and standard deviation of J (1) against the number of evaluations required by se-gPC, standard Weighted Least Squares (i.e.system (8)) and Smolyak Quadrature are compared.It was assumed that one adjoint solution has the same cost as a direct solution (i.e.forward integration of the dynamical system).This is a realistic assumption for the KS system we consider, but generally the cost of obtaining the adjoint solution for a chaotic system is case dependent.In the plot, one adjoint or one forward solution is considered as a single evaluation.The errors are FIG.11: Convergence rates of mean and standard deviation of J (1) against number of evaluations computed with WLS, Smolyak Quadrature and se-gPC for p = 2.The stochastic input is shown in figure 6.
The corresponding plot for J (2) is shown in figure 12. Once again, se-gPC coupled with the adjoint shadowing operator offers a significant computational advantage compared to other approaches.In practice this allows for the efficient quantification of uncertain input variables on the time-average quantities of chaotic systems.

VI. CONCLUSIONS
We propose a framework for efficient uncertainty quantification of time-average quantities of chaotic systems.We derive the adjoint version of the shadowing harmonic operator for sensitivity analysis of chaotic systems in the frequency domain.We subsequently employ the adjoint FIG.12: Convergence rates of mean and standard deviation of J (2) against number of evaluations computed with WLS, Smolyak Quadrature and se-gPC for p = 2.The stochastic input is shown in figure 6.
to compute the sensitivities of the QoI with respect to all uncertain variables and use this information to enrich the weighted least squares system from which the spectral coefficients of polynomial expansion are computed.The adjoint formulation provides all the required sensitivities in a single step, thus significantly increasing the computational efficiency of the method.The sampling points to integrate the dynamical system are obtained by the QR decomposition of an appropriate weighted matrix.The computational cost of the method is independent of the number of stochastic variables for polynomial order p = 1.
The adjoint formulation was applied to the Kuramoto-Sivashinsky equation and found to produce accurate sensitivities with respect to the amplitude of bell-shaped forcing.When these sensitivities were used to augment gPC, the resulting first and second moments computed matched excellently with Monte Carlo results.We then tested the method on a stochastically forced system with 10 independent input variables that determined the actual shape of the forcing.The adjoint was found to produce sensitivities that are in excellent agreement with finite differences and the se-gPC outperformed other UQ methods.
These attributes make the proposed method a promising candidate for application to chaotic systems with a large number of stochastic inputs.

FIG. 4 :
FIG. 4: Contour plots of the (a) state u(x, t) and (b) adjoint λ(x, t) variables for the unforced system, ϕ = 0. Results from a single realisation with a random initial condition.

FIG. 7 :
FIG.7: Contour plots of (a) the state variable and (b) the adjoint variable for the KS system forced with the profile shown in figure6.Results from a random initial condition.