Sensitivity analysis of chaotic systems using a frequency-domain shadowing approach

We present a frequency-domain method for computing the sensitivities of time-averaged quantities of chaotic systems with respect to input parameters. Such sensitivities cannot be computed by conventional adjoint analysis tools, because the presence of positive Lyapunov exponents leads to exponential growth of the adjoint variables. The proposed method is based on the least-square shadowing (LSS) approach [1], that formulates the evaluation of sensitivities as an optimisation problem, thereby avoiding the exponential growth of the solution. However, all existing formulations of LSS (and its variants) are in the time domain and the computational cost scales with the number of positive Lyapunov exponents. In the present paper, we reformulate the LSS method in the Fourier space using harmonic balancing. The new method is tested on the Kuramoto-Sivashinski system and the results match with those obtained using the standard time-domain formulation. Although the cost of the direct solution is independent of the number of positive Lyapunov exponents, storage and computing requirements grow rapidly with the size of the system. To mitigate these requirements, we propose a resolvent-based iterative approach that needs much less storage. Application to the Kuramoto-Sivashinski system gave accurate results with very low computational cost. The method is applicable to large systems and paves the way for application of the resolvent-based shadowing approach to turbulent flows. Further work is needed to assess its performance and scalability.


Introduction
Optimisation of engineering devices is based of the definition of an objective function, usually a time-average quantity J(s), and the evaluation of the problem parameters, s, that minimise or maximise this function, depending on the application. During the optimisation process, the gradient of the objective 5 function with respect to the parameters dJ(s)/ds (also known as sensitivity) is usually required. This is obtained by solving the tangent equation (or the adjoint equations for multiple parameters). In either case, the equations are obtained by linearising the governing non-linear set describing the system around the solution obtained for the reference values of the parameters, s. The tan-10 gent (or adjoint) equations are then integrated forward (or backward) in time respectively to obtain the desired sensitivities.
The aforementioned approach works very well when the governing set of equations describing the system is steady, in which case the solution is a point in phase space. When the evolution is however unsteady, and in particular 15 when the system exhibits chaotic behaviour, this process fails. The reason is that chaotic systems have one or more positive Lyapunov exponents (PLEs), thus two solution trajectories starting from the same initial conditions and evaluated at s and s + δs deviate from each other, leading to exponentially growing sensitivities, as explained in [2]. For this reason for example, model predictive 20 control algorithms for transitional or turbulent flows employ the receding horizon approach, whereby the optimisation is performed over a receding window of finite time [3,4], thus sensitivities remain bounded and reliable.
Several approaches that can compute useful sensitivities in chaotic systems have been proposed. These are based on ensemble schemes [5], the fluctuation 25 dissipation theorem [6], the Fokker-Planck equation [7], cumulant expansions [8], or unsteady periodic orbits [9]. One of the most promising approaches is Least Squares Shadowing (LSS) [1,10,11], which is based on the shadowing lemma [12,13]. For uniformly hyperbolic systems, this lemma guarantees the existence of a solution trajectory evaluated at s + δs that shadows, i.e. remains close to, the reference trajectory evaluated at s. This regularises the problem, avoids the exponential growth, and results in meaningful sensitivities.
This lemma is also central in establishing trust into the statistics of numerical solutions of chaotic systems. Due to round off errors, a computed trajectory will deviate from the true trajectory of the system (starting from the same initial 35 condition). The shadowing lemma guarantees the existence of a true trajectory (with different initial condition) that will shadow the numerical one [14, 15,16], thus the statistics of the computed solution can be trusted. Recent work however [17] has shown that the shadowing trajectories may not be physical, thus casting doubt on this central premise. Several chaotic systems (one dimensional 40 perturbed tent maps) were examined; the shadowing trajectories were found to be physical in one system and non-physical in others. This finding raises impor-  19,20] and the non-intrusive Least Squares Shadowing (NILSS), [21,22]. Both methods can be applied to large systems, but their computational cost scales with the number of PLEs. This limitation restricts the application 50 to systems with relatively small or moderate number of PLE's. For example, NILSS has been applied successfully to 2D flow (backward facing step with 14 PLEs [21]) and to 3D flows (minimal channel flow unit at Re τ = 140 with 150-160 PLEs [23], and flow around a cylinder at Re=525 with less than 30 PLEs [24]). 55 The number of PLEs and how it changes with the system parameters is therefore of critical importance for the application of the method. In the area of turbulent flows, the Lyapunov spectrum has been computed for low-Reynolds channel [25] and weakly chaotic Taylor-Couette flows [26]. A more recent study [27] has investigated the variation of the spectrum with Reynolds number for 60 forced homogeneous isotropic turbulence (HIT). The Reynolds numbers exam-ined (based on the Taylor microscale) were Re λ = 15.5, 21.3, 25.6 (note that these values are considered very small for engineering applications). The number of PLEs were found to be about 25, 60 and 100 respectively (see figure 2 of [27]). A fourth value of Re λ = 37.7 was also considered, but it was not possible to find the number of PLEs because of the slow decay rate of the spectrum.
The latter was found to follow a power-law, λ i ≈ α(i − 1) β + λ 1 , with the exponent β in the region 0.81-0.85 (with the smaller value for the higher Re τ ). The maximum LE, λ 1 , is expected to scale with the inverse of the Kolmogorov time scale, τ η , see theoretical arguments in [28]. This was tested in [29] for HIT and 70 it was found that λ 1 τ η is not constant, but instead grows with Re λ following a power-law, λ 1 τ η ∼ Re γ λ . The above scaling of λ 1 and the fact that the decay rate of the Lyapunov spectrum decreases with Reynolds number means that the number of PLEs grows very rapidly as Reynolds increases. Thus, alternative approaches are 75 required to make LSS and its variants applicable to complex flows of engineering interest. One such approach relies on the understanding of the underlying physical processes. For example, it is well known that momentum transfer is dominated by large scale structures, thus smaller scales (that are responsible for the largest LEs) can be filtered out and their effect modelled, hopefully 80 without loss of accuracy in sensitivity. This is exactly what large eddy simulations (LES) are designed for [30]. In standard LES, the equations are filtered in space, but temporal filtering (with filter time scale ∆) is also possible [31].
In the limit of ∆ → ∞, the temporally-averaged LES (TLES) equations tend to the standard Reynolds-Averaged Navier-Stokes (RANS) equations (section 85 2.3 of [31]). We conjecture therefore that as ∆ increases, application of LSS to the TLES equations will recover the sensitivites predicted by the tangent (or adjoint) method applied to RANS. Thus, parameter ∆ bridges two limits, LSS applied to unfiltered Navier-Stokes (∆ = 0) and the RANS equations (∆ → ∞).
As ∆ increases, the number of PLEs decreases and the problem becomes better 90 conditioned. However, accuracy is traded for computational efficiency, because as ∆ increases the effect of more scales needs to be modelled.
Although the above approach mitigates the rapid growth of the computational cost of LSS (and also provides a useful conceptual framework that bridges two limits), it relies on accurate modelling of the filtered scales. In the present 95 paper, we follow a different approach. All existing formulations of LSS and its variants have been derived in the time domain. If however the LSS is formulated in the frequency domain, the exponential separation of the trajectories for s and s + δs does not appear explicitly. Of course time-and frequency-domain formulations are equivalent, but as will be seen, the latter formulation allows us to 100 gain deep physical insight and also is amenable to iterative solution algorithms that are not possible with the former. Frequency domain approaches have been applied from the perspective of linear [32] and non-linear input-output analysis [33], the frequency response of periodically time-varying base flows [34], or model-based design of transverse wall oscillations for turbulent drag reduction in 105 a channel flow [35]. Similarities and differences with existing frequency-domain approaches are discussed throughout the manuscript.
The paper is organised as follows. Section 2 sets scene and presents the standard LSS algorithm in the time domain. The formulation of the algorithm in the frequency domain is derived in section 3 followed by application to the 110 Kuramoto-Sivasinsky equation in section 4. A resolvent-based iterative algorithm to solve the resulting system is presented in section 5 and the results are further analysed in section 6. We conclude in section 7.
2. Sensitivity analysis of chaotic systems using the shadowing approach 115 Consider a dynamical system governed by a set of ordinary differential equa- where u(t; s) ∈ R Nu is the vector of state variables and s ∈ R Ns is the set of control parameters that define the dynamics of the system. System (1) can arise for example after spatial discretisation of a set of conservation laws that describe mathematically the problem under investigation. We assume that the vector field f (u; s) : R Nu × R Ns → R Nu varies smoothly with u and s.

120
In many applications we are interested in evaluating the sensitivity of a time-averaged quantity J(s) : R Ns → R, to the parameters s. For example, in the area of aerodynamics, J(s) can be the drag coefficient and s the set of variables that describe the shape of an airfoil.
The gradient of J(s) with respect to s is defined as In chaotic systems, the limit and differentiation operations do not commute, i.e.
where v(t) = du ds = lim δs→0 u(t; s + δs) − u(t; s) δs is the sensitivity of the solution u(t; s) to a change δs of s. The reason is that chaotic systems have one or more PLEs. This means that the distance in phase space between u(t; s + δs) and u(t; s), i.e. the Euclidean norm of the vector u(t; s + δs) − u(t; s) that appears in the nominator of (5), grows exponentially at rate ∼ e λ1t , where λ 1 is the maximum of these exponents [2,7,36]. Thus, the 125 quantity 1 T T 0 ∂J ∂u v + ∂J ∂s that appears on the right hand side of (4) diverges as T → ∞. On the other hand, assuming that J(s) varies smoothly with s, the sensitivity dJ ds is finite. If the dynamical system (1) is uniformly hyperbolic, the shadowing lemma [37] guarantees the existence of a solution trajectory evaluated at s + δs that remains always close, i.e. shadows indefinitely, the reference trajectory u(t; s).
We denote this shadowing trajectory as u(τ (t); s + δs), where τ (t) is an appropriate time transformation. The LSS method, proposed in [1], computes the shadowing trajectory by minimising the distance between u(τ (t); s + δs) and u(t; s) in a least squares sense, i.e. min u,τ where α 2 is a constant parameter. Taking the limit δs → 0, leads to the following linear minimisation problem, where v(t) = d ds u(τ (t); s) (8a) The second term within the cost functions (6a) and (7a) penalises the deviation of τ (t) from t. A high value of α 2 results in a small deviation (heavy penalisation), while a small value to light penalisation. The solution of (7) for α 2 = 0, leads to the orthogonality condition between the vectors f (u; s) and v(u; s) at each point along the trajectory, i.e. f (u; s)v(u; s) = 0, a constraint from which η(t) can be obtained [18]. Thus, problem (7) becomes From the solution v lss (t) and η lss (t) of (9), the sensitivity dJ ds can be easily computed [1].

Formulation of the shadowing algorithm in Fourier space
As mentioned in the Introduction, all existing methods solve the minimisation problem (9) in the time domain. In this section, we formulate the problem in the frequency domain, i.e. in Fourier space, and seek a solution that remains bounded.

135
To this end, we consider a reference trajectory u(t; s) of length T and assume that the solution of the minimisation problem (9) is periodic with period T .
Thus, it can be written in terms of Fourier series as, wherev k ,η k denote the Fourier coefficients, ω 0 = 2π T is the fundamental angular frequency and the index k characterises the harmonics with frequencies ω k = kω 0 . We assume similar series expansions for the Jacobian ∂f ∂u (t) and f (t), The matrix-vector product ∂f ∂u (t)v(t) can be written as where which is the convolution sum between the Fourier coefficients of ∂f ∂u (t) and v(t). Similarly, the left hand side of the orthogonality condition (9c) can be expanded where and the notation () ⊤ denotes the transpose operation. In the above two expressions, we have assumed that the weighting matrix associated with the inner product is the identity matrix, but the analysis below can be easily generalised to an inner product defined as f (u; s)v(u; s) = f ⊤ Qv. Finally, In practise, the range of the index k is truncated to lie within the interval For example, if the reference trajectory is sampled every ∆t, q = T 2∆t −1. The frequency spectrum of u(t; s) can also indicate the number of spectral coefficients that must be retained. A finite q amounts to applying a sharp spectral cut-off filter to the above expansions, where all coefficients with |k|, |l| 140 or |k − l| > q are set equal to 0.
Introducing the finite spectral representations to (9b) and (9c), yields the following block set of equations for the k-th pair of coefficientsv k ,η k , where and I Nu is the identity matrix of dimension N u . Each block consists of N u + 1 equations and there in total 2q + 1 blocks, resulting in (N u + 1) × (2q + 1) equations and unknowns.
Due to the sharp spectral cut-off filter mentioned earlier, the starting and final values of index l in the summation is slightly modified when k = 0. For example, for k = −q equation (17) becomes while for k = +q, Stacking the blocks together one below the other results in a linear system of equations that takes the form with the same blocks in each diagonal.

145
System (22) can be also written in expanded matrix form as The above matrix, known known as Hill matrix [38,39], contains square blocks with dimensions (N u + 1) × (N u + 1), and thus has very large storage requirements. The solution of system (22) can be written symbolically as where is the matrix that maps the input R to the output V. As q → ∞, H becomes an operator, termed here shadowing harmonic operator. Note that equations (9b) and (9c) form a linear, time-varying periodic system,  Thus, the shadowing harmonic operator is identical to the standard harmonic operator, defined in [39], applied to the above system. The properties of this operator will be examined in the next section for the Kuramoto-Sivashinsky equation.
The sensitivity dJ ds is computed from [1], which in the frequency domain can be written as, In the above formulation, the reference trajectory u(t; s) was expanded using the same number of Fourier modes as the solution, 2q + 1, and each diagonal of the block Toeplitz matrix T (T m ) corresponds to one harmonic of u(t; s). This is however not necessary. For example, if the spectral content of u(t; s) is concentrated in a few frequencies, then only the relevant diagonals of T (T m ) need to be retained. This can lead to enormous savings in the storage requirements and solution time of system (24). In the limiting case, where only the timeaverage of u(t; s) is retained, the equations decouple and the k-th component In this case, the harmonic balancing method becomes identical to the standard 150 resolvent analysis [32].
Some comments are warranted here to clarify an underlying assumption of the above formulation. Suppose that the reference trajectory u(t; s), and thus ∂f ∂u (t) and f (t), are exactly periodic with period T and the minimisation problem (9) is solved in the time domain using the multiple shooting shadowing 155 method, [18,19]. If the trajectory is sampled at points t 0 . will not necessarily yield a periodic solution, i.e. v(t 0 ) will not necessarily be equal to v(t K ).
In the above formulation, we have not explicitly considered the minimisation of the cost function (9a); instead we closed the system assuming periodicity, i.e. v(t 0 ) = v(t K ) and η(t 0 ) = η(t K ). The same assumption is made in the 165 periodic shadowing method of Lasagna [40], and leads to a sensitivity error that initially decays at a rate 1/T , followed by the asymptotic rate 1/ √ T (the latter dictated by the central limit theorem). There is however an important difference compared to the present method: In [40], the time transformation τ (t) is linear with respect to t, leading to a constant η(t). In the present method, 170 we do not prescribe any form of τ (t), so η(t) is unknown and is obtained by imposing the orthogonality constraint (9c) at every point along the trajectory.
Another difference is that our method is formulated in the frequency domain instead of the time domain. As will be seen later, this can lead to significant simplifications, and allows one to obtain deep physical insight on the dominant factors that determine the sensitivity, which is not possible in the time domain.
An approach closely related to the present one was proposed recently by which is similar to (9b). The authors restrict the harmonic operator to a subspace that is orthogonal to the direction of the phase shift given byf k . This is achieved by projecting out ofĝ ′ k (the Fourier coefficients of the forcing g ′ (t)) the component that would lead to a non-zero projection of the solutionq ′ k tô f k . In the present paper, we seek the sensitivity with respect to s, thus the 185 forcing vector g ′ (t) takes the particular form g ′ (t) = ∂f ∂s (t). This vector is subsequently modified by adding η(t)f (t), where the time dilation η(t) is computed so that the orthogonality constraint is satisfied at all time instants, as already mentioned.
In the above frequency-domain formulation, we seek a harmonic solution that 190 remains bounded, see expansions (10), and does nor suffer from exponentially growing terms. There is a price to pay however; the frequencies are all coupled together leading to large storage requirements for the Hill matrix, see (24).
The properties of the shadowing harmonic operator will be examined in the next section for the Kuramoto-Sivashinsky equation. This test case is small 195 enough that the Hill matrix can be stored and the linear system (24) is solved directly with LU decomposition. In section 5 we propose an iterative method that mitigates the storage and solution time requirements.

Application to the Kuramoto-Sivashinsky equation
We apply the method proposed in the previous section to the Kuramoto Sivashinsky (KS) equation, which displays complex spatio-temporal chaos and is frequently used in the literature as a test case for chaotic systems [41]. The equation takes the form where c is an artificially introduced parameter [42]. The term ∂ 2 u ∂x 2 is responsible 200 for energy production, while ∂ 4 u ∂x 4 adds dissipation to the system. We set L = 128 to generate chaotic solutions [41] and discretise (31) with a second order finite difference scheme with δx = 1. For c = 0, the dynamical system has 16 positive Lyapunov exponents, the largest of which is λ 1 = 0.093 [42].
Two objective functions are considered, the space-time average of the state u(x, t; c) and of the total kinetic energy where an overbar () denotes time-average and a prime () ′ the fluctuation around 205 the average, i.e. u = u + u ′ . We seek their sensitivities with respect to c, i.e.
and As can be seen from figure 3, the sensitivities obtained using the shadowing harmonic operator match very well with those obtained using the preconditioned MSS [19]. A comparison with finite difference (FD) data is also presented. For dJ1 dc , there is a small bias, which has also been observed in the time-domain 230 formulation of the method [19,20,42]; for dJ2 dc the matching with FD is very good.  T . The difference initially decays at a rate 1 T (similarly to [40]) and then at  Thus for larger systems, both storage and solution costs grow fast.
Below we investigate in more detail the properties of the shadowing harmonic operator, while in the next section 5 we explore an approach that can mitigate the aforementioned rapid growth of computational cost and storage requirements for larger systems.

Singular value decomposition of the shadowing harmonic operator
The singular values σ i of the shadowing resolvent matrix H, defined in (26), are obtained from the solution of the eigenvalue problem, The solution maximises the system gain, defined as the ratio of the (squared) 2-norm of the output (response V) to that of the input (forcing R), i.e.
where V 2 2 = V * V and R 2 2 = R * R. Since δx = 1, the norms represent the discrete values of the integrals over the domain, for example turbulent flow around a cylinder see [24]. Note also that as T grows, apart from the largest singular values, the rest start to converge. This indicates that they represent the true behaviour of the system, i.e. they are physically meaningful.  i.e. Φ 1 , Φ 2 and Φ 3 , are shown in panels (b)-(d) respectively. The distributions are more difficult to interpret physically, but note the significant differences with 295 respect to the actual forcing. It is interesting to note for example that they do not exhibit the wavy pattern of df ds (x, t), instead they are relatively smooth but with some local peaks and valleys.
Using the r largest singular values of H, an approximate solution of system (24) can be written as where RΦ i = Φ * i R is the projection of the right hand side R onto the optimal forcing Φ i . For r = (2q + 1)(N u + 1) the approximation is exact. Using V (r) , can be computed from (34) and (35) respectively; these are plotted as functions of r in figure 9. It can be seen that a relatively large value of r is required to obtain an accurate result. Although the first few singular values σ i are large, the component of R along the Φ i direction is weak, thus a substantial number of terms are required to obtain the correct 305 sensitivity. Bearing in mind that R k = d f k ds , 0 ⊤ , this is related to the different patterns between df ds (x, t) and optimal forcings (at least for the first 3 modes) as shown in the previous figure 8.
The linear system (24) has large storage requirements and is time consuming to solve. Below we propose an approach that can mitigate these requirements.

A resolvent-based iterative method for the shadowing direction
Instead of solving directly system (24), an iterative method can be devised, where only a few diagonals are retained in the left hand side and treated implicitly, and the rest are moved to the right hand side and updated at every iteration. For example, by retaining only the blocks of the main diagonal, the Hill matrix becomes diagonal and the blocks decouple. In this case, the iterative method takes the from  where m is the iteration number, andĝ k ,ĥ k denote the explicitly treated terms.
It is instructive to derive the form of vectors g(t) and h(t) in the time domain.
To this end, using Reynolds decomposition, and substituting in (9b) and (9c) we get Taking the time-average we obtain and subtracting the two sets, we get After some rearrangement, Thus we get and taking the Fourier transform leads to (39) for k = 0. Similarly, the timeaverage system (42) can be written as from which we obtain the form that corresponds to k = 0 for system (39).
As can be seen from (44a), g ′ (t) consists of two groups of three terms; the first group involves the fluctuating Jacobian and sensitivity v ′ , and the second the fluctuating η ′ and f ′ . It is possible to get a simplified system by assuming that η is constant (for the physical interpretation see [40]). In this case we get This system requires an additional constraint to obtain the constant η. To this end, we require that f (t) and v(t) are perpendicular in a time-average sense, because f (t) is a real variable. Note that this condition couples together all Fourier components. Taking the Fourier transform of (48) we can form the following iterative method In the first iteration m = 1, we setg (0) k = 0 and we get These two equations can be combined together to obtain η (1) . Denoting the standard resolvent operator as solving forv (1) k and substituting in (51b) we get from which we obtain where Thus the solution of two linear systems that involve the standard resolvent operator, R(kω 0 ), is required. Due to the linearity of (50a), v (1) Applying inverse Fourier transform tov (1) k yields v (1) (t), from whichg(t) can be obtained from (48), and Fourier transformed to findg (1) k . The right hand side of (50a) can then be assembled and the second iteration performed. Note that 315 system (56b) does change with m, thus µ k is computed once. In this process, the key variables λ k , η, and µ k required for the evaluation of the sensitivitŷ v k (also known as shadowing direction) are obtained with the aid of the Resolvent operator R(kω 0 ). This is therefore a resolvent-based iterative method for computing the shadowing direction, called Resolvent-based Shadowing (RbS).
The most efficient approach is to perform LU decomposition of R(kω 0 ) once at the start of the sensitivity analysis and solve systems (56) with forward and backward substitution for each k. Note that the decomposition (57) ofv k into a linear combination of λ k and µ k is valid for all iterations, and of course the final converged solution.

325
In (48a), the term ∂f ′ ∂u v was considered as part ofg ′ (t) and treated explicitly. This is not necessary, but it simplifies the algebra. A better approach would be to obtain v from (48b) and substitute in (48a); this would lead to a form very similar to (51a), albeit more complex. The process to extract η remains the same. Such an approach would couple better the time-average and the very efficiently using existing linear algebra packages, such as MUMPS [43], that exploit the sparse structure of the Jacobian. For three-dimensional inhomogeneous flows this is more challenging, but doable (at least for moderate-size systems).
The proposed method was applied to compute the sensitivities of the objective functions J 1 (c) and J 2 (c), defined in (32) and (33), for the Kuramoto-Sivashinsky equation. Since η is constant, equations (34) and (35) are simplified and In terms of computational cost, it took approximately 0.15s of CPU time to evaluate these sensitivities; this is about two orders of magnitude faster compared to preconditioned MSS and one order faster with respect to the shadow-ing harmonic approach. Again, we stress that these results are case dependent.
More research is need to investigate the performance of the algorithm in other flow cases.
In the next section, we explore the properties of the standard resolvent operator R(kω 0 ). For angular frequencies ω > 2 π 0.3 ≈ 1.9 there is no coupling with the timeaverage Jacobian, leading to σ(ω) ≈ 1/ω (indicated by a black dashed line); again this is consistent with the frequency spectra.  Substituting expression (57) into (58) we get The solutions λ 0 and µ 0 of the linear systems (56a) and (56b) can be written in terms of the optimal forcings and responses evaluated for ω = 0 as where r is the number of retained singular values. This expression is analogous to equation (38) presented earlier for the harmonic balance method. Substitut-ing in (60) we get The sensitivity can therefore be written as the weighted sum of the spatial averages of the optimal responses, Ψ i .

385
Similarly for dJ2 dc , we get from (59), where Expression (63) is useful because it allow us to find the contribution of each frequency on the sensitivity; we investigate this in figure 13a. More specifically, we compute dJ2 dc using frequencies in the interval [−f, f ], where f = ω/(2π), and we plot the result against f . Note the convergence of dJ2 dc to the value predicted by the harmonic resolvent as f increases, i.e. the range of k in (63) expands.

390
It can be seen that only the frequency range [0.01, 0.07] contributes, which is consistent with the spectra shown in figure 2. Outside this range, the spectral content is small, and does not contribute to the sensitivity. Note that this is the result of a single realisation; no averaging over initial conditions has been performed to obtain this plot. Hence, although the maximum singular value is significantly larger compared with the rest as evidenced in 11, keeping just one contribution will not provide 400 accurate results. In order to explain this behaviour for dJ1 dc , the optimal forcing corresponding to σ max at ω = 0 is plotted together with the true forcing ∂f