Stability of Ordinary Differential Equations with Colored Noise Forcing 3 1 0 Timothy Blass1 L.A. Romero 2 2 n January 11, 2013 a J 0 1 Abstract We present a perturbation method for determining the moment sta- ] h bility of linear ordinary differential equations with parametric forcing by p colored noise. In particular, the forcing arises from passing white noise - through an nth order filter. We carry out a perturbation analysis based h on a small parameter ε that gives theamplitude of theforcing. Ourper- t a turbation analysis is based on a ladder operator approach to the vector m Ornstein-Uhlenbeck process. We can carry out our perturbation expan- [ siontoanyorderinε,foralargeclasslinearfilters,andforquitearbitrary linear systems. As an example we apply our results to the stochastically 1 forced Mathieu equation. v 5 Subject Class: Primary: 93E15;Secondary: 60H10, 34D10. 8 0 2 1 Introduction . 1 0 1.1 A Class of Stochastically Forced Linear Equations 3 1 The original goal of this work was to develop a framework for analyzing the : v stability of the stochastically forced Mathieu equation: i X x¨+γx˙ +(ω2+εf(t))x=0, (1) 0 r a where f is a stochastic process, and the stability is determined by the bound- edness of the second moment x2(t) [5, 13]. Here, denotes the sample- hh ii hh·ii average. We wanted to avoid heuristic methods, and consider cases where f(t) is a stochastic process with a realistic power spectral density. In particular, 1DepartmentofMathematicalSciences,CarnegieMellonUniversity,Pittsburgh,PA15213, [email protected] 2ComputationalMathematicsandAlgorithmsDepartment,SandiaNationalLaboratories, MS1320,P.O.Box5800,Albuquerque, NM87123-1320, [email protected] 3Sandia National Laboratories is a multi-program laboratory managed and operated by SandiaCorporation,awhollyownedsubsidiaryofLockheedMartinCorporation,fortheU.S. DepartmentofEnergy’sNationalNuclearSecurityAdministrationundercontractDE-AC04- 94AL85000. 1 we do not want to assume that f is white noise. Hence we want to analyze the case where f(t) is colored noise. However, in order to rigorously derive a Fokker-Planck equation for a stochastic differential equation, the governing equation must include only white noise [5]. We can achieve both goals of rigor and realistic power spectral density by letting f be the output of a linear filter that is forced by a vector white noise ξ. That is, s˙ =Hs+ξ(t), (2) f(t)= a,s(t) , (3) h i whereHisann nreal,diagonalizablematrix,whoseeigenvalueshavenegative real parts, a ×Rn, and , is the standard inner product on Cn. We will ∈ h· ·i take deterministic initial condition s(0) = 0. We assume the noise vector ξ is weighted white noise, meaning ξ =0, ξ(t+τ)ξT(t) =Bδ(τ), (4) hh ii hh ii where B is symmetric and positive semi-definite. Thus, when s(t) solves (2) it is a standard vector-valued Ornstein-Uhlenbeck process. We refer to the scalar process, f(t) = a,s(t) , as colored noise or as an nth-order filter provided h i s(t) solves (2). We will make only mild requirements on the matrices H and B, thereby allowing for wide variability in the power spectral density of the resultingprocess a,s(t) . Thus,inallowingforawiderangeofchoicesofH,B, h i anda,ourapproachaccommodatesabroadclassofcolorednoiseforcingterms. In this paper we will be concerned with the more general problem of linear equations that are being parametrically forced by the function f(t) in equation (3). That is equations of the form x˙ =A x+ε a,s A x, (5) 0 1 h i wheresisthesolutiontothestochasticequation(2),x(t) RN forsomeN 1, ∈ ≥ and A , A are N N constant matrices. 0 1 × The purpose of this paper is to present a perturbation method (assuming ε is small) for determining the stability of the solution x(t) of (5), by which we mean the boundedness of the second moments of x(t). However, our method applies to the pthmoment, so we willnot limit ouranalysisto secondmoments only. Van Kampen has presented a heuristic approach to the case of colored noiseforcing,[24]. Thoughderivedbycompletely differentmeans,his resultfor the Mathieu equation (1) is the same as ours when considering only the first- moments, and without damping. He arrives at his result by truncating some seriesatorderε2, andis only expectedto be valid to this order. Ourmethod is rigorous, can be applied to find solutions to any order in ε, and applies to any moment. We discuss this further in 5.2. § We were originally interested in equation (1) as a model for the response of capillary gravity waves to a time-varying gravitational field arising from ran- dom vertical motions of a container with a free surface (as in [19]). Here f(t) representsthe randomfluctuations in acceleration. Since the Fourier transform 2 of an acceleration should vanish at zero, along with its derivative, the power spectral density of a realistic process f should satisfy S(0) = S′(0) = 0. For example, we can construct a two-dimensional filter using the system (2) that has the power spectral density σβ2ω2 S(ω)= , (6) (ω2+µ2)(ω2+µ2) 1 2 by choosing µ 0 σ 0 β H= −β1 µ , B= 0 0 , a= µ , µ1,µ2,σ >0. 2 2 (cid:18) − (cid:19) (cid:18) (cid:19) (cid:18) − (cid:19) The formula for S(ω) in equation (6) follows from Corollary 7 in Appendix C. The stochastically forced Mathieu equation has been analyzed before, for instancein[2,3,4,9,16,22]but notforthe case(1), orthe generalsetting (5). In[16,22]theyconsideradditiveforcing,andin[2,3,4]theyconsideradifferent typeofparametricforcing. In[9]theyconsideradifferentclassofcolorednoises and study stability by truncating an infinite hierarchy of moment equations. Other studies concern Lyapunov stability, or rely on numerical methods. Our analysis applies to a broad class of equations (5) with a wide variety of forcing terms (2), is semi-analytical (relying only on numerics for the computation of eigenvalues of small matrices), and can be applied to any moment. 1.2 Ladder Operators and the Vector Ornstein-Uhlenbeck Process Our perturbation analysis of the moment stability of equation (5) relies heav- ily on a simple characterization of the eigenvalues and eigenfunctions of the Fokker-Planck equation associated with equation (2). In particular, in 2 and § 3 we characterize the spectrum using ladder operators by generalizing Dirac’s § creationandannihilationoperatorapproachtothequantumharmonicoscillator [11]. Anunderstandingofthespectrumandeigenfunctionsintermsofitsladder operators is crucial to developing the perturbation theory in 4. Though other § authorshaveusedladderoperatorsforOrnstein-Uhlenbeckprocesses,theyhave only considered the scalar case n=1 [20, 21, 23, 25]. We believe the extension to the vector case is by no means trivial, and is interesting in its own right. The probabilitydensity functionP(s ,...,s ,t)associatedwiththe process 1 n s(t) defined by equations (2) and (4) satisfies the Fokker-Planck equation 1 ∂ P = P, with P = div(B P) div(HsP). (7) t D D 2 ∇ − is called the Fokker-Planck operator associated to (2). See [12] for a deriva- D tion of this equation. We note that the Fokker-Planckequation (7) is the same in both the Itoˆ and Stratonovich interpretations because the matrix B is inde- pendent of s (see [12]). The operator will play a crucial role in our stability D analysis. 3 In 2we beginby analyzingthe operator in terms ofits associatedladder § D operators. That is, operators satisfying the commutator equation L [ , ]=µ (8) D L L As in Dirac’s theory of the harmonic oscillator, the significance of the ladder operatorsstems from the fact that if φ is aneigenfunction of with eigenvalue D λ, then the function φ will either vanish, or be an eigenfunction of with L D eigenvalue λ+µ. In 2weshowthatwecanconstructtheladderoperatorsbysolvingamatrix § eigenvalue problem Ty =µy, T=DA, (9) where A is an antisymmetric matrix and D is a symmetric matrix, expressed in terms of H, B. We show there are n raising operators, , which generate k L new eigenfunctions of with an increase in the real part of the eigenvalue, D and n lowering operators that correspondingly decrease the real part of −k L the eigenvalue. We also show that can be expressed in terms of its ladder D operators. In particular, n = µ . (10) k −k k D L L k=1 X where µ is the increment of the ladder operator . That is, [ , ] = µ . k k k k k L D L L This representation is useful for determining the spectrum of . D In 3 we characterize the solutions of § φ=χφ (11) D intermsoftheladderoperators, ,andincrements,µ,solving(8). Inparticular, L we show that any eigenvalue χ of can be written as D n χk = kjµj (12) − j=1 X where µ are the increments of the ladder operators with positive real parts, j and the k are non-negative integers. We will see that the increments µ are j j the negative of the eigenvalues of the matrix H defining the filter in equation (3). We alsoshow that anyeigenfunction of canbe obtained by applying the D ladderoperatortotheeigenfunctionΦ (s)associatedwiththe eigenvalueχ=0 0 of , which is the eigenvalue with the largest real part. D Theresultssummarizedinthelastparagraphrelyonthefactthatrealparts of the eigenvalues of are bounded above (see Lemma 6) , which is proved in D [17, 18, 22], but we give a different and simple proof of this in Appendix B. Here, the domain of is the set of functions that have bounded moments of D anyorder. The spectrumandeigenfunctionsof havebeenstudiedbefore(see D [17, 18, 22]) but not in the context of ladder operators. 4 1.3 Perturbation Expansion for Moment Stability Analy- sis In 4 we use the classical perturbation theory of eigenvalues to carry out an § analysis of the stability of equation (5). Our analysis begins by considering the ODEs for s(t) and x(t) together as a single ODE system. The probability densityfunctionP(s ,...,s ,x ,...,x ,t)forthecombinedsystem(2)and(5) 1 n 1 N solves the Fokker-Planckequation 1 ∂tP = divs(B sP) divs(HsP) divx((A0+ε a,s A1)xP) (13) 2 ∇ − − h i The notationdivs and s refer to divergence and gradientwith respect to only ∇ the sj variables, and similarly divx is divergence in x variables. Equation (13) isthe sameinboththeItoˆandStratonovichinterpretationsbecausethe matrix B is independent of s and x (see [12]). Wecanderiveanequationforthepthmarginalmomentsbymultiplying(13) by monomials xα and integrating with respect to dx, where α is a multi-index of order p. The result is an equation for m(s,t), a vector of the pth marginal moments, which is of the form ∂ m= m+Γ m+ε a,s Γ m. (14) t 0 1 D h i Note that is a differential operator in the s variables only, D 1 ϕ= divs(B sϕ) divs(Hsϕ). (15) D 2 ∇ − Inequation(14)eachcomponentofm(s,t)isoftheform xαP(x,s,t)dxfor RN some multi-index α with α = p. m indicates applied to each component | | D D R ofm. For much of our analysiswe can assume that the matrices Γ and Γ are 0 1 given to us, but we illustrate how to obtain these matrices from the matrices A and A for the particular case of the Mathieu equation in 5. The matrices 0 1 § Γ , Γ in (14) are constant and depend on which moments one is considering 0 1 N +p 1 (see example in equation (53)). There are J = − distinct pth p (cid:18) (cid:19) order monomials in N variables, therefore Γ and Γ are J J matrices. 0 1 × As in a standard stability analysis, in order to determine the stability of (14), we look for solutions of the form m(s,t) = eλtm(s). Our equation for m(s) becomes λm= m+Γ0me+ε a,s Γ1m. (16) D h i That is, the equation for the pth marginal moments of x(t) can be written as an eigenvalue problem, and stability is decided by the sign of the real part of the largest eigenvalue. We doa perturbationanalysisassumingthatthe magnitude ε ofthe forcing is small. Our analysis relies on the fact that that when ε = 0 the eigenfunc- tions of equation(16) are the direct product of the eigenfunctions of and the D eigenvectors of the matrix Γ . 0 5 A key observation (see Lemma 9) for the perturbation analysis is that for any vector a we can determine constants α and β such that a,s can be k k h i written as n a,s = (α +β ), (17) k k k −k h i L L k=1 X where are the ladder operators satisfying (8). The proof of Lemma 9 is ±k L given in Appendix D. In 4 we show that when ε = 0, the eigenvalue of equation (16) with the § largest real part is the same as the largest eigenvalue of the matrix Γ . If λ is 0 0 this unperturbed eigenvalue, then λ(ε)=λ +λ ε2+... with 0 2 J λ = ψ ,Γ φ ψ ,Γ φ G(ν ν ) (18) 2 h 1 1 jih j 1 1i 1− j j=1 X whereν =λ ,ν aretheeigenvaluesofΓ ,andφ ,ψ aretheeigenvectorsand 1 0 j 0 j j normalized adjoint eigenvectors of Γ . Equation (18) uses the extended power 0 spectral density G(z), which is defined for a general stationary random process in (79), and is given explicitly in (81) for the filter a,s(t) . h i The form of λ in (18) is derived for forcing terms that have the form (3), 2 however, the fact that this is simply a weighted sum of values of G, whose coefficients depend only onΓ ,Γ (which do notdepend onthe filter), suggests 0 1 such a formula could hold for any process with a well-defined extended power spectral density. We have carried the perturbation analysis to higher orders, butthe higher-ordercoefficientsdo notappearto havesuchasimple formas in equation (18). The methodin 4involvesconstructingmatricesΓ andΓ , which, asmen- 0 1 § tioned earlier, depend on A ,A , and the representation of the pth marginal 0 1 moments as a vector. We use the stochastic Mathieu equation as a specific ex- amplein 5. In 5.3wediscussanumericalmethodfordeterminingthestability § § of (5) without assuming that ε is small. We compare these numerical results to our perturbation results up to both secondandfourth order for the Mathieu equation, and show they are in excellent agreement. In 5.4 we give a second § representation(whosederivationisgivenin[8])forλ thatdoesnotinvolvethe 2 matrices Γ and Γ , but deals directly with the matrices A and A . We have 0 1 0 1 found that this representation simplifies numerical computations. 2 Existence and Properties of Ladder Operators Inthis sectionwe define the notionofa ladderoperator,showhow to construct these operators, and prove some basic lemmas about them. Lemma 1 shows howladderoperatorscanbeusedtogenerateneweigenfunctionsthathavetheir eigenvalue changed by the increment of the ladder operator. Lemma 2 shows howtofindtheladderoperators andtheirincrementsµ bysolvingamatrix k k L eigenvalueproblem. Lemma3showsthattheincrementsoftheladderoperators 6 are zero and µ ,where µ are the eigenvalues of the matrix H defining our k k ± − filter(seeequation(2)). Lemma4givesthecommutatorrelationsfortheladder operators, and Lemma 5 shows that the operator can be expressed as a D weighted sum of , where are the ladder operators. Throughout this −k k k L L L section (and the rest of the paper) the operator is that defined in (15). D The basic lemmas in this section are crucial to the rest of this paper, and hence we have tried to write this section so that the lemmas stand out clearly. Though the lemmas are all easily stated, the proofs of some of the lemmas are quite technical, especially when attention is given to ensuring that they apply for complex eigenvalues of the matrix T. For this reason, we have relegated many of the proofs to Appendix A. Before discussing ladder operators it should be noted that we define the domain of as the set of functions that have bounded moments of any order. D Thus, our definition of the domain of differs from that given in [17]. In D that paper they defined the domain based on the exponential decay of the eigenfunction Φ that we define in equation (32) and discuss in 3. The two 0 § definitions of the domain give the identical eigenfunctions, but we believe ours is morenaturalsince itdoes notrequireknowingthe solutionaheadoftime. In [17] they discuss a continuous spectrum that arises if the domain is defined so that the eigenfunctions of are only required to be square integrable (or some D similarly less restrictive condition). An examination of these eigenfunctions shows that they have a power law decay as s goes to infinity, and hence do not have moments of all orders. Hence our definition of the domain also excludes this continuous spectrum. We now give a definition of a ladder operator of . D Definition 1. An operator is a ladder operator for with increment µ if [ , ] = µ for some µ CL, where [ , ] denotes theDcommutator [ , ] = D L L ∈ · · D L . DL−LD The following lemma shows that ladder operators can be used to generate new eigenfunctions from ones that we already know. Lemma 1. Suppose is a ladder operator such that [ , ] = µ . Let φ be L D L L an eigenfunction of with eigenvalue χ. Then either φ = 0, or φ is an D L L eigenfunction of with eigenvalue χ+µ. D Proof. Wehave φ φ =µ φ. Sinceφisaneigenfunctionof ,thisgives DL −LD L D us φ=(χ+µ) φ. DL L We definedthe domainof to bethe setoffunctions thathavemomentsof D allorders. Itshouldbe notedthatLemma1wouldnotapply ifthe domainhad been(forexample)thesetofallsquareintegrablefunctions. Inthatcaseathird possibility wouldexist. It couldbe that the function φ is squareintegrable,but the function φ is not. Thus φ would not generate a new eigenfunction. L L We will show that has 2n+1 ladder operators , k = 0,...,n. We ±k D L begin by decomposing into simple differential and multiplicative operators. D 7 Definition 2. We define the operators L ,j =1,...,2n+1 as follows. j L φ=∂ φ for j =1,...,n j sj L φ=s φ for j =1,...,n (19) j+n j L φ= φ for j =2n+1 2n+1 I Here is the identity operator. Note that [L ,L ] = 0 unless j k = n, j k I | − | and [L ,L ]= . In particular, we have j j+n I [L ,L ]=δ j,k =1,...,n (20) j k+n jk I We note that can be expressed in the operators L as j D 2n+1 1 = d L L . (21) jk j k D 2 k,j=1 X WeletDdenotethesymmetricmatrixwithcomponentsd in(21). Thechoice jk ofd in(21)isnotunique, butwemakeanexplicitchoicethatmakesthis ma- jk trix symmetric. If we let A denote the antisymmetric matrix with components a given by jk [L ,L ]=a , (22) j k jk I then we have explicit expressions for A and D 0 I 0 B H 0 n n − . . A= In 0n .. , D= HT 0n .. . (23) − − 0 0 0 ... tr(H)  ···   −      For details regarding the construction of D see Lemma 14 in Appendix A. Just as has a representationinterms ofthe operatorsL , its ladder oper- j D ators will also be expressed in terms of the L . Consider an operator j 2n+1 = y L . (24) j j L j=1 X We write y for the vector of coefficients of . From the representations (21) L and(24), we see thatthe commutator [ , ] involvessums ofterms of the form D L L L L L L L ,whichdonotatfirstsightappeartobelinearintheoperators i j k k i j − L , m = 1,...,2n+1. However, by twice applying the commutator relations m inequation(20)wecanshowthatL L L L L L isin factasumofthe L . i j k k i j m − Determining the coefficient vector y and increment µ thus becomes a matrix eigenvalue problem. The details of how we arrive at this form are given in Appendix A. Here we will merely state the result of these manipulations. Lemma2. Ifyisthevectorofcoefficientsfor , asdefinedasinequation(24), L then the equation [ , ] = µ can be written as a matrix eigenvalue problem D L L Ty=µy, where T=DA, and D and A are defined in (23). 8 We makethe assumptionthatthe eigenvaluesofHhavenegativerealparts, and the eigenvectors form a complete set. For simplicity of the arguments, we will also assume that the eigenvalues of H aresimple. By explicitly writing out the eigenvalue problem Ty =µy we can determine the eigenvalues µ in terms k of the eigenvalues of the matrix H. We will give the details of the proof in Appendix A. Lemma 3. The eigenvalues of T=DA are 0, µ , k=1,...,n, where µ k k { ± } − are the eigenvalues of the matrix H. Note that is the identity operator with increment 0. Thus, our analysis 0 L only involves the 2n ladder operators for k =1,...,n. ±k L In doing the perturbation expansion it will be necessary to have the com- mutator relations of the operators . Finding the commutator relations for i L [ , ] can be turned into a linear algebra problem involving the eigenvectors j k L L of the matrix T. In particular, using equations (24) and (22) we get [ , ]= yTAy . (25) Lj Lk j k I From equation (23) it is easily seen tha(cid:0)t (cid:1) AD= (DA)T = TT (26) − − If DAy = µy, then multiplying both sides of this by A and using equation (26) we see that Ay is an eigenvector of TT with eigenvalue µ. With this in − mind, the right hand side of equation (25) can be written as the inner product between the vector y and the adjoint eigenvector of T associated with µ . j k − Using the fact that the eigenvectors and adjoint eigenvectors of a matrix form a bi-orthogonalset, we can arrive at a simple expression for the commutators. When dealing with complex quantities, the notation in this argument gets to be a bit tedious, and we will leave the details to Appendix A. The final commutator result is given by the following lemma. Lemma 4. For j,k 1 we have [ , ]=0 and [ , ]=δ . j k −j k jk ≥ L L L L I In Dirac’s theory of the harmonic oscillator,he shows that the Hamiltonian operator can be written as the product of the raising and lowering operators. We now generalize this result to the vector case. In this case the operator D can be written as a weighted sum of the products of the raising and lowering operators. The next lemma shows that the weights are in fact the eigenvalues µ of the matrix T. We leavethe proofof this lemma to Appendix A, but note k that its proof is probably the most subtle one in this paper. Lemma 5. The differential operator = 2n+1 1d L L can be written as D i,j=1 2 i,j i j n P = µ . (27) k −k k D L L k=1 X An important feature of the decomposition (27) is that only terms of the form , k >0, appear (there are no terms of the form for k >0). −k k k −k L L L L 9 3 Eigenvalues and Eigenfunctions of D In this section we will use the ladder operator formalism to completely charac- terize the eigenvalues and eigenfunctions of the operator . We note that the D spectrum of has already been studied and characterized[17, 18, 22], but not D in terms of ladder operators. We include another proofof those results because the characterization in terms of ladder operators is used in the perturbation analysis in 4. § As with Dirac’s theory of the quantum harmonic oscillator, the analysis of the spectrum usingladder operatorsrequiresthatthe realpartofthe spectrum is bounded above. We will now state this as a lemma, but leave the proof to Appendix B. Lemma 6. The real part of spectrum of the operator , as defined in (15), is D bounded above. The following theorem will allow us to characterize the eigenfunction asso- ciated with the eigenvalue with the largest real part. Theorem 1. LetΦ(s)beaneigenfunction of (as inequation (15))associated D with the eigenvalue having the largest real part. We must have Φ = 0 for k L k =1,...,n. Proof. Suppose Φ(s) is an eigenfunction of with eigenvalue χ. If Ψ= Φ= k D L 6 0, then Ψ(s) will be an eigenfunction of with eigenvalue χ+µ . This will k D giveusaneigenvaluewithalargerrealpartthanχ. Henceifχistheeigenvalue with the largest real part, then Φ=0 for all k. k L Remark 1. The system Φ = 0 for each k = 1,...,n is an over-determined k L system of first order differential equations. The fact that a solution exists is non-trivial. However, the fact that [ , ] = 0 implies that the Frobenius k j L L Theorem applies (see [1, 6]), which guarantees the system is solvable. Remark 2. As in the comment following Lemma 1, we should note that the domain of is defined as the set of functions that have moments of all orders. D If the domain of were defined using the less stringent requirement that the D eigenfunctions were square integrable, it would not be necessary that Φ =0 k L for all k. This is because in this case Φ does not have to generate a new k L eigenfunction. Itcouldinsteadproduceafunctionthatisnotsquareintegrable. By Theorem 1, the “top” eigenfunction Φ (s) (i.e. the eigenfunction associ- 0 ated to the largest eigenvalue of ) must satisfy the equations Φ =0. If yk k 0 D L is the eigenvector of T associated with the eigenvalue µ , and if µ = 0, then k k 6 the last component of yk vanishes (see the proof of Lemma 3 in Appendix A). That is, we can write p k yk = q . (28) k   0   10