Leading order response of statistical averages of a dynamical system to small stochastic perturbations
LLEADING ORDER RESPONSE OF STATISTICAL AVERAGES OF ADYNAMICAL SYSTEM TO SMALL STOCHASTIC PERTURBATIONS
RAFAIL V. ABRAMOVA bstract . The classical fluctuation-dissipation theorem predicts the average response ofa dynamical system to an external deterministic perturbation via time-lagged statisticalcorrelation functions of the corresponding unperturbed system. In this work we developa fluctuation-response theory and test a computational framework for the leading orderresponse of statistical averages of a deterministic or stochastic dynamical system to anexternal stochastic perturbation. In the case of a stochastic unperturbed dynamical sys-tem, we compute the leading order fluctuation-response formulas for two different cases:when the existing stochastic term is perturbed, and when a new, statistically indepen-dent, stochastic perturbation is introduced. We numerically investigate the effectivenessof the new response formulas for an appropriately rescaled Lorenz 96 system, in both thedeterministic and stochastic unperturbed dynamical regimes.
1. I ntroduction
Under suitable conditions, the fluctuation-dissipation theorem (FDT) [28, 29, 41] fur-nishes an approximation to the statistical response of a dynamical system to a deter-ministic external perturbation via statistical correlations of the unperturbed dynamics.The FDT offers more insight into statistical properties of dynamical processes near equi-librium in various scientific applications [13, 15, 16, 18, 21–25, 30, 32, 37, 38]. In the pastworks [1–3, 7, 8, 10–12], a computational framework predicting the average response ofboth deterministic and stochastic dynamical systems to a small deterministic externalperturbation was developed and extensively studied.In the current work, we develop the fluctuation-response theory and numerically testthe computational framework of the response of statistical averages to a stochastic exter-nal perturbation, for both the deterministic and stochastic unperturbed dynamics. Forthe deterministic unperturbed dynamics, our set-up is similar to the one used in [35],however, the resulting formula we arrive at is different from the one in [35], on whichwe will comment below. For the stochastic unperturbed dynamics, we consider two dif-ferent types of perturbations: first, where the existing stochastic term is perturbed, and,second, when a new, statistically independent, stochastic term is introduced. We test thecomputational framework of the stochastic response on the Lorenz 96 system [33, 34],which we used as a test-bed nonlinear chaotic system with forcing and dissipation forvarious purposes in the past [1–11, 38]. D epartment of M athematics , S tatistics and C omputer S cience , U niversity of I llinois at C hicago , 851 S. M organ st ., C hicago , IL 60607 E-mail address : [email protected] . Date : September 17, 2018. a r X i v : . [ n li n . C D ] D ec RAFAIL V. ABRAMOV
Before going into the details, here we start by explaining the basic idea of the averageresponse and how it can be expressed via the statistical properties of the underlyingunperturbed dynamical system, and also what problems one runs into while consideringwhat otherwise seems to be a rather simple dynamical set-up.1.1.
Deterministic dynamics.
We start by considering a system of ordinary differentialequations of the form(1.1) d x t d t = f ( x t ) ,where t is a scalar time variable, x t is an N -dimensional vector in the Euclidean space R N ,and f : R N → R N is a smooth vector field. Observe that the solution x t can be specifiedin the form of a semigroup φ t ,(1.2) x t = φ t x = x + (cid:90) t f ( x s ) d s ,where x is the initial condition. We assume that any solution x t = φ t x of (1.1) is attracted,as t → ∞ , to a compact set M ⊂ R N , on which it possesses a unique invariant ergodicmeasure µ . We will say that M is the global attractor of (1.1). Here we assume that thesystem in (1.1) is chaotic and mixing, that is, it has positive first Lyapunov exponent anddecaying time autocorrelation functions.Let A ( x ) be a twice-differentiable function on R N , then we denote its µ -average as(1.3) (cid:104) A (cid:105) = (cid:90) M A ( x ) d µ ( x ) .Observe that even though A ( x t ) changes with time t , its µ -average (cid:104) A ( x t ) (cid:105) is fixed in t ,due to the fact that µ is invariant on M under φ t ,(1.4) (cid:104) A ( x t ) (cid:105) = (cid:90) M A ( φ t x ) d µ ( x ) = (cid:90) M A ( x ) d µ ( x ) = (cid:104) A (cid:105) .1.2. The concept of the average response to a stochastic perturbation.
Consider thesituation where the average (cid:104) A (cid:105) is computed across a statistical ensemble of solutionsof (1.1), which is distributed according to the invariant measure µ above. As we alreadypointed out, this average (cid:104) A (cid:105) is constant in time. However, assume that, at t = µ is nolonger invariant for the new, modified dynamics. Because of that, the statistical average (cid:104) A (cid:105) with respect to µ becomes time-dependent for the perturbed dynamics. Here we donot assume that µ is necessarily the Sinai-Ruelle-Bowen measure [45] as the finite timeresponse to an external perturbation does not require it; an SRB measure is, however,necessary for the infinite time response to be differentiable with respect to a deterministicexternal perturbation [42, 43].The difference between the new time-dependent ensemble average of A and its previ-ous stationary (for the unperturbed dynamics) value is then called the “response”:(1.5) Response of A = δ (cid:104) A (cid:105) ( t ) = (cid:104) A (cid:105) new ( t ) − (cid:104) A (cid:105) old . EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 3
Observe that above the average is taken not only with respect to the statistical ensembleof solutions, but also over all possible realizations of the external stochastic perturba-tion. Our goal here is to derive the leading order term in the response, which dependsonly on statistics of the unperturbed dynamical system, under the assumption that theperturbation is sufficiently small.The more obvious, “brute force” approach, would be to do the following:1. Start with a point x on M , and emit two trajectories out of x : the unperturbed one,given by φ t x , and the perturbed one, given by the corresponding solution of the per-turbed system.2. Clearly, both trajectories, perturbed and unperturbed, are generally nonlinear func-tions of elapsed time t and initial condition x . So, assuming that the unperturbedsolution φ t x is “known”, figure out a suitable way to “expand” the perturbed solutionaround the unperturbed one in small increments, and keep only the leading orderterm.3. Recall that this has to be done for every x ∈ M , so, average the above result out withrespect to the invariant measure µ , and over all possible realizations of the externalstochastic perturbation. Provided that the leading order response from the previousitem is somehow expressed in terms of trajectories of the unperturbed system, the µ -average can be replaced with the long-term time average, with help of Birkhoff’stheorem [14].With the exception of averaging across realizations of stochastic perturbations, this iswhat was previously done for the deterministic perturbations of chaotic and stochasticdynamical systems [1–3, 7, 10–12]. It is a long and cumbersome way of deriving theresponse, and the result involves the (computationally expensive) tangent map T t x , givenby(1.6) T t x = ∂∂ x φ t x .The situation is further complicated by the fact that, for chaotic dynamical systems, T t x grows exponentially fast in t , which causes a numerical instability for moderate andlong response times. Thus, this approach can only be practically used for rather shortresponse times (although it is usually quite precise, provided that the response time issufficiently short [1–3, 7]).1.3. The forward Kolmogorov equation.
Another way to compute the average responseis to employ the concept of the probability density p of a statistical ensemble distribution.The key idea here is to use that fact that, while x t is governed by nonlinear dynamics, thepartial differential equation for p (called the forward Kolmogorov equation [19], and alsothe Fokker-Planck equation [41]) is linear. In particular, for the deterministic dynamicalsystem in (1.1), the forward Kolmogorov equation for p is given by(1.7) ∂∂ t p ( t , x ) = − D · ( p ( t , x ) f ( x )) , RAFAIL V. ABRAMOV where D is the differentiation operator with respect to the vector-argument of the func-tion it acts upon. Above, the dot-product of D with a vector-function a ( x ) refers to(1.8) D · a ( x ) = ∂ a i ( x ) ∂ x i ,with the usual summation convention. Observe that in order for the solution p of (1.7)to remain a probability density, its integral over R N must remain equal to 1 (which,together with the non-negativity of p , implies that p must vanish at infinity), even if p itself changes with time. Indeed, one can verify that the integral of p over R N ispreserved by (1.7), which is necessary for p to remain a probability density.The Kolmogorov equation above is an extremely useful tool for working with thestatistical properties of the system in (1.1), since it describes the statistical distribution ofthe system in a direct fashion. Unfortunately, it cannot be used directly to compute theresponse of the deterministic system in (1.1), for the following reason.Since, as stated earlier, any solution of (1.1) attracts to M as t becomes infinite, itwould be natural to think that, in the limit as t → ∞ , p becomes the density of the ergodicinvariant measure µ on M . However, here lies the fundamental “incompatibility” of theKolmogorov equation in (1.7), and the limiting dynamics of (1.1) on its global attractor M : for many applied dynamical systems, especially those with dissipation and forcing[17, 42, 43, 45], the invariant measure µ on M is not differentiable in x (it is also saidthat it is not continuous with respect to the Lebesgue measure on M ). In this situation,the (non-stationary) solution p ( t , x ) of (1.7) contracts exponentially rapidly along certaindirections of the phase space (while appropriately expanding transversally, so that itsintegral over R N remains 1), becoming singular in the infinite time limit.Observe that above we considered arguably the most simple setup for a dynamicalsystem, which describes a wide class of applied problems. Yet, we cannot make use ofthe Kolmogorov equation (1.7) to statistically describe dynamics near the attractor of thesystem in (1.1), which is necessary for understanding of how the system responds to anexternal perturbation. Therefore, in order to be able to use the Kolmogorov equationin (1.7), we must be willing to consider a suitable modification of (1.1), which renders itsinvariant measure µ continuous with respect to the Lebesgue measure on R N . Arguably,the simplest such modification is achieved via a stochastic noise added into the otherwisedeterministic dynamical system in (1.1).1.4. Stochastic dynamics.
Here we are going to consider a stochastic modification of (1.1),achieved via introducing an additional noise term via a Wiener process W t of dimension K :(1.9) d x t = f ( x t ) d t + G ( x t ) d W t ,where G ( x t ) is a smooth N × K matrix. For convenience, here we interpret the resultingintegral of the solution(1.10) x t = x + (cid:90) t f ( x s ) d s + (cid:90) t G ( x s ) d W s EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 5 in the sense of It ˆo [26,27]. The forward Kolmogorov equation for the differential equationin (1.9) is given by [19, 39](1.11) ∂ p ∂ t = − D · ( p f ) + D : ( p GG T ) ,where “:” denotes the Frobenius product of two matrices, so that(1.12) D : A ( x ) = ∂ A ij ( x ) ∂ x i ∂ x j for an N × N matrix A .To ensure the smoothness of solutions of (1.11), here we follow [40] and assume thatboth f and G have bounded derivatives of all orders, and that the matrix product GG T is uniformly positive definite in R N . The latter automatically means that the columns of G span R N for any x ∈ R N , implying K ≥ N . Additionally, we will assume that there is aunique smooth stationary probability density p which sets the right-hand side of (1.11)to zero.Observe that the solution x t of (1.9) cannot be represented by a t -semigroup likein (1.2), since W t depends on t explicitly. However, instead a similar representationcan be done for the Kolmogorov equation in (1.11) with help of the transitional proba-bility density p ∗ . Indeed, let p ∗ ( t , x , x ) denote the solution of (1.11), for which the initialcondition at t = δ ( x − x ) . Then, assuming that at time t thesolution is p ( t , x ) , its extension p ( t + s , x ) for s ≥ p ∗ as follows:(1.13) p ( t + s , x ) = (cid:90) R N p ∗ ( s , x , y ) p ( t , y ) d y def = P s p ( t , x ) .Now, let p denote the stationary smooth probability density of (1.11), such that(1.14) 12 D : ( p GG T ) − D · ( p f ) = p = P t p for any t ≥ A ( x ) is given by(1.16) (cid:104) A (cid:105) = (cid:90) R N A ( x ) p ( x ) d x ,where we assume that A ( x ) is such that the integral above is finite. As before, if astatistical ensemble of solutions x t of (1.9) is distributed according to p , then the ensem-ble average of A is constant in time, even though each individual solution in such anensemble is in general not stationary.1.5. The layout of the paper.
In what is to follow, we arrange the presentation in thereverse order (relative to what was presented above), due to the fact that, as was men-tioned previously, it turns out to be easier to start with a stochastic differential equationof the form in (1.9) and derive the leading order response via the Kolmogorov equa-tion (1.11), which we do in Section 2. Then, in Section 3 we return to the deterministicunperturbed dynamics, and derive the response formula in the “brute force” fashion,
RAFAIL V. ABRAMOV sketched above. In Section 4 we show that, if one formally replaces the invariant mea-sure µ of the deterministic system in (1.1) with a smooth density approximation, then theresponse formulas for the deterministic and stochastic unperturbed dynamics becomeequivalent. In Section 5 we derive simplified response formulas for both the determin-istic and stochastic unperturbed dynamics under the assumption that the probabilitydensity of the unperturbed state is Gaussian, as was previously done in [9–12, 38] for de-terministic perturbations. In Section 6 we present the numerical experiments with boththe deterministic and stochastically forced Lorenz 96 models to verify the computedresponse formulas. Section 7 summarizes the results of the work.2. L eading order response of a stochastic dynamics to a stochasticperturbation We start by considering the response of the stochastic dynamics in (1.9), as due to thefact that the invariant state p of the unperturbed dynamics in (1.9) is a smooth stationarysolution of the Kolmogorov equation in (1.11), it is in fact much easier technically toderive the corresponding leading order response formula (as opposed to the situationwith a deterministic unperturbed dynamics). Here we will consider two different typesof perturbation: first, when the existing stochastic matrix is perturbed, and, second,when a new, statistically independent stochastic perturbation is added.2.1. Perturbing the existing stochastic term.
First, we are going to assume that thealready present in (1.9) stochastic diffusion matrix G ( x ) is perturbed by a small time-dependent term starting at t = x t = f ( x t ) d t + ( G ( x t ) + εη ( t ) H ( x t )) d W t ,where 0 < ε (cid:28) H ( x ) is a matrix of the same dimension and smoothness properties as G , and η ( t ) is abounded, piecewise continuous and square-integrable function, which is zero for neg-ative values of t . Then, the corresponding perturbed Kolmogorov equation is obtainedfrom (1.11) by replacing G with G + εη H :(2.2) ∂ p ε ∂ t = − D · ( p ε f ) + D : (cid:16) p ε ( G + εη ( t ) H )( G + εη ( t ) H ) T (cid:17) .We assume that the solution p ε of the perturbed Kolmogorov equation above dependssmoothly on ε for sufficiently small ε , and admits the expansion(2.3) p ε = p + ε p + ε p + . . . ,where p is the stationary solution of the unperturbed Kolmogorov equation, while theexpansion terms p i , i > ε and have zero initial conditions. We seekthe perturbed solution in the leading order of ε , which leads directly to(2.4) ∂ p ∂ t = − D · ( p f ) + D : ( p GG T ) + η ( t ) D : (cid:16) p ( GH T + HG T ) (cid:17) .One can verify directly that the solution for p is given by(2.5) p ( t ) = (cid:90) t P t − s (cid:16) D : ( p ( GH T + HG T )) (cid:17) η ( s ) d s , EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 7 where P t is defined in (1.13). The response of A in the leading order of ε is thus givenby(2.6) δ (cid:104) A (cid:105) ( t ) = (cid:90) R N A ( p − p ) d x = ε (cid:90) R N Ap d x + O ( ε ) = ε (cid:90) t R ( t − s ) η ( s ) d s + O ( ε ) ,with(2.7) R ( t ) = (cid:90) R N A ( x ) P t (cid:16) D : ( p ( GH T + HG T )) (cid:17) ( x ) d x .At this point, we use the definition of P t in (1.13) to obtain R ( t ) in terms of the time-correlation function [40, 41](2.8) R ( t ) = (cid:90) R N (cid:90) R N A ( x ) p ∗ ( t , x , y ) B ( y ) p ( y ) d y d x ,with(2.9) B ( y ) = p − ( y ) D : ( p ( GH T + HG T ))( y ) .Above, the division by p is allowed since it is the solution to an elliptic equation whichvanishes at infinity, and thus is never zero for finite y .For the practical computation of the time correlation function in (2.8), we use theBirkhoff-Khinchin theorem [20] and replace the spatial integrals in (2.8) with the timeaverage along a long-term trajectory as(2.10) R ( t ) = lim r → ∞ r (cid:90) r A ( x t + s ) B ( x s ) d s ,which results, after substituting the expression for B ( x s ) , in(2.11) R ( t ) = lim r → ∞ r (cid:90) r A ( x t + s ) p − ( x s ) D : (cid:16) p ( x s )( G ( x s ) H T ( x s ) + H ( x s ) G T ( x s )) (cid:17) d s .The long-term trajectory x s above is computed via a numerical simulation of (1.9), froman arbitrary initial condition.2.2. Adding a new stochastic term.
Here we assume that a new small stochastic term isadded to (1.9),(2.12) d x t = f ( x t ) d t + G ( x t ) d W t + εη ( t ) H ( x t ) d W (cid:48) t ,where the Wiener process W (cid:48) t is independent of W t . In order to derive the correspondingKolmogorov equation, we rewrite the above equation in the form(2.13) d x t = f ( x t ) d t + (cid:101) G ( t , x t ) d (cid:102) W t ,where(2.14) (cid:102) W t = (cid:18) W t W (cid:48) t (cid:19) , (cid:101) G ( t , x t ) = (cid:16) G ( x t ) (cid:12)(cid:12)(cid:12) εη ( t ) H ( x t ) (cid:17) .The corresponding perturbed Kolmogorov equation is, obviously, given by(2.15) ∂ p ε ∂ t = − D · ( p ε f ) + D : ( p ε (cid:101) G (cid:101) G T ) . RAFAIL V. ABRAMOV
Further observing that(2.16) (cid:101) G (cid:101) G T = GG T + ε η ( t ) HH T ,we arrive at(2.17) ∂ p ε ∂ t = − D · ( p ε f ) + D : ( p ε GG T ) + ε η ( t ) D : ( p ε HH T ) .Observe that there is no first-order term in ε , so we can expand p ε near p in even powersof ε as(2.18) p = p + ε p + ε p + . . . ,where, as before, p i for i > ε and have zero initial condition, andsimilarly obtain the equation for p as(2.19) ∂ p ∂ t = − D · ( p f ) + D : ( p GG T ) + η ( t ) D : ( p HH T ) .Again, one can verify that p is given by(2.20) p ( t ) = (cid:90) t P t − s D : ( p HH T ) η ( s ) d s .The response of A in the leading order of ε (which is now ε ) is thus given by(2.21) δ (cid:104) A (cid:105) ( t ) = (cid:90) R N A ( p − p ) d x = ε (cid:90) R N Ap d x + O ( ε ) = ε (cid:90) t R ( t − s ) η ( s ) d s + O ( ε ) ,with(2.22) R ( t ) = (cid:90) R N A ( x ) P t (cid:16) D : ( p HH T ) (cid:17) ( x ) d x .Following the same steps as above for R ( t ) , we express R ( t ) via the time average as(2.23) R ( t ) = lim r → ∞ r (cid:90) r A ( x t + s ) p − ( x s ) D : (cid:16) p ( x s ) H ( x s ) H T ( x s ) (cid:17) d s .Observe that in this case the response is quadratic in εη ( t ) (which is unlike the previouscase, where the existing diffusion matrix was perturbed).3. L eading order response of a deterministic dynamics to a stochasticperturbation Now, we consider a small external stochastic perturbation of the deterministic systemin (1.1) of the form(3.1) d x t = f ( x t ) d t + εη ( t ) H ( x t ) d W t ,where the perturbation term has the same properties as in the previous section, while,for the purposes of the derivation, we additionally require f to be uniformly Lipschitzin R N [19, 40] to ensure the existence of solutions to (3.1) (recall that the unperturbedsystem back in (1.1) does not necessarily require it for global existence). Here, however,we cannot use an expansion of Kolmogorov equation near the stationary unperturbedstate, since this state may not necessarily be continuous with respect to the Lebesgue EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 9 measure. Instead, we will have to employ the differentiability of the resulting stochasticflow with respect to the perturbation [31].We denote the solution to the perturbed system in (3.1) by x ε t = φ ε t x , and rewrite (3.1)in the integral form as(3.2) φ ε t x = x + (cid:90) t f ( φ ε s x ) d s + ε (cid:90) t η ( s ) H ( φ ε s x ) d W s ,where the stochastic integral is computed in the sense of It ˆo. Let A ( x ) be a twice differ-entiable function, then one can write its second-order Taylor expansion in ε as(3.3) A ( φ ε t x ) − A ( φ t x ) = ε D A ( φ t x ) · ∂ ε φ t x ++ ε (cid:104) D A ( φ t x ) · ∂ ε φ t x + D A ( φ t x ) : ( ∂ ε φ t x ⊗ ∂ ε φ t x ) (cid:105) + o ( ε ) ,where “ ⊗ ” is the outer product of two vectors, that is,(3.4) x ⊗ y = x i y j .Also, the following notation is used above:(3.5) ∂ ε φ t x def = ∂φ ε t x ∂ε (cid:12)(cid:12)(cid:12)(cid:12) ε = .For the ε -derivative of φ ε t x (which we assume to exist almost surely for finite t accordingto [31]) we compute(3.6) ∂∂ε φ ε t x = (cid:90) t D f ( φ ε s x ) ∂∂ε φ ε s x d s + ε (cid:90) t η ( s ) D H ( φ ε s x ) ∂∂ε φ ε s x d W s ++ (cid:90) t η ( s ) H ( φ ε s x ) d W s ,which results, by setting ε =
0, in(3.7) ∂ ε φ t x = (cid:90) t D f ( φ s x ) ∂ ε φ s x d s + (cid:90) t η ( s ) H ( φ s x ) d W s .At this point, we need to solve the It ˆo integral equation above. Applying the It ˆo differ-entiation formula to both sides of (3.7) results in(3.8) d ( ∂ ε φ t x ) = D f ( φ t x ) ∂ ε φ t x d t + η ( t ) H ( φ t x ) d W t .At the same time, it is easy to verify that the tangent map T t x from (1.6) satisfies(3.9) ∂∂ t T t x = D f ( φ t x ) T t x ,which further yields(3.10) d ( ∂ ε φ t x ) = (cid:18) ∂∂ t T t x (cid:19) (cid:0) T t x (cid:1) − ∂ ε φ t x d t + η ( t ) H ( φ t x ) d W t . Now we multiply both sides of the above identity by the inverse of T t x on the left, whichresults, after taking into account the identity = A − ∂ I ∂ t = A − ∂∂ t ( AA − ) = A − ∂ A ∂ t A − + ∂∂ t ( A − ) for an arbitrary matrix A , in(3.11) (cid:0) T t x (cid:1) − d ( ∂ ε φ t x ) = − (cid:18) ∂∂ t (cid:0) T t x (cid:1) − (cid:19) ∂ ε φ t x d t + η ( t ) (cid:0) T t x (cid:1) − H ( φ t x ) d W t .Pulling the first term in the right-hand side above to the left, combining the terms andintegrating from 0 to t , we arrive at(3.12) (cid:0) T t x (cid:1) − ∂ ε φ t x = (cid:90) t η ( s ) ( T s x ) − H ( φ s x ) d W s + (cid:104)(cid:0) T t x (cid:1) − , ∂ ε φ t x (cid:105) t ,where the last term is the quadratic covariation of the processes (cid:0) T t x (cid:1) − and ∂ ε φ t x :(3.13) (cid:104)(cid:0) T t x (cid:1) − , ∂ ε φ t x (cid:105) t = (cid:90) t d (cid:16) ( T s x ) − (cid:17) d ( ∂ ε φ s x ) = (cid:90) t ∂∂ s (cid:16) ( T s x ) − (cid:17) d s d ( ∂ ε φ s x ) .However, since we have assumed above that ∂ ε φ t x = (cid:82) t d ( ∂ ε φ s x ) is almost surely finitefor finite t , the quadratic covariation above in (3.13) is almost surely zero. Further mul-tiplying (3.12) by T t x on the left and taking into account its cocycle property, we finallyarrive at(3.14) ∂ ε φ t x = (cid:90) t η ( s ) T t − s φ s x H ( φ s x ) d W s .For the second ε -derivative we further obtain by the differentiation of (3.6)(3.15) ∂ ∂ε φ ε t x = (cid:90) t (cid:20) D f ( φ ε s x ) ∂ ∂ε φ ε s x + D f ( φ ε s x ) : (cid:18) ∂∂ε φ ε s x ⊗ ∂∂ε φ ε s x (cid:19)(cid:21) d s ++ (cid:90) t (cid:20) η ( s ) D H ( φ ε s x ) ∂∂ε φ ε s x + εη ( s ) ∂∂ε (cid:18) D H ( φ ε s x ) ∂∂ε φ ε s x (cid:19)(cid:21) d W s ,which becomes, upon setting ε = ∂ ε φ t x = (cid:90) t (cid:104) D f ( φ s x ) ∂ ε φ s x + D f ( φ s x ) : ( ∂ ε φ s x ⊗ ∂ ε φ s x ) (cid:105) d s ++ (cid:90) t η ( s ) D H ( φ s x ) ∂ ε φ s x d W s .Now recall that, according to the definition of the average response in (1.5), we shouldcompute the expectation (that is, the average) of the second-order Taylor expansionin (3.3) over all realizations of W t , which leads to(3.17) E A ( φ ε t x ) − A ( φ t x ) = ε D A ( φ t x ) · E ∂ ε φ t x ++ ε (cid:104) D A ( φ t x ) · E ∂ ε φ t x + D A ( φ t x ) : E ( ∂ ε φ t x ⊗ ∂ ε φ t x ) (cid:105) + o ( ε ) . EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 11
One immediately observes that(3.18a) E ∂ ε φ t x = E (cid:90) t η ( s ) T t − s φ s x H ( φ s x ) d W s = E (cid:90) t η ( s ) D H ( φ s x ) ∂ ε φ s x d W s = W t . After some computation whileremembering Duhamel’s principle and It ˆo’s isometry, we arrive at(3.19a) E ( ∂ ε φ t x ⊗ ∂ ε φ t x ) = (cid:90) t T t − s φ s x H ( φ s x ) H T ( φ s x ) (cid:16) T t − s φ s x (cid:17) T η ( s ) d s ,(3.19b) E ∂ ε φ t x = (cid:90) t (cid:90) s T t − s φ s x D f ( φ s x ) : (cid:20) T s − τφ τ x H ( φ τ x ) H T ( φ τ x ) (cid:16) T s − τφ τ x (cid:17) T (cid:21) η ( τ ) d τ d s .Now we recall that the expectations above are taken under the condition that the sto-chastic flows start at x , which, in turn, belongs to the attractor of (1.1). Therefore, wefurther need to average the result above over the invariant measure µ of the unperturbedsystem. The result is, after discarding the higher-order terms,(3.20a) δ (cid:104) A (cid:105) ( t ) = ε (cid:90) t R ( t − s ) η ( s ) d s + o ( ε ) ,(3.20b) R ( t ) = (cid:90) M (cid:20) D A ( φ t x ) : (cid:16) T t x H ( x ) H T ( x ) (cid:0) T t x (cid:1) T (cid:17) ++ D A ( φ t x ) (cid:90) t T t − s φ s x D f ( φ s x ) : (cid:16) T s x H ( x ) H T ( x ) ( T s x ) T (cid:17) d s (cid:21) d µ ( x ) .4. T he equivalence of the response formulas for the deterministic andstochastic dynamics While the response formula in (3.20b) is rather complicated for a practical use, one canactually show that the response operator R ( t ) in (3.20b) can be written in a more conciseway:(4.1) R ( t ) = (cid:90) M ∂ A ( φ t x ) ∂ x : (cid:0) H ( x ) H T ( x ) (cid:1) d µ ( x ) .Indeed, observe that, first,(4.2a) ∂ A ( φ t x ) ∂ x = D A ( φ t x ) T t x ,(4.2b) ∂ A ( φ t x ) ∂ x = D A ( φ t x ) : (cid:0) T t x ⊗ T t x (cid:1) + D A ( φ t x ) ∂∂ x T t x ,where the combination of the Frobenius and outer product for matrices denotes(4.3) (cid:2) A : ( B ⊗ C ) (cid:3) ij = A kl B ki C lj . Next, differentiating (3.9) with respect to x yields(4.4) ∂∂ t (cid:18) ∂∂ x T t x (cid:19) = D f ( φ t x ) ∂∂ x T t x + D f ( φ t x ) : (cid:0) T t x ⊗ T t x (cid:1) .Duhamel’s principle then yields(4.5) ∂∂ x T t x = (cid:90) t T t − s φ s x D f ( φ s x ) : ( T s x ⊗ T s x ) d s .Combining the results, we obtain(4.6) ∂ A ( φ t x ) ∂ x = D A ( φ t x ) : (cid:0) T t x ⊗ T t x (cid:1) + D A ( φ t x ) (cid:90) t T t − s φ s x D f ( φ s x ) : ( T s x ⊗ T s x ) d s ,which leads to the above claim.It is interesting that a different response formula was obtained in [35] for a stochasticperturbation of a deterministic dynamics. The derivation in [35] was also different:instead, a stochastic perturbation was used directly in the second-order response formulafor deterministic perturbations [44], which was further scaled by a factor of one-half.As we mentioned above, generally one cannot assume that the invariant measure ofthe deterministic dynamics of the form (1.1) possesses a density, since most often thecompact set on which the solution of (1.1) lives has a complicated structure. However,let us assume that there exists a smooth probability density p ( x ) >
0, such that p d x is a suitable, in an appropriate for our purposes sense, approximation for the invariantmeasure d µ . Under such a hypothetical assumption, one writes (4.1) in the form(4.7) R ( t ) = (cid:90) M ∂ A ( φ t x ) ∂ x : (cid:0) H ( x ) H T ( x ) (cid:1) p ( x ) d x .Now that µ is replaced by the density p ( x ) , one can integrate the above expression byparts, obtaining(4.8) R ( t ) = (cid:90) M A ( φ t x ) D : (cid:16) p ( x ) H ( x ) H T ( x ) (cid:17) d x .Replacing the measure averages with the time averages with help of Birkhoff’s theorem,we obtain the same formula as in (2.23):(4.9) R ( t ) = lim r → ∞ r (cid:90) r A ( x t + s ) p − ( x s ) D : (cid:16) p ( x s ) H ( x s ) H T ( x s ) (cid:17) d s .In other words, under the assumption of a differentiable approximation to the invariantstate, the time-averaged response formula for the deterministic unperturbed dynamics isidentical to the response formula for the stochastic dynamics in (2.23), where the externalstochastic perturbation is statistically independent to the unperturbed noise term.Below we will see that this approach allows to obtain a sensible approximation tothe response operator even in a situation where the unperturbed dynamics is purelydeterministic, similar to what was observed for the deterministic perturbations in [10–12, 38]. The author privately communicated that his response formula applies to a perturbation noise in theStratonovich form.
EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 13
5. T he quasi -G aussian approximation for the response operator Observe that the response formulas for the stochastic unperturbed dynamics in (2.11)and (2.23) are not computable directly, since the equilibrium density p ( x ) of (1.9) is notgenerally known explicitly. It is, theoretically, possible to compute the response in (4.1)by computing the tangent map T t x in parallel with x t (for more details, see [1–3, 7, 10–12])and using (3.20b) expressed as a time-lagged autocorrelation function over the time-series average. However, the latter option is very expensive from the computationalstandpoint, and, for a chaotic unperturbed dynamical system, it will only remain com-putationally stable for short response times.Instead, we are going to use a simplified method to compute the response, calledthe quasi-Gaussian FDT (qG-FDT) approximation [11]. The main idea of the qG-FDTapproximation is that p ( x ) in (2.11), (2.23) and (4.9) is replaced with its Gaussian ap-proximation, which has the same mean state and covariance matrix as does p ( x ) . Forthat, first observe that (2.11) and (2.23) (which is identical to (4.9)) can be written as(5.1a) R ( t ) = lim r → ∞ r (cid:90) r A ( x t + s ) D : (cid:16) G ( x s ) H T ( x s ) + H ( x s ) G T ( x s ) (cid:17) d s ++ lim r → ∞ r (cid:90) r A ( x t + s ) (cid:16) D · (cid:16) G ( x s ) H T ( x s ) + H ( x s ) G T ( x s ) (cid:17)(cid:17) · D p ( x s ) p ( x s ) d s ++ lim r → ∞ r (cid:90) r A ( x t + s ) (cid:16) G ( x s ) H T ( x s ) + H ( x s ) G T ( x s ) (cid:17) : D p ( x s ) p ( x s ) d s .(5.1b) R ( t ) = lim r → ∞ r (cid:90) r A ( x t + s ) D : (cid:16) H ( x s ) H T ( x s ) (cid:17) d s ++ lim r → ∞ r (cid:90) r A ( x t + s ) (cid:16) D · (cid:16) H ( x s ) H T ( x s ) (cid:17)(cid:17) · D p ( x s ) p ( x s ) d s ++ lim r → ∞ r (cid:90) r A ( x t + s ) (cid:16) H ( x s ) H T ( x s ) (cid:17) : D p ( x s ) p ( x s ) d s .Now, let us denote the mean state of p ( x ) as m , and its covariance matrix as C :(5.2a) m = (cid:90) R N x p ( x ) d x = lim r → ∞ r (cid:90) r x s d s ,(5.2b) C = (cid:90) R N ( x − m ) ⊗ ( x − m ) p ( x ) d x = lim r → ∞ r (cid:90) r ( x s − m ) ⊗ ( x s − m ) d s .Then, the Gaussian approximation p G ( x ) for p ( x ) is given by the explicit formula(5.3) p G ( x ) = (cid:112) ( π ) N det C exp (cid:18) − ( x − m ) · C − ( x − m ) (cid:19) ,which results in(5.4a) D p G ( x ) p G ( x ) = − C − ( x − m ) , (5.4b) D p G ( x ) p G ( x ) = C − (cid:0) ( x − m ) ⊗ ( x − m ) (cid:1) C − − C − .The approximations above are then inserted directly into (5.1a) and (5.1b), resulting inexplicit time-lagged autocorrelation functions, computed along the long-term time se-ries of solutions of the unperturbed system in (1.9). Observe that the autocorrelationsin (5.1a) and (5.1b), computed via the Gaussian approximations in (5.4), simplify some-what further if the matrices G and H are constant, as only one term out of the threeremains in such a case. 6. C omputational experiments In this section we investigate the validity of the qG-FDT response formulas in (5.1a)and (5.1b) for both the deterministic and stochastic unperturbed dynamics. For the test-bed dynamical system, we choose the well-known Lorenz 96 model.6.1.
The rescaled Lorenz 96 model.
The test system used to study the new method ofstochastic parameterization here is the rescaled Lorenz 96 system with stochastic forcing.Without the stochastic forcing, it was previously used in [4–9, 36] to study the determin-istic reduced model parameterization, as well as the average response to deterministicexternal perturbations. The rescaled Lorenz 96 system with stochastic forcing is givenby(6.1) d x i = (cid:20) x i − ( x i + − x i − ) + β ( ¯ x ( x i + − x i − ) − x i ) + F − ¯ x β (cid:21) d t + G d W it ,with 1 ≤ i ≤ N . Above, W it denotes the family of N mutually independent Wienerprocesses, indexed by i , with G being a constant stochastic forcing parameter, so that, interms of the notations in (1.9), G ( x ) is a constant multiple of the identity matrix,(6.2) G ( x ) = G I .The model has periodic boundary conditions: x i + N = x i . The parameters ¯ x and β are the statistical mean and the standard deviation, respectively, for the correspondingunrescaled Lorenz 96 model [33, 34](6.3) d x i d t = x i − ( x i + − x i − ) − x i + F ,with the same periodic boundary conditions. The rescaling above ensures that, in theabsence of the stochastic forcing (that is, G =
0) the Lorenz 96 model in (6.1) has zeromean state and unit standard deviation, and that the time scale of evolution of its so-lution is roughly the same for different values of F . In the current work, we test theresponse theory, developed above, for two values of F , F =
24 and F =
8, and two val-ues of G , G = G =
0, with the latter corresponding to the purely deterministicunperturbed dynamics. As in the original paper [34], we set N = R N , which means that theexistence of strong solutions to (6.1) is not guaranteed [19, 40]. Nonetheless, below we EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 15
Skewness Kurtosis F = G = · − F = G = · − F = G = F = G = · − able
1. The skewness and kurtosis.demonstrate via the numerical simulations that, for the chosen parameters of the system,numerical solutions exist for a long enough time to allow the reliable time-averaging forthe response computation.6.2.
Long-term statistics of the unperturbed dynamics.
The computational settings forthe numerical simulations were chosen as follows: • Forward integration scheme: Runge-Kutta 4th order for the deterministic part ofthe time step, forward Euler for the stochastic part of the time step; • Time discretization step: ∆ t = • Time averaging window: T av = • Spin-up time window (time skipped between the initial condition and the begin-ning of the time averaging window): T skip = • Initial condition: each initial state x i , 1 ≤ i ≤ N , is generated at random usingnormal distribution with zero mean and unit standard deviation.In Figure 1 we show the histograms of the probability density functions (PDFs), com-puted by the standard bin-counting, as well as the simplest time-lag autocorrelationfunctions of the solution with itself, the latter computed numerically as(6.4) C ( t ) = T av (cid:90) T skip + T av T skip x ( s ) x ( t + s ) d s ,where x ( t ) denotes one of the N variables of (6.1). Obviously, due to the translationalinvariance of (6.1), both the PDFs and correlation functions are identical across differentvariables. Observe that the PDFs look close to Gaussian, and the time autocorrelationfunctions decay rather rapidly within the first five units of time. There are two reasonswhy we need to check the decay of the time autocorrelation functions: first, we need toensure that the time averaging window T av is much longer than the decay time scale of C ( t ) for the adequate statistical averaging; and, second, we need to estimate the timescale of development of the response, since it is directly connected to the time scale ofthe autocorrelation functions according to (5.1a) and (5.1b).For more precise estimates of how close the PDFs on Figure 1 are to Gaussian, inTable 1 we show the skewness (third moment) and kurtosis (fourth moment) of thePDFs, nondimensionalized by the appropriate powers of the variance. For the purelyGaussian distribution, the skewness is zero, and the kurtosis is 3. Observe that, in thisrespect, the PDFs for the dynamical regimes with greater F and greater G are closerto the Gaussian, and thus we may expect generally better performance of the qG-FDTapproximations in (5.1a) and (5.1b) for those regimes. F igure
1. Probability density functions and time autocorrelation functions.6.3.
The external perturbations and the response function.
For testing the responsetheory developed above, we use a rather simple set-up. We set the stochastic perturba-tion matrix H entirely to zero, except for its single upper-left corner element, which isset to 1:(6.5) H ( x ) = · · ·
00 0 · · · · · · .The product εη ( t ) is set to a small constant η at zero response time:(6.6) εη ( t ) = (cid:26) t < η , t ≥ η : η = η = x , and constitutes a scalar Wiener noise η W t .As far as the choice of the response function A ( x ) is concerned, it is obvious thatmonitoring a single scalar quantity (as presented in the theory above) is not sufficient toevaluate the detailed impact of the external stochastic perturbation, even of the simplesttype we chose above, on the model. Thus, we choose to monitor the response to thestochastic perturbation of each model variable instead. More precisely, instead of mon-itoring one response function for the whole system, we monitor N of those, separatelyfor each model variable:(6.7) A i ( x ) = x i , 1 ≤ i ≤ N ,where the square of a variable x i is chosen (rather than the variable itself) because it islikely to respond more substantially to a stochastic perturbation. If we denote the setof all A i ( x ) as the vector A ( x ) , the latter can be expressed concisely as the Hadamard EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 17 product of x with itself:(6.8) A ( x ) = x ◦ x .Note that in our previous works [1–3, 7, 9–12, 38], where the deterministic perturbationswhere studied, the typical choice of A ( x ) was x itself. Here, however, we prefer (6.8)because x by itself is unlikely to respond to a stochastic perturbation in a well-articulatedfashion.The set-up above allows us to compute the response of all x i , 1 ≤ i ≤ N , separately, toa small constant stochastic forcing at the first variable, x , which is switched on at zerotime. Due to the constant nature of the stochastic perturbation, the response formulasare simplified as(6.9) δ (cid:104) A (cid:105) ( t ) = R ( t ) η , R ( t ) = (cid:90) t R ( s ) d s ,or(6.10) δ (cid:104) A (cid:105) ( t ) = R ( t ) η , R ( t ) = (cid:90) t R ( s ) d s ,depending on the type of stochastic forcing, where R ( s ) and R ( s ) are the correspond-ing vectorized leading order qG-FDT response operators from (5.1a) and (5.1b) for thevector response function in (6.8). In what follows, we display the operators R ( t ) or R ( t ) (again, depending on the type of stochastic perturbation) computed at differenttimes t = η (for comparison with R ) or η (forcomparison with R ), respectively. The ensemble simulations are performed with theensemble size of 20000 members (which were sampled from the same long term trajec-tory of the unperturbed system as was used to compute the statistics and the qG-FDTresponse), with 1000 realizations of the Wiener process carried out for each member.Each such ensemble simulation took several hours on a 16-processor Intel Xeon server,running fully in parallel. In contrast, each qG-FDT response computation took only fewminutes using a single CPU core.In addition to the plots of the response at different times, we show the relative errorsand collinearity correlations between the actually measured response, and the responsepredicted by the qG-FDT formulas in (5.1a) or (5.1b), depending on the type of noiseperturbation. The relative error is defined as the ratio of the Euclidean norm of thedifference between the qG-FDT response and the normalized measured response, overthe Euclidean norm of the qG-FDT response:(6.11) Relative error = (cid:13)(cid:13)(cid:13) δ (cid:104) A (cid:105) measured η −R (cid:13)(cid:13)(cid:13) (cid:107)R (cid:107) when perturbing the existing noise term, (cid:13)(cid:13)(cid:13)(cid:13) δ (cid:104) A (cid:105) measured η −R (cid:13)(cid:13)(cid:13)(cid:13) (cid:107)R (cid:107) when perturbing with a new noise term. The collinearity correlation is defined as the Euclidean inner product of the qG-FDTresponse with one of the measured responses, normalized by the product of their corre-sponding Euclidean norms:(6.12) Collinearity correlation = δ (cid:104) A (cid:105) measured · R(cid:107) δ (cid:104) A (cid:105) measured (cid:107)(cid:107)R(cid:107) .It is easy to see that the collinearity correlation achieves its maximum value of 1 if andonly if one response is the exact multiple of the other.6.4. Perturbing the existing noise term.
In this section we study the response of two dy-namical regimes of the stochastically forced ( G = F =
24, and is, according to Figure 1and Table 1, the closest to the Gaussian regime of all examined in Section 6.2. In Figure 2we show the response of the function (6.8) to the perturbation of the existing stochasticmatrix G = I by the perturbations described above in (6.5)–(6.6), with η set to 0.05 and0.1. With help of the periodicity of (6.1), the response variables in Figure 2, as well asall subsequent figures, are displayed so that the variable x (on which the perturbationoccurs) is at the center of the plot, with x immediately to the right, and x N to the left.Observe that the precision of the qG-FDT response prediction in this case is trulystriking – there is hardly any visual difference between the qG-FDT prediction and thedirectly measured responses for both η = η = F =
24 and G = G = I as before, by the constant deterministic forcing F is set to F =
24. Accordingto Figure 1 and Table 1, this regime is the second closest to the Gaussian, and it isclearly manifested in the difference between the qG-FDT prediction and the directlymeasured response with η = η = Perturbing with a new noise.
Here we show the results of numerical simulationswhere the external stochastic perturbation is introduced into the system via a separatenoise realization. In this situation, we can study both the fully deterministic and thestochastically forced dynamics of the rescaled Lorenz 96 model in (6.1); together withtwo different values of the forcing F , this constitutes four different combinations of pa-rameters: ( F = G = F = G = F = G = F = G = EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 19 -0.100.10.20.30.40.5 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0.5, time = 0.5qG-FDT η =0.05 η =0.1 -0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0.5, time = 1qG-FDT η =0.05 η =0.1-0.200.20.40.60.81 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0.5, time = 2qG-FDT η =0.05 η =0.1 -0.200.20.40.60.811.21.4 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0.5, time = 4qG-FDT η =0.05 η =0.1 F igure
2. Response to the existing noise perturbation, F = G = F = G = F = G = η = η = t = · − t = t = t = · − η = η = t = t = t = t = able
2. Relative errors and collinearity correlations in the response to theexisting noise perturbation, F = G = G =
0, since there is no guarantee -0.100.10.20.30.40.5 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0.5, time = 0.5qG-FDT η =0.05 η =0.1 -0.100.10.20.30.40.50.60.7 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0.5, time = 1qG-FDT η =0.05 η =0.1-0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0.5, time = 2qG-FDT η =0.05 η =0.1 -0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0.5, time = 4qG-FDT η =0.05 η =0.1 F igure
3. Response to the existing noise perturbation, F = G = F = G = F = G = η = η = t = t = t = t = η = η = t = t = t = t = able
3. Relative errors and collinearity correlations in the response to theexisting noise perturbation, F = G = µ of the unperturbed dynamics is even continuous with re-spect to the Lebesgue measure, let alone possesses a Gaussian density. However, practiceshowed in the past with deterministic perturbations [9–12, 38] that the quasi-Gaussian EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 21
Rel. Errors, F = G = F = G = η = η = t = t = t = t = η = η = t = t = t = t = able
4. Relative errors and collinearity correlations in the response to theperturbation with a new noise, F = G = F = G = F = G = η = η = t = t = t = t = η = η = t = t = t = t = able
5. Relative errors and collinearity correlations in the response to theperturbation with a new noise, F = G = G = G = I to each of the 20000 ensemble members, while modelingthe stochastic perturbation via another 1000 independent realizations of the Wiener noisefor each ensemble simulation as before. This is done to retain the same computationalexpense as for the other studied cases, which, of course, leads to statistical undersam-pling of the expectation over the noise realizations. Indeed, in the case of stochasticperturbations independent from the unperturbed noise, the expectation must be com-puted over the comparable number of realizations separately for the unperturbed noiseand for the stochastic perturbation to retain comparable averaging fidelity. Still, we findthat even in this simplified setup the qG-FDT approximation shows good agreementwith the measured response, at least for sufficiently short response times.In Figure 4 and Table 4 we show the qG-FDT prediction together with the directlymeasured response with perturbations η = η = F =
24 and G = t ≤
2. However, what we can also ob-serve is the large discrepancy between the two directly measured responses for differentperturbation magnitudes, which develops at t =
4. Note that no such discrepancy was -0.100.10.20.30.40.5 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0, time = 0.5qG-FDT η =0.05 η =0.1 -0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0, time = 1qG-FDT η =0.05 η =0.1-0.200.20.40.60.81 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0, time = 2qG-FDT η =0.05 η =0.1 -1-0.500.511.5 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0, time = 4qG-FDT η =0.05 η =0.1 F igure
4. Response to the perturbation with a new noise, F = G = η = η = F =
8, and allother parameters set as above. Here, observe that the quality of the qG-FDT predictiondeteriorates even further, as the relative errors increase to 36-53%, and the collinearitycorrelations drop to 93-87% for the response times t ≤
2. This is to be expected, however,since this regime is the farthest from the Gaussian, according to Figure 1 and Table 1.The discrepancy between the directly measured responses for perturbations of differentmagnitudes, which could be attributed to the developing structural instability of theattractor, also manifests itself at time t =
4, as for the previously studied deterministicregime with F = EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 23 -0.2-0.100.10.20.30.4 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0, time = 0.5qG-FDT η =0.05 η =0.1 -0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0, time = 1qG-FDT η =0.05 η =0.1-0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0, time = 2qG-FDT η =0.05 η =0.1 -0.4-0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0, time = 4qG-FDT η =0.05 η =0.1 F igure
5. Response to the perturbation with a new noise, F = G = F = G = F = G = η = η = t = · − t = · − t = t = η = η = t = t = t = t = able
6. Relative errors and collinearity correlations in the response to theperturbation with a new noise, F = G = η = η = F =
24 and G = -0.100.10.20.30.40.5 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0.5, time = 0.5qG-FDT η =0.05 η =0.1 -0.100.10.20.30.40.50.60.7 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0.5, time = 1qG-FDT η =0.05 η =0.1-0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0.5, time = 2qG-FDT η =0.05 η =0.1 -2-1012345 26 31 36 1 6 11 16 R e s pon s e VariableF=24, G=0.5, time = 4qG-FDT η =0.05 η =0.1 F igure
6. Response to the perturbation with a new noise, F = G = F = G = F = G = η = η = t = t = t = t = η = η = t = t = t = t = able
7. Relative errors and collinearity correlations in the response to theperturbation with a new noise, F = G = t = t =
1, with relative errors about 8-16%, and collinearity correlations about
EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 25 -0.100.10.20.30.4 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0.5, time = 0.5qG-FDT η =0.05 η =0.1 -0.100.10.20.30.40.50.60.7 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0.5, time = 1qG-FDT η =0.05 η =0.1-0.200.20.40.60.8 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0.5, time = 2qG-FDT η =0.05 η =0.1 -1.5-1-0.500.51 26 31 36 1 6 11 16 R e s pon s e VariableF=8, G=0.5, time = 4qG-FDT η =0.05 η =0.1 F igure
7. Response to the perturbation with a new noise, F = G = t =
2. This is likely the manifestation of statisticalundersampling mentioned above, since the size of the statistical ensemble is unchangedeven though an additional independent stochastic forcing is introduced into the dynam-ics.In Figure 7 and Table 7 we show the results for the regime with F = G = t = t = F = G = the errors for longer times, t = t =
4, are smaller than those for the regime with F =
24, again, likely due to the fact that the statistical undersampling does not manifestas strongly in a less chaotic regime. The collinearity correlations appear to follow thesame trend as the relative errors. 7. S ummary
In this work we develop a fluctuation-response theory and test a computational al-gorithm for the leading order response of chaotic and stochastic dynamical systems tostochastic perturbations. The key property of this approach is that it allows to estimatethe average response to an external stochastic perturbation from a certain combinationof the time-lagged averages of the unperturbed system. For dynamical systems, whichare already stochastic, we consider two cases: first, where the existing stochastic term isperturbed; and, second, where a new stochastic perturbation is introduced, which cor-respondingly leads to different leading order average response formulas. We also showthat, under appropriate assumptions, the resulting formulas for leading order responseto a stochastic perturbation for the deterministic and stochastic unperturbed dynamicsare equivalent. For the practical computation of the leading order response approxima-tion, we derive the approximate quasi-Gaussian response formulas, where the probabil-ity density of the unperturbed statistical state is assumed to be Gaussian. We numericallyinvestigate the validity of the quasi-Gaussian response formulas for stochastic perturba-tions of both deterministic and stochastically forced Lorenz 96 system. We find that thequasi-Gaussian response formulas appear to more effective for the regimes where theunperturbed dynamics is already stochastic. Additionally, for the stochastic perturba-tions of the deterministic Lorenz 96 model, we observe what seems to be a manifestationof structural instability of the system’s attractor under a stochastic perturbation.
Acknowledgment.
The work was supported by the Office of Naval Research grantN00014-15-1-2036. R eferences [1] R.V. Abramov. Short-time linear response with reduced-rank tangent map.
Chin. Ann. Math. ,30B(5):447–462, 2009.[2] R.V. Abramov. Approximate linear response for slow variables of deterministic or stochastic dynamicswith time scale separation.
J. Comput. Phys. , 229(20):7739–7746, 2010.[3] R.V. Abramov. Improved linear response for stochastically driven systems.
Front. Math. China ,7(2):199–216, 2012.[4] R.V. Abramov. A simple linear response closure approximation for slow dynamics of a multiscalesystem with linear coupling.
Multiscale Model. Simul. , 10(1):28–47, 2012.[5] R.V. Abramov. Suppression of chaos at slow variables by rapidly mixing fast dynamics through linearenergy-preserving coupling.
Commun. Math. Sci. , 10(2):595–624, 2012.[6] R.V. Abramov. A simple closure approximation for slow dynamics of a multiscale system: Nonlinearand multiplicative coupling.
Multiscale Model. Simul. , 11(1):134–151, 2013.[7] R.V. Abramov. Linear response of the Lyapunov exponent to a small constant perturbation.
Comm.Math. Sci , 14(4):1155–1167, 2016.[8] R.V. Abramov. A simple stochastic parameterization for reduced models of multiscale dynamics.
Fluids , 1(1), 2016.
EADING ORDER RESPONSE OF STATISTICAL AVERAGES TO STOCHASTIC PERTURBATIONS 27 [9] R.V. Abramov and M. Kjerland. The response of reduced models of multiscale dynamics to smallexternal perturbations.
Commun. Math. Sci. , 14(3):831–855, 2016.[10] R.V. Abramov and A.J. Majda. Blended response algorithms for linear fluctuation-dissipation forcomplex nonlinear dynamical systems.
Nonlinearity , 20:2793–2821, 2007.[11] R.V. Abramov and A.J. Majda. New approximations and tests of linear fluctuation-response forchaotic nonlinear forced-dissipative dynamical systems.
J. Nonlin. Sci. , 18(3):303–341, 2008.[12] R.V. Abramov and A.J. Majda. New algorithms for low frequency climate response.
J. Atmos. Sci. ,66:286–309, 2009.[13] T. Bell. Climate sensitivity from fluctuation dissipation: Some simple model tests.
J. Atmos. Sci. ,37(8):1700–1708, 1980.[14] G.D. Birkhoff. Proof of the ergodic theorem.
Proc. Natl. Acad. Sci. , 17(12):656–660, 1931.[15] G. Carnevale, M. Falcioni, S. Isola, R. Purini, and A. Vulpiani. Fluctuation-response in systems withchaotic behavior.
Phys. Fluids A , 3(9):2247–2254, 1991.[16] B. Cohen and G. Craig. The response time of a convective cloud ensemble to a change in forcing.
Quart. J. Roy. Met. Soc. , 130(598):933–944, 2004.[17] J.-P. Eckmann and D. Ruelle. Ergodic theory of chaos and strange attractors.
Rev. Mod. Phys. , 57(3):617–656, 1985.[18] D. Evans and G. Morriss.
Statistical Mechanics of Nonequilibrium Liquids . Academic Press, New York,1990.[19] I.I. Gikhman and A.V. Skorokhod.
Introduction to the Theory of Random Processes . Courier Dover Publi-cations, 1969.[20] I.I. Gikhman and A.V. Skorokhod.
The Theory of Stochastic Processes I . Classics in Mathematics.Springer, 2004.[21] A. Gritsun. Fluctuation-dissipation theorem on attractors of atmospheric models.
Russ. J. Numer. Math.Modeling , 16(2):115–133, 2001.[22] A. Gritsun and G. Branstator. Climate response using a three-dimensional operator based on thefluctuation-dissipation theorem.
J. Atmos. Sci. , 64:2558–2575, 2007.[23] A. Gritsun, G. Branstator, and V. Dymnikov. Construction of the linear response operator of an atmo-spheric general circulation model to small external forcing.
Num. Anal. Math. Modeling , 17:399–416,2002.[24] A. Gritsun, G. Branstator, and A.J. Majda. Climate response of linear and quadratic functionals usingthe fluctuation dissipation theorem.
J. Atmos. Sci. , 65:2824–2841, 2008.[25] A. Gritsun and V. Dymnikov. Barotropic atmosphere response to small external actions. theory andnumerical experiments.
Atmos. Ocean Phys. , 35(5):511–525, 1999.[26] K. It ˆo. Stochastic integral.
Proc. Imperial Acad. Tokyo , 20:519–524, 1944.[27] K. It ˆo. On stochastic differential equations.
Mem. Amer. Math. Soc. , 4:1–51, 1951.[28] R. Kubo. Statistical mechanical theory of irreversible processes I: General theory and simple applica-tions to magnetic and conduction problems.
J. Phys. Soc. Jpn. , 12:507–586, 1957.[29] R. Kubo. The fluctuation-dissipation theorem.
Rep. Prog. Phys. , 29:255–284, 1966.[30] R. Kubo, M. Toda, and N. Hashitsume.
Statistical Physics II: Nonequilibrium Statistical Mechanics .Springer-Verlag, New York, 1985.[31] H. Kunita.
Stochastic flows and stochastic differential equations . Cambridge University Press, 1997.[32] C. Leith. Climate response and fluctuation-dissipation.
J. Atmos. Sci. , 32:2022–2025, 1975.[33] E. Lorenz. Predictability: A problem partly solved. In
Proceedings of the Seminar on Predictability , Shin-field Park, Reading, England, 1996. ECMWF.[34] E. Lorenz and K. Emanuel. Optimal sites for supplementary weather observations.
J. Atmos. Sci. ,55:399–414, 1998.[35] V. Lucarini. Stochastic perturbations to dynamical systems: A response theory approach.
J. Stat. Phys ,146:774–786, 2012.[36] V. Lucarini and S. Sarno. A statistical mechanical approach for the computation of the climatic re-sponse to general forcings.
Nonlin. Processes Geophys. , 18:7–28, 2011. [37] A.J. Majda, R.V. Abramov, and B. Gershgorin. High skill in low frequency climate response throughfluctuation dissipation theorems despite structural instability.
Proc. Natl. Acad. Sci. , 107(2):581–586,2010.[38] A.J. Majda, R.V. Abramov, and M.J. Grote.
Information Theory and Stochastics for Multiscale Nonlin-ear Systems , volume 25 of
CRM Monograph Series of Centre de Recherches Math´ematiques, Universit´e deMontr´eal . American Mathematical Society, 2005. ISBN 0-8218-3843-1.[39] B. Øksendal.
Stochastic Differential Equations: An Introduction with Applications . Universitext. Springer,6th edition, 2010.[40] G. Pavliotis.
Stochastic Processes and Applications , volume 60 of
Texts in Applied Mathematics . Springer,2014.[41] H. Risken.
The Fokker-Planck Equation . Springer-Verlag, New York, 2nd edition, 1989.[42] D. Ruelle. Differentiation of SRB states.
Comm. Math. Phys. , 187:227–241, 1997.[43] D. Ruelle. General linear response formula in statistical mechanics, and the fluctuation-dissipationtheorem far from equilibrium.
Phys. Lett. A , 245:220–224, 1998.[44] D. Ruelle. Nonequilibrium statistical mechanics near equilibrium: Computing higher order terms.
Nonlinearity , 11:5–18, 1998.[45] L.-S. Young. What are SRB measures, and which dynamical systems have them?