Universal Upper Estimate for Prediction Errors under Moderate Model Uncertainty
UUniversal Upper Estimate for Prediction Errors under Moderate ModelUncertainty
Bálint Kaszás and George Haller Institute for Mechanical Systems, ETH Zürich, Leonhardstrasse 21, 8092 Zürich, Switzerland
We derive universal upper estimates for model-prediction error under moderate but otherwise unknown model uncer-tainty. Our estimates give upper bounds on the leading order trajectory-uncertainty arising along model trajectories,solely as functions of the invariants of the known Cauchy-Green strain tensor of the model. Our bounds turn out tobe optimal, which means that they cannot be improved for general systems. The quantity relating the leading-ordertrajectory-uncertainty to the model uncertainty is the Model Sensitivity, which we find to be a useful tool for a quickglobal assessment of the impact of modeling uncertainties in various domains of the phase space. Examining the expec-tation that Finite-Time Lyapunov Exponents capture sensitivity to modeling errors, we show that this does not generallyfollow. However, we find indication that certain ridges of FTLE may persist in the MS field.
We present a method of sensitivity analysis for generaldynamical systems, subjected to deterministic or stochas-tic modeling uncertainty. Using the properties of the un-perturbed dynamics, we derive a universal bound for theleading-order prediction error. This bound motivates thedefinition of the Model Sensitivity, a scalar quantity, de-pending on the initial condition and time. We demon-strate, using nonlinear numerical models that the ModelSensitivity provides both a global view over the phasespace of the dynamical system and in some situations, alocalized, time-dependent predictor of uncertainties alongtrajectories. We find that the phase-space structure ofthe Model Sensitivity (MS) is related, but not identical tothat of the Finite-Time Lyapunov Exponents (FTLE). Weformulate conditions under which robust features of theFTLE field are expected to also be seen in the MS field.
I. INTRODUCTION
One of the challenges in modeling real-world phenomena isuncertainties that enter the model formulation. Depending onthe context, these can arise as a result of incomplete or noisydata, uncertainty in the mathematical model, or even the errorintroduced by numerical algorithms. Here we seek to boundthe impact of these uncertainties on specific model trajecto-ries utilizing minimal information on the modeling errors butsubstantial information on the internal dynamics of the knownmodel.A range of approaches exist to assess the impact of modeluncertainty. One such approach, response theory, originatesfrom statistical physics, where a central question was an equi-librium system’s response to infinitesimal perturbations. Un-der a (possibly time-dependent) perturbation to a Hamiltoniansystem, Kubo’s formula establishes a link between the ex-pected value of the linear-order response and certain quanti-ties of the unperturbed system. This linear response theorywas generalized to systems with uniform hyperbolicity. Withthis assumption, Ruelle’s work put the theory on a rigorousfoundation, providing formulas for the asymptotic expansionof the invariant measure of the perturbed system. Numericalevidence shows that linear response can be observed even in systems that are not strictly uniformly hyperbolic. In partic-ular, in the field of climate science, response theory has beensuccessfully used to assess the various possible scenarios ofanthropogenic climate change .For general dynamical systems, an additional source ofuncertainty is also present: sensitivity to initial conditions.This means that a small error in the system’s initial condi-tion grows exponentially, governed by the largest Lyapunovexponent . A common illustration of this phenomenon isweather prediction, in which long term predictions are impos-sible due to the exponential error-growth. Considerable efforthas gone into assessing these difficulties, for example, by as-similating real-world observations into model prediction ,or by using ensemble methods to obtain a statisticalcharacterization.Another important question is the sensitivity of model pre-dictions to slight changes in the model parameters . Thissensitivity is often characterized by the derivatives of an ob-servable (a function of the model variables) with respect tothose parameters . While response theory studies the sys-tem’s behavior under external influence, sensitivity analysisseeks to assess prediction errors under modeling uncertainty.To reduce the computational cost, the sensitivity is oftencomputed from the linearized dynamics (tangent method) along a reference trajectory.Usually, the observed quantity is an infinite-time averagecomputed along trajectories. Direct calculations need to uti-lize sufficiently long Monte Carlo simulations of the fullmodel and finite difference approximations for the deriva-tive. In this case, the tangent method generally results inasymptotically unbounded sensitivities, which do not matchthe bounded ones computed directly . As noted in Ref. 22,the issue comes from exchanging two limits: the sensitivity ofan infinite-time average is the derivative of the infinite-timeaverage, while the tangent method calculates the infinite-timeaverage of a derivative.It has been suspected that similarly to sensitivity with re-spect to initial conditions, sensitivity with respect to param-eters is also governed by the largest Lyapunov exponent ofthe underlying trajectory. This is supported by numericalresults , but can also be intuitively understood: the differ-ential equation that describes the growth of perturbations toinitial conditions has the same homogeneous part as the one a r X i v : . [ n li n . C D ] J u l describing error growth due to parameter changes .To circumvent this problem of unbounded averages, theensemble method calculates the sensitivity over shorter timeintervals for several randomly selected trajectories usingthe linearized dynamics. Then, the true sensitivity of theinfinite-time average can be approximated by the ensembleaverage .For ergodic systems, the least-squares shadowingmethod offers an alternative calculation of the trueparameter sensitivity of an infinite-time average. Instead ofsolving the tangent equation, the method looks for a nearbyshadowing trajectory which has a uniformly bounded distancefrom the reference trajectory. Practically, this means that anonlinear optimization problem has to be solved. Solvingthe linearized version of this problem, it is possible to obtainmeaningful sensitivities , even for chaotic systems . Themethod was also implemented in turbulent fluid dynamicalsimulations . Further improvements on calculating sensitiv-ities for chaotic systems take advantage of unstable periodicorbits .In contrast to the methods mentioned above, we focus hereon finite-time predictions and their uncertainties, as opposedto infinite-time averages. This is motivated by the fact thatcertain models may not be defined for infinite times, or theinfinite time averages may not be accurate representations ofthe system .We derive universal bounds on the sensitivity of model pre-dictions under small modeling errors. Our estimates only as-sume the knowledge of a general bound on the model errors,yet yield trajectory-specific bounds for the model-predictionerrors. These bounds provide a granular assessment of theimpact of modeling errors, depending only on the known lo-cal dynamics of the phase space in the absence of model un-certainties. We relate the arising model sensitivities to FiniteTime Lyapunov Exponents (FTLE) and their ridges , andhence to Lagrangian Coherent Structures (LCS) , which areorganizing structures in the idealized model’s phase space. Wefind that in general, the ridges FTLE do not signal the pres-ence of a ridge in the model sensitivity. However, we for-mulate a plausible condition under which a correspondence isexpected.We also extend the analysis to cases when both determin-istic and stochastic uncertainties are present. We show thatassuming multiplicative Gaussian noise, the expectation valueof the observation error can be bounded by an asymptotic for-mula, analogous to the purely deterministic case. All theseestimates even turn out to be optimal: we give examples inwhich the inequalities become equalities. In addition, throughnumerical examples of models that represent various levels ofcomplexity, we show that the bounds developed for the obser-vation error hold for surprisingly large modeling uncertaintiestoo. II. SET-UP
We first consider a parametrized family of deterministic dif-ferential equations˙ x = f ( x , t )+ ε g ( x , t , ε ) x ∈ U ⊂ R n , t ∈ [ t , t ] , ≤ ε (cid:28) f and g are assumed to be smooth functions oftheir arguments. Trajectories of this equation are of the form x ( t ; t , x , ε ) , which are as smooth in their arguments as f is.We can think of ε g ( x , t ; ε ) as a family of perturbations repre-senting errors to a known model system˙ x = f ( x , t ) , (II.2)our “best understanding” of the given problem. The perturba-tions of the form ε g ( x , t ; ε ) represent the modeling uncertaintyof the underlying problem, such as a systematic bias with spa-tial and temporal dependence. We assume that this term isbounded in norm.We are interested in how trajectories change under changesin the parameter ε . While the exact nature of the family ε g ( x , t ; ε ) is generally unknown for ε >
0, we still seek toassess the leading-order uncertainty of trajectories in case anoverall bound on ε | g ( x , t ; ε ) | is available. We call this leading-order uncertainty the model-sensitivity of the trajectory withrespect to the parameter ε .We will show that even for completely general systems,there exists a bound on the leading-order uncertainty, whichcan, in practice, be even used to bound the proper trajectoryuncertainty.Next, we will assume stochastic model uncertainty byadding a white-noise-driven stochastic process in the pertur-bation to the known vector field f . This translates into astochastic differential equation of the form dx t = f ( x t , t ) dt + ε g ( x t , t , ε ) dt + εσ ( x t , t ) d W t , (II.3) x ∈ U ⊂ R n , t ∈ [ t , t ] , ≤ ε (cid:28) . The SDE is understood in terms of the Itô interpretation,where W t is an n -dimensional Wiener-process and σ ( x t , t ) isthe covariance matrix. The coefficient functions in (II.3) willbe assumed to satisfy additional assumptions that guaranteethe existence of a strong solution to the equation. These typesof stochastic perturbations represent either random errors inthe model or unresolved effects. III. DETERMINISTIC MODEL SENSITIVITY
Traditionally, model sensitivity refers to the change in anobservable under changes in the model equations, usuallythrough parameters. Here, we take the observable to be simplythe model trajectory itself.
A. Influence of deterministic uncertainty
Let us first assume that (II.1) holds, generating a flow-mapF tt : U → U , x (cid:55)→ x ( t , t , x ; ε ) . (III.1)To assess the effect of a slight change in ε , we consider thenorm of the difference between the idealized model trajectory x ( t ) : = x ( t , t , x ; 0 ) (III.2)and the real one (with ε (cid:54) = x ε ( t ) : = x ( t , t , x ; ε ) (III.3)starting from the same initial condition. This trajectory un-certainty is given by | x ( t , t , x ; 0 ) − x ( t , t , x ; ε ) | = | x ( t ) − x ε ( t ) | . (III.4)By classic results on ordinary differential equations, theflow map F tt is as smooth in the parameter ε as is the vec-tor field f + ε g , and hence can also be Taylor-expanded in ε .This gives the leading-order trajectory uncertainty as ε (cid:12)(cid:12)(cid:12)(cid:12) ∂ x ε ( t ) ∂ ε (cid:12)(cid:12)(cid:12)(cid:12) ε = = ε | η ( t , t , x ) | . (III.5)The vector η , which is the derivative of the flow-mapwith respect to ε , obeys the (inhomogeneous) equation ofvariations , also called the tangent model of II.1:˙ η = ∇ f (cid:0) x ( t ) (cid:1) η + g (cid:0) x ( t ) , t ; 0 (cid:1) , η ( t ; t , x ) = . (III.6)The solution of this initial value problem is η ( t ; t , x ) = (cid:90) tt φ ts (cid:0) x ( s ) (cid:1) g (cid:0) x ( s ) , s ; 0 (cid:1) ds , (III.7)where the deformation gradient, φ tt ( x ) = ∇ F tt ( x ) , is thenormalized fundamental matrix solution of the equation ofvariations ˙ η = ∇ f (cid:0) x ( t ) , t (cid:1) η , (III.8)i.e., the homogeneous part of the linear system of ordinarydifferential equations (III.6).Therefore, the leading-order change to a trajectory x ( t ) due to changes in the model is ε | η ( t ; t , x ) | = ε (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) tt φ ts (cid:0) x ( s ) (cid:1) g (cid:0) x ( s ) , s ; 0 (cid:1) ds (cid:12)(cid:12)(cid:12)(cid:12) . (III.9) This quantity can be bounded from above as ε | η ( t ; t , x ) | ≤ ε (cid:90) tt (cid:12)(cid:12) φ ts (cid:0) x ( s ) (cid:1) g (cid:0) x ( s ) , s ; 0 (cid:1)(cid:12)(cid:12) ds ≤ (cid:90) tt (cid:13)(cid:13) φ ts (cid:0) x ( s ) (cid:1)(cid:13)(cid:13) (cid:12)(cid:12) ε g (cid:0) x ( s ) , s ; 0 (cid:1)(cid:12)(cid:12) ds ≤ ε (cid:90) tt (cid:13)(cid:13) φ ts (cid:0) x ( s ) (cid:1)(cid:13)(cid:13) ds (cid:13)(cid:13) g (cid:0) x ( s ) , s ; 0 (cid:1)(cid:13)(cid:13) ∞ ≤ ε (cid:90) tt (cid:113) Λ ts ( x ( s )) ds (cid:13)(cid:13) g (cid:0) x ( s ) , s ; 0 (cid:1)(cid:13)(cid:13) ∞ , (III.10)where || · || ∞ refers to the supremum norm and Λ ts (cid:0) x ( s ) (cid:1) de-notes the largest eigenvalue of the (right) Cauchy–Green straintensor C ts (cid:0) x ( s ) (cid:1) = (cid:2) φ ts (cid:0) x ( s ) (cid:1)(cid:3) T φ ts (cid:0) x ( s ) (cid:1) . In other words, (cid:112) Λ ts ( x ( s )) is the largest singular value of φ ts (cid:0) x ( s ) (cid:1) .Let ∆ ∞ ( x , t ) : = ε (cid:13)(cid:13) g (cid:0) x ( · ) , · ; 0 (cid:1)(cid:13)(cid:13) ∞ = (III.11) = ε max s ∈ [ t , t ] (cid:12)(cid:12) g (cid:0) x ( s ) , s ; 0 (cid:1)(cid:12)(cid:12) denote the maximal leading-order model uncertainty alongthe trajectory x ( t ) of the idealized model (II.2). With this no-tation, let us define the leading-order trajectory uncertainty atany time instant t ∈ [ t , t ] as δ ( x , t ) : = ε | η ( t ; t , x ) | ≤ (cid:90) tt (cid:113) Λ ts ( x ( s )) ds ∆ ∞ ( x , t ) . (III.12)For any finite k ∈ N + , we also define the corresponding timeaveraged leading-order trajectory uncertainty as the temporal L k norm of δ ( x , t ) : δ k ( x ) : = (cid:107) δ ( x , t ) (cid:107) L k = ε k (cid:115) (cid:90) t t [ η ( t ; t , x )] k dt . (III.13)To obtain a uniform bound for δ k ( x ) over the time interval [ t , t ] , we can simply let k → ∞ and find that δ ∞ ( x ) : = (cid:107) δ ( x , t ) (cid:107) ∞ = ε max t ∈ [ t , t ] | η ( t , t , x ; 0 ) | . (III.14)Similarly, for the maximal leading-order model uncertainty,we can set ∆ ∞ ( x ) : = ε max s ∈ [ t , t ] ∆ ∞ ( x , s ) = ε max s ∈ [ t , t ] (cid:12)(cid:12) g (cid:0) x ( s ) , s ; 0 (cid:1)(cid:12)(cid:12) . (III.15)Then, by Eq. (III.10), the trajectory uncertainty δ k ( x ) can beestimated from above as δ k ( x ) = ε k (cid:115) (cid:90) t t [ η ( t ; t , x )] k dt ≤ ∆ ∞ ( x , t ) (cid:13)(cid:13)(cid:13)(cid:13) (cid:90) tt (cid:113) Λ ts ( x ( s )) ds (cid:13)(cid:13)(cid:13)(cid:13) L k . (III.16)Taking the supremum norm of both sides gives δ ∞ ( x ) ≤ ∆ ∞ ( x , t ) max t ∈ [ t , t ] (cid:90) tt (cid:113) Λ ts ( x ( s )) ds . (III.17)Note that, the upper bound (cid:82) tt (cid:112) Λ ts ( x ( s )) ds ∆ ∞ ( x , t ) isgenerally not a monotone function of t . This implies that inorder to evaluate the right hand side of (III.17), one needs tocompute the integral involved for all time instants in [ t , t ] .In conclusion, the estimate (III.12) shows that the leading-order trajectory uncertainty under modeling errors can be es-timated from above by a measure of sensitivity with respect toinitial conditions within the model. Remark . Recall that sensitivity with respect to initial condi-tions over a time interval [ s , t ] is typically characterized by thefinite-time Lyapunov exponent (or FTLE) , given byFTLE ts (cid:0) x ( s ) (cid:1) = t − s log (cid:113) Λ ts ( x ( s )) . (III.18)With this in mind, our uncertainty estimate can be rewrittenas a functional of the FTLE field as follows: δ ( x , t ) ≤ ∆ ∞ ( x , t ) (cid:90) tt exp (cid:2) ( t − s ) FTLE ts (cid:0) x ( s ) (cid:1)(cid:3) ds . (III.19) Remark . The calculation of the integral in (III.12) can alsobe done in backward time, which is sometimes more conve-nient. Following the results on the smallest eigenvalue of theCauchy-Green strain tensor , we note that (cid:113) Λ ts ( x ( s )) = (cid:112) λ min [ C st ( x ( t ))] (III.20)and hence (cid:90) tt (cid:113) Λ ts ( x ( s )) ds = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) t t (cid:112) λ min [ C st ( x )] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (III.21)Numerically, formula (III.21) requires the evaluation ofthe integral of the square-root of the largest eigenvalue ofbackward-time Cauchy–Green strain tensor C st (cid:0) x ( t ) (cid:1) , com-puted over [ t , s ] , with s decreasing from t to t . This canbe computed by finite-differencing along backward-time tra-jectories starting from a regular grid at time t , back to time t . The advantage of this approach is that the Cauchy-Greenstrain tensor is always calculated for the same initial point dur-ing integration. However, this point is the time- t position ofthe idealized model trajectory. To obtain the bound as a func-tion of the time- t position, one needs to map the values backfrom time t to time t with the idealized flow-map.With the above estimates, we can now bound the leading-order trajectory uncertainty of the dynamical system. Mostmethods currently available for calculating sensitivity mea-sures need additional assumptions, such as the existence of aninvariant measure , ergodicity or a specific form of themodeling errors (to run direct simulations). While these oftengive precise predictions on the value of the sensitivity, they are not applicable to typical dynamical systems. In contrast, theinequality (III.12) holds for all dynamical systems of the form(II.1). In addition, we can also use it to formulate a bound onthe proper (not only leading order) uncertainty in the dynami-cal system’s trajectories. Theorem 1 (Universal bound on trajectory uncertainty) . Con-sider the dynamical system defined over a finite time interval [ t , t ] , and on a compact domain U ⊂ R n , by (II.1) . Denoteby x ( t ) the idealized model’s solution (III.2) , starting from x at t . Similarly, let x ε ( t ) be the true solution (III.3) , belongingto an arbitrary ε (cid:54) = , starting from the same initial condition.Then, for any δ > small enough, there exists ε > , suchthat for ε < ε the following inequality holds for all t ∈ [ t , t ] and x ∈ U | x ε ( t ) − x ( t ) | ≤ (cid:18) (cid:90) tt (cid:113) Λ ts ( x ( s )) ds + δ (cid:19) ∆ ∞ ( x , t ) . (III.22) Proof.
See Appendix VIII A.Our Theorem 1 provides a bound for small values of ε ,that is computable numerically and holds for any time instantin the time interval and any initial condition in the domain.At first, the dependence on a finite δ may seem problematic.However, we note that for small enough ε , the size- δ correc-tion can be made arbitrarily small. Our numerical findingsindicate that (III.22) tends to be satisfied even for δ = eitherrequire knowledge of the perturbed trajectory itself or intro-duce Gronwall-type estimates that vastly overestimate the er-ror, due to their universality in space and time.For example, assume that in system (II.1), f satisfies theLipschitz condition with Lipschitz constant L and the per-turbation ε g ( x , t ) is uniformly bounded by a constant M = max ∆ ∞ . We then obtain | x ε ( t ) − x ( t ) | ≤ ML ( e L ( t − t ) − ) fromthe classic Gronwall-lemma . This is a rigorous but highlyconservative upper bound on the trajectory uncertainty, asseen from a direct comparison with (III.22).To illustrate the difference between the two estimates, con-sider the classic damped-forced Duffing-oscillator˙ x = y , (III.23)˙ y = x − x − δ y + A cos t , with δ = .
15 and A = .
3. For these parameter values,the system has a global attractor , contained in the region U = [ − . , . ] × [ − . , . ] . Using this fact, we choose theglobal Lipschitz-constant L = ε g ( x , t ) = ( , ε sin ( ω p t )) , representing a high-frequencydeterministic perturbation to system (III.23). For this choiceof perturbation, we can select the uniform bound M = ε forthe model error. . . . . . . . . . . . . . . E rr o r Trajectory uncertaintyGronwall-estimateLeading-ordertrajectory uncertainty estimate (III.12)
FIG. 1. Comparison of uncertainty estimates, applied to systemIII.23. The dashed-dotted curve shows the phase-space distance be-tween the idealized model solution and the real one, started from thesame initial condition ( . , . ) , with ε = .
01 and ω p =
10. Thered curve is the bound obtained from the leading-order bound in in-equality (III.12). The blue curve is the Gronwall-type, rigorous upperbound, defined by ML ( e Lt − ) , with M = ε and L = As seen in Fig 1, both the Gronwall-type estimate and theleading-order bound of (III.12) substantially overestimate theactual distance between the true and the idealized model tra-jectories. However, while the Gronwall-estimate suggests anoverall exponential increase for all trajectories starting in U ,our leading-order bound (IV.12) depends on the unperturbedtrajectory, providing a tighter estimate on the trajectory uncer-tainty. IV. STOCHASTIC MODEL SENSITIVITY
It is often reasonable to assume a stochastic model error asone of the sources of uncertainty in the model. In that case,trajectories obey the following stochastic differential equation(SDE) ˙ x = f ( x , t ) + ε g ( x , t , ε ) + εσ ( x , t ) ξ ( t ) . (IV.1)The uncertainty comes from a white-noise process ξ and f , g , σ are smooth functions. Equation (IV.1) can be inter-preted in the Itô-sense as dx t = f ( x t , t ) dt + ε g ( x t , t , ε ) dt + εσ ( x t , t ) dW t , (IV.2)on a probability space ( Ω , F , P ) , with W t being an n -dimensional Wiener process, f : R n × [ t , t ] × R → R n is thedeterministic part of the SDE, and σ ( · , · ) : R n × [ t , t ] → R n × n is the covariance matrix of the noise. Both f and σ are as-sumed to be measurable, smooth functions of their arguments.Our goal is to characterize the leading-order deviation ofthe solution process x t of (IV.2) from the solution of the ide-alized model ( ε = ω ∈ Ω of thenoise. To achieve such a characterization, we develop an up-per estimate similar to (III.17). We first state the necessaryand sufficient conditions for the existence of a solution pro-cess x t , derive the SDE governing the leading-order trajectoryuncertainty (a stochastic analog to the equation of variations),and give bounds on the expected value of the norm of its so-lutions.Assume that there exist constants C , D >
0, such that for all x , y ∈ R n , t ∈ [ t , t ] and small enough ε >
0, we have | f ( x , t ) + ε g ( x , t , ε ) | + | εσ ( x , t ) | ≤ C ( + | x | ) , | f ( x , t ) + ε g ( x , t , ε ) − f ( y , t ) − ε g ( y , t , ε ) | + | εσ ( x , t ) − εσ ( y , t ) | ≤ D | x − y | . (IV.3)Then, Equation (IV.2) along with the deterministic ini-tial condition x t = t = x has a unique solution x t which isadapted to the filtration generated by W s for s ≤ t . In addi-tion, E (cid:16) (cid:82) t t | x t | dt (cid:17) < ∞ holds and the sample paths of thesolution x t ( ω ) are continuous .The following theorem pro-vides an analogue of the equation of variations (III.8) in thestochastic setting. Theorem 2 (Small noise expansion) . Assume that the coeffi-cients in (IV.2) have bounded and measurable partial deriva-tives up to second order. Then, there exists ¯ ε > , such thatfor ε < ¯ ε the solution x ε t can be written asx ε t = x t + εη t + ε R ( t , ε ) , (IV.4) with the same notation as we had in (III.6), but now with η t denoting a stochastic process. The remainder term, R ( t , ε ) ,is bounded in the mean-squared sense, i.e., there exists K > ,such that sup t ∈ [ t , t ] (cid:104) E | R ( t , ε ) | (cid:105) ≤ K . (IV.5) The coefficients x t and η t satisfy the system of stochasticdifferential equationsdx t = f (cid:0) x t , t (cid:1) dt , x t = t = x , (IV.6) d η t = ∇ f (cid:0) x t , t (cid:1) η t dt + g (cid:0) x t , t ; 0 (cid:1) dt + σ (cid:0) x t , t (cid:1) dW t , η t = t = . (IV.7) Proof.
This result is the application of the small-noise expan-sion of stochastic differential equations , which is anal-ogous to the equation of variations for ordinary differentialequations. The proof is essentially the extension of the knownresult for the vector-valued, autonomous case , to also allowfor nonautonomous and parameter-dependent SDE-s. For de-tails, see Appendix VIII B. Remark . The zeroth-order SDE in ε , Eq. (IV.6), is preciselythe idealized model. Hence, the solution process x t is deter-ministic and could be also written as x t ≡ x ( t ) . Theorem 3.
Let φ tt ( x ) be the normalized fundamental ma-trix solution to (III.8). Then, η t defined as the solution to thelinear SDE (IV.7) , is an Ornstein-Uhlenbeck process that canbe written as η t = (cid:90) tt φ ts (cid:0) x s (cid:1) g (cid:0) x s , s (cid:1) ds + (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s , s (cid:1) dW s . (IV.8) Proof.
This result is well-known for scalar stochastic differen-tial equations. The extension to our multi-dimensional settingis given in Appendix VIII C.Following this result, let N ( t ) = || η t || denote the norm ofthe vector valued stochastic process η t , which measures theleading-order trajectory uncertainty arising from both deter-ministic and stochastic modeling errors. The leading-ordertrajectory uncertainty is then ε N ( t ) . Using formula (IV.8) for η t , we can define the deterministic term ( N d ), the stochasticterm ( N s ) and the mixed term ( N m ) of this leading-order tra-jectory uncertainty. To remain consistent with the notation ofSection III, we have the deterministic term δ ( x , t ) = ε N d ( t ) .The full expression for N ( t ) is: N ( t ) (IV.9) = (cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) g (cid:0) x s , s (cid:1) ds + (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s , s (cid:1) dW s (cid:19) = (cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) g (cid:0) x s , s (cid:1) ds (cid:19) + (cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s , s (cid:1) dW s (cid:19) + (cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s , s (cid:1) dW s (cid:19) (cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) g (cid:0) x s , s (cid:1) ds (cid:19) = N d ( t ) + N s ( t ) + N m ( t ) Formula (IV.9) allows us to formulate a stochastic extensionof Theorem 1, which applies even in the stochastic setting.The quantity to be estimated is now the mean-square of theleading-order trajectory uncertainty.
Theorem 4 (Bound on the mean-squared leading-order tra-jectory uncertainty) . The leading-order trajectory uncertaintycan be bounded in the mean-square sense as ε E (cid:2) N ( t ) (cid:3) ≤ (cid:18) (cid:90) tt (cid:113) Λ ts ( x s ) ds (cid:19) ∆ ∞ ( x , t ) (IV.10) + (cid:90) tt tr (cid:2) C ts (cid:0) x s (cid:1)(cid:3) ds ∆ σ ∞ ( x , t ) , where we have introduced the notation ∆ σ ∞ ( x , t ) = ε max s ∈ [ t , t ] tr (cid:104) σ (cid:0) x s , s (cid:1) T σ (cid:0) x s , s (cid:1)(cid:105) .Proof. The proof consists of a computation of the expectedvalues of N s and N m , since N d is purely deterministic and wasalready computed before. The details of the proof are given inAppendix VIII D.Note that if the model has no stochastic error, i.e., σ ( x , t ) ≡ , Theorem 4 gives N ( t ) = N d and ∆ σ ∞ ( x , t ) ≡ ε E [ N ( t )] = ε N ( t ) = δ ( x , t ) ≤ (cid:82) tt (cid:112) Λ ts ( x s ) ds ∆ ∞ ( x , t ) . This is consistent with the upperbound derived in Section III.Rearranging expression (IV.10), we obtain a quantity, com-puted in terms of the idealized model and the relative strengthof errors (deterministic or stochastic). We refer to this quan-tity as Model Sensitivity (MS) , defined asMS tt ( x ; r ) : = (cid:18) (cid:90) tt (cid:113) Λ ts ( x s ) ds (cid:19) + r (cid:90) tt tr [ C ts (cid:0) x s (cid:1) ] ds , (IV.11)where r : = ∆ σ ∞ ( x , t ) / ∆ ∞ ( x , t ) is the ratio characterizing therelative importance of the stochastic modeling errors. By cal-culating MS tt for several initial conditions in a phase-spaceregion of interest, we can quickly identify locations of highsensitivity to modeling errors. By Theorem 4, these locationsare expected to show higher uncertainty.Moreover, by Theorem 4, the leading order trajectory un-certainty is related to MS, in the mean-square sense, by ε E [ N ( t ) ] ≤ MS tt ( x ; r ) ∆ ∞ ( x , t ) . (IV.12)In other words, MS is the coefficient relating the leading or-der mean-squared trajectory uncertainty to the modeling un-certainty.As in the purely deterministic case, we obtain a theoremthat relates ( MS tt ) to the proper trajectory uncertainty. Theorem 5 (Universal bound on the mean-squared trajectoryuncertainty) . Consider the stochastic dynamical system de-fined over a finite time interval [ t , t ] and on a compact do-main U ⊂ R n , by the SDE (IV.2) . Then, for any δ > thereexists an ε > , such that for ε < ε the following inequalityholds for all t ∈ [ t , t ] and x ∈ U: (cid:113) E ( | x ε t − x ( t ) | ) ≤ ∆ ∞ ( x , t ) (cid:16)(cid:113) MS tt ( x , r ) + δ (cid:17) . (IV.13) Proof.
See Appendix VIII E.By Theorem 5, the bound on the mean-squared leading-order trajectory uncertainty is extended to the actual mean-squared trajectory uncertainty, for small enough ε . Then, theMS can be used to calculate a time-dependent upper bound onthe trajectory uncertainty, which will be true for any perturba-tion of size ∆ ∞ , assuming a ratio of r between stochastic anddeterministic modeling errors.We also note that in practice, the bound (IV.13) tends to besatisfied even without including the size- δ correction (simi-larly to Theorem 1). This means that the much simpler ex-pression of Theorem 4 can be used to assess the mean-squaredtrajectory uncertainty. In the next section, we demonstrate thisfact on a few examples. V. COMPUTATION OF TRAJECTORY UNCERTAINTYESTIMATES
We start by an explicit calculation of MS for linear systems.Within this class of systems, we can find examples provingthe optimality of our estimates. Consider the constant co-efficient linear stochastic differential equation, driven by an n − dimensional Wiener-process W t , d x t = Ax t dt + ε b dt + εσ d W t , x , b ∈ R n , A , σ ∈ R n × n . (V.1)Here, b is a (constant) deterministic perturbation vector, σ is the covariance matrix of the noise and ε ≥ ∆ ∞ = ε | b | ∆ σ ∞ = ε || σ || F . (V.2)The equation of variations of system (V.1) is simply ˙ φ tt = A φ tt , which gives φ tt = e A ( t − t ) for the flow-map gradient.Then, by formula (IV.11), MS isMS tt = (cid:32) (cid:90) tt (cid:114) Λ (cid:104)(cid:0) e A ( t − s ) (cid:1) T e A ( t − s ) (cid:105) ds (cid:33) + || σ || F | b | (cid:90) tt tr (cid:20)(cid:16) e A ( t − s ) (cid:17) T e A ( t − s ) (cid:21) ds . (V.3)From this, we can obtain the bound on the leading-ordertrajectory uncertainty after multiplying by ε | b | .On the other hand, we can calculate the trajectory uncer-tainty directly. The idealized system (with ε =
0) has the gen-eral solution x t = e A t x , while the solution to the perturbedproblem is the stochastic process x t = e A t x + ε (cid:90) tt e A ( t − s ) b ds + ε (cid:90) tt e A ( t − s ) σ d W s . (V.4)The mean-square of the difference between the idealizedmodel solution, and the real solution is E ( | x ε t − x t | )= ε (cid:18) (cid:90) tt e A ( t − s ) b ds (cid:19) + ε E (cid:18) (cid:90) tt e A ( t − s ) σ d W s (cid:19) = ε (cid:18) (cid:90) tt e A ( t − s ) b ds (cid:19) + ε (cid:90) tt || e A ( t − s ) σ || F ds . (V.5)Here, we used that the expected value of the mixed term iszero and the expression for the second integral follows fromItô’s isometry.An immediate consequence of this calculation is the opti-mality of Theorem 4. If system (V.1) is a scalar equation, x t ∈ R , A , σ , b ∈ R , then once we evaluate the integrals, weobtain E ( | x ε t − x t | ) = ε b A (cid:16) e A ( t − t ) − (cid:17) + ε σ A ( e A ( t − t ) − )= ε b MS tt . (V.6)This shows that Theorem 4 is optimal: the bound it providescannot be strengthened for general systems. A. Numerical examples
Example 1.
The Duffing oscillatorTo illustrate our main results, we apply formula (IV.12) totwo models of differring complexity. First, let us consideronce again the damped-driven Duffing oscillator, defined by(III.23). In the presence of a deterministic, time-periodic per-turbation, the trajectory uncertainty was already shown in Fig.1. To assess the sensitivity to general, possibly stochastic per-turbations, we first calculate MS tt and display it on a uniformgrid over the domain U = [ − . , . ] × [ − . , . ] , for twotime intervals of interest, [ , π ] and [ , π ] . This calculationonly requires knowledge of the idealized system and the rela-tive magnitude of modeling errors. For the calculation of theCauchy-Green strain tensor, we use finite differences, over asecondary grid to increase accuracy.MS fields are shown in Fig. 2. We now assume a spe-cific modeling error that contains both a deterministic and astochastic component. The equations then are SDEs, whichread as dx t = ydt , (V.7) dy t = ( x t − x t − δ y t + A cos t ) dt + ε sin ( ω p t ) dt + ε dW t . In this case, both types of errors are assumed to be of norm ε ,that is, ∆ ∞ = ∆ σ ∞ = ε , with ω p = ).As noted earlier for the calculation of MS tt ( x , r ) , as well asfor the leading-order trajectory uncertainty, we did not makeany assumptions on the form of the modeling errors. For thisreason, given one specific instance of modeling errors and twopoints x ( ) and x ( ) , the relationMS tt (cid:16) x ( ) , r (cid:17) < MS tt (cid:16) x ( ) , r (cid:17) (V.8)does not imply that the actual trajectory uncertainty will begreater in x ( ) than in x ( ) . Instead, what we can concludefrom (V.8), is that the dynamics at x ( ) is such, that it canallow higher trajectory uncertainty than at x ( ) . FIG. 2. Model sensitivity (MS) for the Duffing oscillator under bothdeterministic and stochastic modeling errors. The value of the MSis obtained from formula (IV.11), applied to system (III.23) withparameters δ = . , A = .
3. Both the deterministic model errorand the noise is assumed to have amplitude ε , that is ∆ ∞ ( x , t ) = ε , ∆ σ ∞ ( x , t ) = ε , with ε = .
01. In the upper panel, the time interval ofinterest is [ , π ] , while in the lower panel, it is [ , π ] . Light bluedots mark the starting points of the trajectories relevant for Fig. 3. Example 2.
The Charney- DeVore modelNext, we turn to a higher-dimensional model . It isdemonstrated in Ref. 42, that a six dimensional reduced ordermodel for barotropic flow over topography admits multipleequilibria, and can even exhibit tipping transitions betweenthem. Therefore, the Charney-DeVore model is expectedto show highly unstable transient behavior , which results inhigh sensitivity with respect to perturbations. The dynamicalequations are ˙ x = (cid:101) γ x − C ( x − x ∗ ) , ˙ x = − ( α x − β ) x − Cx − δ x x , ˙ x = ( α x − β ) x − γ x − Cx + δ x x , ˙ x = (cid:101) γ x − C ( x − x ∗ ) + λ ( x x − x x ) , ˙ x = − ( α x − β ) x − Cx − δ x x , (V.9)˙ x = ( α x − β ) x − γ x − Cx + δ x x . The coefficients α m , β m , γ m , δ m are defined by α m = √ π m m − b + m − b + m , β m = β b b + m , δ m = √ π b − m + b + m , (cid:101) γ m = γ m m − √ b π , λ = √ π , γ m = γ m m − √ b π ( b + m ) . (V.10)As in Refs. 42 and 43, we set the parameters, tocorrespond to the multistable regime: ( x ∗ , x ∗ , C , β , γ , b ) =( . , − . , . , . , . , . ) .The MS tt field is shown in Fig. 4 along a few slices of phasespace. Similarly to the low-dimensional Duffing oscillator, thephase space of the Charney-DeVore model also exhibits highvariability for MS.Although we cannot show the complete MS field in thishigh-dimensional phase space, this example demonstrateshow our method remains applicable in higher-dimensionalsystems. Even in this lower-dimensional representation, wecan distinguish structures, with particularly high sensitivity toperturbations over the chosen time scale.Next, we fix a modeling error to the equations (V.9) in theform g ( x , t ) = b sin ( k · x ) cos ( ω p t ) , | b | = , σ = √ I . (V.11)This represents a deterministic modeling error that is periodicin both time and space. In addition, each coordinate is per-turbed by an independent Wiener-process.The vector b is of unit length and has components b =( . , . , . , . , . , . ) . The wave vector is k = ( . , . , . , . , . , . ) , ω p =
10. Withthis choice of the parameters, the magnitude of the pertur-bations is once again ∆ ∞ ( x , t ) = ε | b | = ε and ∆ σ ∞ = ε tr σ T σ = ε , r = ε .While the actual mean-squared trajectory uncertainty is over-estimated for the interval [0,15] (by a factor of around 10),the insets of Fig. 5 show that the bound is relatively tight overshorter time intervals. . . . . . . S q u a r e d E rr o r Estimate from Model SensitivityMean-Squared Error 0 5 10 15Time0 . . . . . . S q u a r e d E rr o r . . . . . . S q u a r e d E rr o r FIG. 3. Square of the difference between the idealized model solutions and the perturbed solutions to the Duffing-system (V.7). The greycurves show the error along a few sample paths, the black curve is the mean-squared error computed from 2000 sample paths. The redcurve is the upper bound on the mean-squared leading-order trajectory uncertainty, defined by MS tt ( x , r ) ∆ ∞ ( x , t ) . The modeling errorsare detailed in the text, ∆ ∞ = ε with ε = .
01 and r =
1. The initial conditions are: left panel: x = ( − . , − . ) , middle panel: x = ( − . , − . ) , right panel: x = ( . , − . ) FIG. 4. Model sensitivity (MS), computed for the Charney-DeVore model (V.9). The parameter values used are given in the text, the timeinterval is t = , t = . The ratio of the importance of stochastic and deterministic modeling errors was set to r =
1. The figures showdifferent two-dimensional slices of the six-dimensional phase space. Light blue dots mark the starting points of the trajectories relevant forFig. 5.
VI. GEOMETRIC STRUCTURE OF THE MODELSENSITIVITY
Geometric descriptions of uncertainty in dynamical sys-tems involve the finite-time Lyapunov exponent. This quan-tity describes the growth rate of infinitesimal perturbations toinitial conditions.
Ridges of the FTLE field often signal re-pelling material surfaces in the phase space . Under furtherassumptions , one can rigorously conclude the presence of arepelling hyperbolic LCS from an FTLE ridge.Our results show, that the FTLE field in itself is not suffi-cient to characterize sensitivity to modeling errors in dynam-ical systems. By Theorem 4, one needs to integrate the timedependent FTLE field over the time interval of interest, to ob-tain MS. Regardless, the two fields are clearly related.For purely deterministic perturbations, let us denote themaximal eigenvalue of the Cauchy-Green strain tensor, com- puted over [ s , t ] at the point x by Λ s , t ( x , t ) = Λ ts (cid:0) F st ( x ) (cid:1) , (VI.1)where F st ( x ) is the flow-map of the idealized model (II.2).Specifically, with this notation, we can write the FTLE asFTLE tt ( x ) = log (cid:112) Λ t , t ( x , t ) t − t . (VI.2)We will use quantity Λ s , t to connect features of the MSfield to those of the FTLE field. Such a connection is alreadysuggested by Fig. 6, which compares the FTLE field of theCharney-DeVore model to the field12 ( t − t ) log MS tt ( x , r ) (VI.3)0 S q u a r e d E rr o r ε = 0 . × ε = 0 . . . . . . S q u a r e d E rr o r ε = 0 . ε = 0 . . . . . . S q u a r e d E rr o r × − ε = 0 . × − ε = 0 . FIG. 5. Square of the difference between the idealized and thereal solutions to the Charney-DeVore model. Grey lines showthe squared error along sample paths, the black curve is the av-erage computed from 1000 sample paths. The red curve is theupper bound on the mean-squared leading-order trajectory uncer-tainty. The time interval [ , ] is magnified in the insets. The tra-jectories start from the point x ( ) = ( , − . , , − . , , ) [ x ( ) = ( . , , , − . , , ) ] in the left [right] column. Themagnitude of the perturbations is ε , which is indicated above thepanels. r = of the same model. That is, we display MS on a similar scaleas the FTLE for a better comparison. This scale will be justi-fied later.The figure shows that the main features of the FTLE fieldare also found in the MS field, if they are compared over thesame time interval. Specifically, the main organizers of thedynamics, the FTLE-ridges, tend to persist in the MS field. Acloser look reveals, however, that this is not always the case.For example, in the region around ( x = . , x = − . ) , finerridges of the FTLE field disappear in the MS field.To analyze this phenomenon, we adopt the following defi-nition of a ridge from Ref. 45. Definition 1.
Let f : R n → R be a smooth function and M ⊂ R n a compact, codimension-one manifold with a bound-ary ∂ M . The manifold M is a ridge of the scalar field f if both M and ∂ M are normally attracting invariant manifolds for thegradient system ˙ x = ∇ f ( x ) . (VI.4)The term normally attracting invariant manifold refers toan invariant manifold for which contraction along the mani-fold is dominated by contraction normal to it. This allows theuse of results that guarantee the persistence of ridges undersmall perturbations to the scalar field f .To find a condition relating the MS- and FTLE-ridges, weassume that the deterministic modeling error in (IV.11) is theonly contributor to MS, that is, MS tt ( x , r ) with r =
0. In thiscase, we can writeMS tt ( x , ) = (cid:18) (cid:90) tt (cid:113) Λ s , t ( x , t ) ds (cid:19) = Λ t , t ( x , t ) (cid:32) (cid:90) tt (cid:115) Λ s , t ( x , t ) Λ t , t ( x , t ) ds (cid:33) . (VI.5)Taking the logarithm and using expression (III.18) for theFTLE field, we obtainlog MS tt ( x , ) ( t − t ) = FTLE tt ( x ) + t − t log (cid:90) tt (cid:115) Λ s , t ( x , t ) Λ t , t ( x , t ) ds . (VI.6)Equation (VI.6) shows that we are able to write the appro-priately scaled MS field as a perturbation of the FTLE field.The difference between the MS and FTLE fields is shown inFig. 7, which suggests the values of the two fields differ sub-stantially, even at the location of the persisting ridge. We con-clude that in general, ridges of the MS field are different fromthose of the FTLE field. As seen from (VI.6), the MS fieldmust be treated as a finite-size-perturbation to the FTLE field:the persistence results of Ref. 45 do not apply.The general results on persistence of normally hyperbolicinvariant manifolds state that for a ridge of a scalar field f ( x ) to persist in the field f ( x ) , the appropriate gradientvector fields in (VI.4) must be C - θ close, for θ small enough.In our setting, this translates into a condition on the gradientof the difference field, defined by (VI.6). This indicates thatridges of the FTLE field, along whichsup x ∈ V (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∇ t − t log (cid:90) tt (cid:115) Λ s , t ( x , t ) Λ t , t ( x , t ) ds (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ θ , sup x ∈ V (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ∇ t − t log (cid:90) tt (cid:115) Λ s , t ( x , t ) Λ t , t ( x , t ) ds (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ θ , (VI.7)holds for θ sufficiently small, are expected to be close toridges of the scalar field logMS tt ( x , ) ( t − t ) . VII. CONCLUSIONS
We have investigated the effect of modeling uncertaintieson trajectories of a dynamical system. Under general smooth-ness assumptions, in a deterministic setting, we derived a1
FIG. 6. Comparison between FTLE and MS fields for the Charney-DeVore model. The scalar fields are shown for the time interval [ , ] ,over the x − x plane. In the left panel, the FTLE field, while in the right panel, the field ( t − t ) logMS tt ( x , ) is plotted.FIG. 7. Difference between the MS and the FTLE fields in theCharney-DeVore model. The quantity logMS tt ( x , ) ( t − t ) − FTLE tt ( x ) isshown over the x − x plane, for the time interval [ , ] . bound on the leading order trajectory uncertainty which canbe computed using the idealized model dynamics and assum-ing a bound on the magnitude of the modeling error. Ourupper bound depends on the eigenvalues of an appropriateCauchy-Green strain tensor of the idealized model, allowingfor a location-specific assessment of trajectory uncertainty inthe phase space.We have also generalized our result to the case of stochas-tic modeling errors. In that setting, we have introduced the Model Sensitivity (MS) , a coefficient relating the modeling un-certainty to the bound of the mean-squared leading-order tra-jectory uncertainty. This MS is computed solely in terms ofthe idealized, deterministic dynamics as a general functional of the invariants of the Cauchy-Green strain tensor. As a con-sequence, we do not need to assume any specific form for themodeling errors to quickly assess their effects. Contrary toprior, statistical and data based methods, MS quantifies tra-jectory sensitivity based on the dynamical properties of theknown model.We have also shown that our bounds on the leading ordertrajectory uncertainty are optimal. Specifically, for a classof linear systems, we gave an example in which the mean-squared trajectory uncertainty was exactly equal to the prod-uct of the MS and the modeling uncertainty. Therefore, forgeneral systems satisfying our smoothness assumptions, ourbounds cannot be improved.On numerical examples, we showed that the MS field canexhibit complex structure in the phase space, which allows usto distinguish particularly sensitive regions. We have shownthat the MS fields, in general, are similar to the FTLE fields,which are often used to characterize instability in phase space.This is in line with the usual reasoning, that the instabilitieswithin a dynamical system typically grow with the rate of thelargest Lyapunov exponent. Our results make this argumentprecise, by pointing out the exact relationship between MSand FTLE. In particular, we find that not all features of theFTLE field persist in the MS field.The sensitivity analysis developed here becomes more com-putationally intensive for higher-dimensional dynamical sys-tems. Indeed, the calculation of the Cauchy-Green straintensor becomes problematic for even a few hundred dimen-sions. This could be improved by using approximate methodssuch as OTD (Optimally Time Dependent) modes which en-able the calculation of the dominant Cauchy-Green eigenvaluewith much less effort. This approach could give useful resultseven for certain climate models, for which sensitivity analysisis critical.2 ACKNOWLEDGEMENTS
We are grateful to Hessam Babaee for his help with the im-plementation of the Charney-DeVore model.
DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openlyavailable in the repository Ref. 47. It contains the numericalimplementation of the examples discussed in the text.
VIII. APPENDIXA. Proof of Theorem 1
The trajectory uncertainty, based on a power-series expan-sion, is (cid:12)(cid:12) x ( t ) − x ε ( t ) (cid:12)(cid:12) = | εη ( t , t , x ) + O ( ε ) | . (VIII.1)This can be bounded by | εη ( t , t , x ) + O ( ε ) | ≤ ε | η ( t , t , x ) | + | O ( ε ) | , (VIII.2)where the remainder O ( ε ) term can be bounded by M ε ,when ε ≤ ¯ ε for some ¯ ε >
0. This bound depends on x , butassuming a compact domain U within R n , we can choose theconstants M > ε such that this bound is satisfied for all x ∈ U .Next, we simply substitute the bound (III.12) for theleading-order trajectory uncertainty into Eq. (VIII.2) to ob-tain (cid:12)(cid:12) x ( t ) − x ε ( t ) (cid:12)(cid:12) ≤ (cid:18) (cid:90) tt (cid:113) Λ ts ( x ( s )) ds (cid:19) ∆ ∞ ( x , t ) + M ε = (cid:18) (cid:90) tt (cid:113) Λ ts ( x ( s )) ds (cid:19) ∆ ∞ ( x , t ) + ∆ ∞ ( x , t ) M ε ∆ ∞ ( x , t )= (cid:18) (cid:90) tt (cid:113) Λ ts ( x ( s )) ds + M ε ∆ ∞ ( x , t ) (cid:19) ∆ ∞ ( x , t ) (VIII.3)Comparing Eq. (VIII.3) to Eq. (III.22) we obtain that if δ < min t ∈ [ t , t ] x ∈ U M ¯ ε ∆ ∞ ( x , t ) = ¯ δ , (VIII.4)then we can set ε : = max t ∈ [ t , t ] x ∈ U (cid:114) ∆ ∞ ( x , t ) δ M . (VIII.5)Otherwise, if δ ≥ ¯ δ , inequality (VIII.3) is satisfied for ε : = ¯ ε .Hence, for all δ > ε : = min max t ∈ [ t , t ] x ∈ U (cid:114) ∆ ∞ ( x , t ) δ M , ¯ ε (VIII.6)as claimed. B. Proof of Theorem 2
The statement follows from an asymptotic expansions forstochastic differential equations of the form dX ε t = µ ( X ε t ) dt + Σ ε ( X ε t ) dW t . (VIII.7)Assume that the coefficient functions µ and Σ ε have boundedand measurable derivatives up to order m . Then, there exists¯ ε >
0, such that for ε < ¯ ε one can recursively obtain stochasticdifferential equations for the stochastic variables X ( k ) t . For any k < mX ε t = X t + ε X t + ε X t + ... + ε k − X k − t + ε k R k ( t , ε ) , (VIII.8)with the remainder term being bounded in the mean-squaredsense. For a proof of the n dimensional case, see Ref. 39.To generalize this result for time- and ε -dependent drift anddiffusion coefficients, we proceed by introducing an SDE onthe extended phase space R n × [ t , t ] × [ , ¯ ε ] of Eq. (IV.2). Let X ε t ∈ R n × [ t , t ] × [ , ¯ ε ] , X ε t = ( x t , , x t , , ..., x t , n , t , ε ) and µ ( X ε t ) = f ( x t , t ) + ε g ( x t , t , ε ) f ( x t , t ) + ε g ( x t , t , ε ) ... f ( x t , t ) n + ε g ( x t , t , ε ) n , Σ ε ( X ε t ) = εσ ( x t , t ) · · · εσ n ( x t , t ) εσ n ( x t , t ) · · · εσ nn ( x t , t ) · · · · · · . (VIII.9)The resulting SDE has the desired form of Eq. (VIII.7) andthe coefficients retain the analytic properties of the functions f and σ , i. e. they remain measurable and have boundedderivatives. This means we can apply the result of Ref. 39. toobtain the following first-order expansion X ε t = X t + ε X t + ε R ( t , ε ) (VIII.10)for the solutions. The coefficients in the expansion are gov-erned by the following set of linear SDEs: dX t = µ ( X t ) dt , X t = t =( x , t , ) , dX t = ∇ µ ( X t ) X t dt + ∂ Σ ε ∂ ε (cid:12)(cid:12)(cid:12)(cid:12) X t dW t , X t = t =( , t , ) . Setting X t = ( x t , t , ) , X t = ( η t , t , ) and keeping only thethe first n entries of the vectors yields the following expansion3for the nonautonomous system (IV.2): dx t = f (cid:0) x t , (cid:1) dt , x t = t = x , d η t = ∇ f (cid:0) x t , t (cid:1) η t dt + g (cid:0) x t , t ; 0 (cid:1) dt + σ (cid:0) x t , t (cid:1) dW t , η t = t = , (VIII.11)as claimed. C. Proof of Theorem 3
We seek a solution of the inhomogeneous, linear SDE(IV.7) using the method of ’variation of constants’ on the so-lution of the homogeneous equation. Let the solution of thecorresponding homogeneous equation be x H ( t , x ) = ϕ ( t ) x . (VIII.12)Here, ϕ ( t ) is the (linear) flow map of Eq. (III.8), map-ping initial conditions at time t to their position at time t .By the method of variation of constants, let η t be of the form η t = ϕ ( t ) x t for some random variable x t . We now computethe differential d η t , keeping in mind that η t is a vector-valuedstochastic process, requiring the use of Itô’s formula. How-ever, since η t = η ( t , x ) is only linear in the x variable, wesimply have d η t = ˙ ϕ ( t ) x t dt + ϕ ( t ) dx t . (VIII.13)Substituting Eq. (IV.7) into Eq. (VIII.13) and noting that ϕ isthe fundamental matrix-solution to the equation of variations(III.8) yields d η t = ˙ ϕ ( t ) x t dt + ϕ ( t ) dx t = ∇ f (cid:0) x t , t (cid:1) x t dt + ϕ ( t ) dx t , ϕ ( t ) dx t = g (cid:0) x t , t ; 0 (cid:1) dt + σ (cid:0) x t , t (cid:1) dW t , (VIII.14) dx t = ϕ ( t ) − (cid:2) g (cid:0) x t , t ; 0 (cid:1) dt + σ (cid:0) x t , t (cid:1) dW t (cid:3) . The last expression in Eq. (VIII.14) is an Itô-integral, whichcan be evaluated as x t = (cid:90) tt ϕ ( s ) − g (cid:0) x s , s ; 0 (cid:1) ds + (cid:90) tt ϕ ( s ) − σ (cid:0) x s , s (cid:1) dW s . (VIII.15)Using the form of η t and observing that φ ts (cid:0) x s (cid:1) = ϕ ( t ) ϕ ( s ) − is the normalized fundamental matrix solution to Eq. (III.8),we obtain η t = (cid:90) tt ϕ ( t ) ϕ ( s ) − g (cid:0) x s , s ; 0 (cid:1) ds + (cid:90) tt ϕ ( t ) ϕ ( s ) − σ (cid:0) x s , s (cid:1) dW s (VIII.16) = (cid:90) tt φ ts (cid:0) x s (cid:1) g (cid:0) x s , s ; 0 (cid:1) ds + (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s , s (cid:1) dW s , which proves the statement of Eq. (IV.8). D. Proof of Theorem 4
First, we compute E ( N ) . By the properties of the Itô-integral, the expected value of the mixed term in Eq. (IV.9)is 0, and hence E ( N )= E ( N d ) + E ( N s ) + E ( N m )= N d + E ( N s )+ (cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) g (cid:0) x s , s ; 0 (cid:1) ds (cid:19) E (cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s s (cid:1) dW s (cid:19) = N d + E ( N s ) . (VIII.17)For the stochastic part of the mean-squared leading-ordertrajectory uncertainty, we utilize Itô’s isometry component-wise to obtain E ( N s )= E (cid:34)(cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s , s (cid:1) dW s (cid:19) (cid:35) = E (cid:20)(cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s , s (cid:1) dW s (cid:19) (cid:18) (cid:90) tt φ ts (cid:0) x s (cid:1) σ (cid:0) x s , s (cid:1) dW s (cid:19)(cid:21) = E (cid:34) ∑ i , j , k , l , m (cid:18) (cid:90) tt ( φ ts ) i j σ jk ( dW s ) k (cid:19) (cid:18) (cid:90) tt ( φ ts ) il σ lm ( dW s ) m (cid:19)(cid:35) = ∑ i , j , k , l , m E (cid:20)(cid:18) (cid:90) tt ( φ ts ) i j σ jk ( φ ts ) il σ lm [( dW s ) k , ( dW s ) m ] (cid:19)(cid:21) (VIII.18)The notation [( dW s ) k , ( dW s ) m ] refers to the quadraticcovariation of the processes ( dW s ) k and ( dW s ) m . Sincethe components of the n -dimensional Wiener-process are as-sumed to be independent, we have (by Itô’s isometry), E (cid:18) (cid:90) tt ( φ ts ) i j σ jk ( φ ts ) il σ lm [( dW s ) k , ( dW s ) m ] (cid:19) = E (cid:18) (cid:90) tt ( φ ts ) i j σ jk ( φ ts ) il σ lm δ km ds (cid:19) , (VIII.19)where δ km is the Kronecker-delta. Denoting the Frobenius-norm by || · || F : R n × n → R + , we have || A || F = ∑ i , j | A i j | = tr ( A T A ) . Therefore, E ( N s ) = (cid:90) tt || φ ts σ || F ds ε E (cid:2) N s ( t ) (cid:3) ≤ (cid:90) tt ε || φ ts || F ds max s ∈ [ t , t ] (cid:13)(cid:13) σ (cid:0) x s , s (cid:1)(cid:13)(cid:13) F = (cid:90) tt tr (cid:2) C ts (cid:0) x s (cid:1)(cid:3) ds ∆ σ ∞ ( x , t ) . (VIII.20)In Section III. A, we also concluded in Eq. (III.12)that δ ( x , t ) = ε N d ( t ) ≤ (cid:82) tt (cid:112) Λ ts ( x s ) ds ∆ ∞ ( x , t ) . Substituting4Eqs. (III.12) and (VIII.20) into Eq. (VIII.17) implies ε E (cid:2) N ( t ) (cid:3) ≤ (cid:18) (cid:90) tt (cid:113) Λ ts ( x s ) ds (cid:19) ∆ ∞ ( x , t )+ (cid:90) tt tr (cid:2) C ts (cid:0) x s (cid:1)(cid:3) ds ∆ σ ∞ ( x , t ) , (VIII.21)as claimed. E. Proof of Theorem 5
Using the small-noise expansion (IV.4) for the mean-squared trajectory uncertainty, for ε < ¯ ε , we obtain E (cid:0) | x ε t − x ( t ) | (cid:1) = E (cid:0) | εη t + ε R ( t , ε ) | (cid:1) . (VIII.22)Using the Minkowski-inequality, we also find that (cid:113) E ( | x ε t − x ( t ) | ) ≤ (cid:113) ε E ( | η t | ) + (cid:113) ε E ( | R ( t , ε ) | ) . (VIII.23)Since the second order remainder term in Eq. (IV.4) isbounded in the mean-squared sense, we have, for some K < ∞ , sup t ∈ [ t , t ] E ( | R ( t , ε ) | ) ≤ K . (VIII.24)To bound E ( | η t | ) , we use Theorem 4 in the form of Eq. IV.11to obtain ε E ( | η t | ) ≤ MS tt ( x , r ) ∆ ∞ ( x , t ) . (VIII.25)Substituting bounds (VIII.24) and (VIII.25) into the origi-nal expression (VIII.22), we have (cid:113) E ( | x ε t − x ( t ) | ) ≤ (cid:113) MS tt ( x , r ) ∆ ∞ ( x , t ) + (cid:113) K ε (VIII.26)Since x is taken from a compact domain U ⊂ R n , we canchoose the constant K to be independent of x . After rear-ranging the terms, we obtain (cid:113) E ( | x ε t − x ( t ) | ) ≤ (cid:113) MS tt ( x , r ) ∆ ∞ ( x , t ) + ∆ ∞ ( x , t ) K ε ∆ ∞ ( x , t )= (cid:18)(cid:113) MS tt ( x , r ) + K ε ∆ ∞ ( x , t ) (cid:19) ∆ ∞ ( x , t ) . (VIII.27)Comparing (IV.13) to (VIII.27), we obtain the statement ofTheorem 5 after setting ε : = min max t ∈ [ t , t ] x ∈ U (cid:115) ∆ ∞ ( x , t ) δ K , ¯ ε . (VIII.28) R. Kubo, Journal of the Physical Society of Japan , 570 (1957). D. Ruelle, Nonlinearity , 855 (2009). R. V. Abramov and A. J. Majda, Nonlinearity , 2793 (2007). A. J. Majda and X. Wang, Communications in Mathematical Sciences ,145 (2010). V. Lucarini, Journal of Statistical Physics , 1698 (2018). C. E. Leith, Journal of Atmospheric Sciences , 2022 (1975). A. Gritsun and V. Lucarini, Physica D: Nonlinear Phenomena , 62(2017). V. Lembo, V. Lucarini, and F. Ragone, Scientific Reports , 8668 (2020). T. Bódai, V. Lucarini, and F. Lunkeit, Chaos: An Interdisciplinary Journalof Nonlinear Science , 023124 (2020). J. M. Nese, Physica D: Nonlinear Phenomena , 237 (1989). F. Ginelli, H. Chaté, R. Livi, and A. Politi, Journal of Physics A: Mathe-matical and Theoretical , 254005 (2013). E. Kalnay,
Atmospheric modeling, data assimilation and predictability (Cambridge University Press, 2002). E. Ott, B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza,E. Kalnay, D. Patil, and J. A. Yorke, Tellus A: Dynamic Meteorology andOceanography , 415 (2004). C. Grudzien, A. Carrassi, and M. Bocquet, SIAM-ASA Journal on Uncer-tainty Quantification , 1335 (2018). E. Hawkins and R. Sutton, Bulletin of the American Meteorological Society , 1095 (2009). N. Maher, S. Milinski, L. Suarez-Gutierrez, M. Botzet, L. Kornblueh,Y. Takano, J. Kröger, R. Ghosh, C. Hedemann, C. Li, et al. , Journal ofAdvances in Modeling Earth Systems , 2050 (2019). T. Tél, T. Bódai, G. Drótos, T. Haszpra, M. Herein, B. Kaszás, andM. Vincze, Journal of Statistical Physics (2019). M. A. Laughton, Journal of Electronics and Control , 577 (1964). D. G. Cacuci, M. Ionescu-Bujor, and I. M. Navon,
Sensitivity and uncer-tainty analysis, volume II: applications to large-scale systems (CRC press,2005). D. J. Lea, M. R. Allen, and T. W. Haine, Tellus, Series A: Dynamic Mete-orology and Oceanography , 523 (2000). D. J. Lea, T. W. Haine, M. R. Allen, and J. A. Hansen, Quarterly Journalof the Royal Meteorological Society , 2587 (2002). J. Thuburn, Quarterly Journal of the Royal Meteorological Society , 73(2005). Q. Wang, Journal of Computational Physics , 1 (2013). Q. Wang, R. Hu, and P. Blonigan, Journal of Computational Physics ,210 (2014). A. Ni and Q. Wang, Journal of Computational Physics , 56 (2017). D. Lasagna, A. Sharma, and J. Meyers, Journal of Computational Physics , 119 (2019). M. Mathur, G. Haller, T. Peacock, J. E. Ruppert-Felsot, and H. L. Swinney,Physical Review Letters , 144502 (2007). G. Haller, Annual Review of Fluid Mechanics , 137 (2015). I. Arnold, Vladimir,
Ordinary Differential Equations (Springer-Verlag,Berlin, 1992). G. Haller and T. Sapsis, Chaos , 023115 (2011). F. Brauer, Journal of Mathematical Analysis and Applications , 198(1966). U. Kirchgraber, Celestial Mechanics , 351 (1976). J. Guckenheimer and P. Holmes,
Nonlinear oscillations, dynamical systems,and bifurcations of vector fields (Springer-Verlag, New York, 1983). A. Hadjighasem, M. Farazmand, and G. Haller, Nonlinear Dynamics ,689 (2013). B. Øksendal,
Stochastic Differential Equations (Springer-Verlag, NewYork, 2010). Y. N. Blagoveshchenskii, Theory of probability and its applications (1962). M. I. Freidlin and A. D. Wentzell,
Random Perturbations of DynamicalSystems (Springer, 2012). C. W. Gardiner,
Handbook of Stochastic Methods (Springer, 2004). S. Albeverio and B. Smii, Stochastic Processes and their Applications ,1009 (2015). K. Onu, F. Huhn, and G. Haller, Journal of Computational Science , 26(2015). J. G. Charney and J. G. DeVore, Journal of the Atmospheric Sciences , D. T. Crommelin, J. D. Opsteegh, and F. Verhulst, Journal of the Atmo-spheric Sciences , 1406 (2004). H. Babaee, M. Farazmand, G. Haller, and T. P. Sapsis, Chaos (2017). G. Haller, Physica D: Nonlinear Phenomena , 574 (2011). D. Karrasch and G. Haller, Chaos , 043126 (2013). N. Fenichel, Indiana University Mathematics Journal , 193 (1971). https://github.com/balintkaszas/ModelSensitivity . O. Kallenberg,