Linear stochastic differential algebraic equations with constant coefficients

LINEAR STOCHASTIC DIFFERENTIAL-ALGEBRAIC EQUATIONS WITH CONSTANT COEFFICIENTS AURELI ALABERT∗AND MARCO FERRANTE† Abstract. We consider linear stochastic differential-algebraic equations with constant coefficients and additive white noise. Due to the nature of this class of equations, the solution must be defined as a generalised process (in the sense of Dawson and Fernique). We provide sufficient conditions for the law of the variables of the solution process to be absolutely continuous with respect to Lebesgue measure. 1. Introduction A Differential-Algebraic Equation (DAE) is, essentially, an Ordinary Differential Equation (ODE) F(x, ˙x) = 0 that cannot be solved for the derivative ˙x. The name comes from the fact that in some cases they can be reduced to a two-part system: A usual differential system plus a “nondifferential” one (hence “algebraic”, with some abuse of language), that is ½˙x=f(x, y) 0 = g(y, z) for some partitioning of xinto variables x, y, z. In general, however, such a splitting need not exist. In comparison with ODE’s, these equations present at least two major difficulties: the first lies in the fact that it is not possible to establish general existence and uniqueness results, due to their more complicate structure; the second one is that DAE’s do not regularise the input (quite the contrary), since solving them typically involves differentiation in place of integration. At the same time, DAE’s are very important objects, arising in many application fields; among them we mention the simulation of electrical circuits, the modelling of multibody mechanisms, the approximation of singular perturbation problems arising e.g. in fluid dynamics, the discretisation of partial differential equations, the analysis of chemical processes, 2000 Mathematics Subject Classification. 60H10, 34A09. ∗Supported by grants 2001SGR99-00174 of CIRIT and BFM2003-0261 of MCYT. †Partially supported by a grant of the CRM, Bellaterra, Spain and a grant of the GNAMPA, Italy. 1 2 AURELI ALABERT AND MARCO FERRANTE and the problem of protein folding. We refer to Rabier and Rheinboldt [9] for a survey of applications. The class of DAE’s most treated in the literature is, not surprisingly, that of linear equations, which have the form A(t) ˙x(t) + B(t)x(t) = f(t), with x, f :R+→Rnand A, B :R+→Rn×n. When Aand Bare constant matrices the equation is said to have constant coefficients. Recently, there has been some incipient work (Schein and Denk [11] and Winkler [12]) on Stochastic Differential-Algebraic Equations (SDAE). In view to incorporate to the model a random external perturbation, an additional term is attached to the differential-algebraic equation, in the form of an additive noise (white or coloured). The solution will then be a stochastic process instead of a single function. Since the focus in [11] and [12] is on numerical solving and the particular applications, some interesting theoretical questions have been left aside in these papers. Our long-term purpose is to put SDAE into the mainstream of stochastic calculus, developing as far as possible a theory similar to that of stochastic differential equations. In this first paper our aim is to investigate the solution of linear SDAE with constant coefficients and an additive white noise, that means A˙x(t) + Bx(t) = f(t)+Λξ(t), where ξis a white noise and A, B, Λ are constant matrices of appropriate dimensions. We shall first reduce the equation to the so-called Kronecker Canonical Form (KCF), which is easy to analyse, and from whose solution one can recover easily the solution to the original problem. Unfortunately, it is not possible to extend this approach to the case of linear SDAE with varying coefficients, just as happens in the deterministic case, where a bunch of different approaches have been proposed. Among these, the most promising in our opinion is that of Rabier and Rheinboldt [8]. Due to the simple structure of the equations considered here, it is not a hard task to establish the existence of a unique solution in the appropriate sense. However, as mentioned before, a DAE does not regularise the input f(t) in general. If white noise, or a similarly irregular noise is used as input, then the solution process to a SDAE will not be a usual stochastic process, defined as a random vector at every time t, but instead a “generalised process”, the stochastic analogous of a Schwartz generalised function. The paper is organised as follows: in the next section we shall provide a short introduction to linear DAE’s and to generalised processes. In the third section we shall define what we mean by a solution to a linear SDAE and in Section 4 we shall provide a sufficient condition for the existence of LINEAR STOCHASTIC DIFFERENTIAL-ALGEBRAIC EQUATIONS 3 a density of the law of the solution. In the final Section 5 we shall discuss a simple example arising in the modelling of electrical circuits. Superscripts in parentheses mean order of derivation. All function and vector norms throughout the paper will be L2norms. 2. Preliminaries on DAE and generalised processes In this section we briefly introduce two topics: the (deterministic) differential-algebraic equations and the generalised processes. An exhaustive introduction on the first topic can be found in Rabier and Rheinboldt [9], while the basic theory of generalised processes can be found in Dawson [1], Fernique [2], or Chapter 3 in Gel’fand and Vilenkin [3]. 2.1. Differential-Algebraic Equations. Consider an implicit autonomous ODE, (2.1) F(x, ˙x) = 0 , where F:= F(x, p) : Rn×n→Rnis a sufficiently smooth function. If the partial differential DpF(x, p) is invertible at every point (x0, p0), one can easily prove that the implicit ODE is locally reducible to an explicit ODE. If DpF(x0, p0) is not invertible, two cases are possible: either the total derivative DF (x0, p0) is onto Rnor it is not. In the first case, and assuming that the rank of DpF(x, p) is constant in a neighbourhood of (x0, p0), (2.1) is called a differential-algebraic equation, while in the remaining cases one speaks of an ODE with a singularity at (x0, p0). Alinear DAE is a system of the form (2.2) A(t) ˙x+B(t)x=f(t), t ≥0, where A(t), B(t)∈Rn×nand f(t)∈Rn. The matrix function A(t) is assumed to have a constant (non-full) rank for any tin the interval of interest. (Clearly, if A(t) has full rank for all tin an interval, then the DAE reduces locally to an ODE.) In the simplest case, when Aand Bdo not depend on t, we have a linear DAE with constant coefficients, and an extensive study of these problems has been developed. The theory starts with the definition of a regular matrix pencil: Definition 2.1. Given two matrices A, B ∈Rn×n, the matrix pencil (A, B) is the function λ7→ λA +B, for λ∈R. It is called a regular matrix pencil if det(λA +B)6= 0 for some λ. A classical result, due to Weierstrass and Kronecker, states that the matrices of a regular matrix pencil can be simultaneously transformed into a convenient canonical form, as stated in the following proposition (see e.g. Gripentrog and M¨arz [4] for the proof): 4 AURELI ALABERT AND MARCO FERRANTE Proposition 2.2. Given a regular matrix pencil (A, B), there exist nonsingular matrices Pand Qand integers 0≤d, q ≤n, with d+q=n, such that PAQ =µId0 0N¶and PBQ =µJ0 0Iq¶ where Id, Iqare identities of dimensions dand q,N=blockdiag(N1, . . . , Nr), with Nithe qi×qimatrix Ni=         0 1 0 . . . 0 0 0 1 · · · 0 . . ........ . . . . ....1 0· · · · · · · · · 0         , and Jis in Jordan canonical form. In what follows, we shall always assume that (H.1) (A, B) is a regular matrix pencil so that the proposition above applies. Notice that the matrix Nis nilpotent, with nilpotency index given by the dimension of its largest block. This nilpotency index of the matrix N in this canonical form is a characteristic of the matrix pencil and we shall call it the index of the equation (2.2). From a given DAE with constant coefficients Aand Bthat satisfy (H.1), multiplying from the left by Pand defining the new variables y=Q−1x, we get a new linear DAE with matrices PAQ and PBQ in place of Aand B, which can be easily solved (see Section 3). The regularity of the solution depends directly on the index of the equation. Remark 2.3. Without hypothesis (H.1), a linear DAE may possess an infinity of solutions or no solution at all, depending on the right-hand side. Consider for instance, µ0 0 1 1¶˙x(t) + µ1 1 0 0¶x(t) = µf1(t) f2(t)¶ with any fixed initial condition. 2.2. Generalised processes. Let D0be the space of distributions (generalised functions) on some open set U⊂R, that is, the dual of the space D=C∞ c(U) of smooth functions with compact support defined on U. A random distribution on U, defined on the probability space (Ω,F, P), is a measurable mapping X: (Ω,F)→(D0,B(D0)), where B(D0) denotes the Borel σ-field, relative to the strong dual topology (equivalently, the weak-? LINEAR STOCHASTIC DIFFERENTIAL-ALGEBRAIC EQUATIONS 5 topology). Denoting by hX(ω), φithe action of the distribution X(ω)∈ D0 on the test function φ∈ D, it holds that the mapping ω7→ hX(ω), φiis measurable from (Ω,F) into (R,B(R)), hence a real random variable hX, φion (Ω,F, P). The law of Xis determined by the law of the finite-dimensional vectors (hX, φ1i, . . . , hX, φdi), φi∈ D,d∈N. The sum of random distributions Xand Yon (Ω,F, P), defined in the obvious manner, is again a random distribution. The product of a real random variable αand a random distribution, defined by hαX, φi:= αhX, φi, is also a random distribution. The derivative of a random distribution, defined by hX0, φi:= −hX, φ0i, is again a random distribution. Given a random distribution X, the mapping X:D → L0(Ω) defined by φ7→ hX, φiis called a generalised stochastic process. This mapping is linear and continuous with the usual topologies in Dand in the space of all random variables L0(Ω). Notice that we can safely overload the meaning of the symbol X. The mean functional and the correlation functional of a random distribution are the deterministic distribution φ7→ E[hX, φi] and the bilinear form (φ, ψ)7→ E[hX, φihX, ψi], respectively, provided they exist. A simple example of random distribution is white noise ξ, characterised by the fact that hξ, φiis centred Gaussian, with correlation functional E[hξ, φihξ, ψi] = RRφ(s)ψ(s)ds. In particular, hξ, φiand hξ, ψiare independent if the supports of φand ψare disjoint. Whenever this property holds true for a process, we will say that it takes independent values on disjoint sets. In this paper we will use as the base set the half-line U=]0,+∞[. White noise on Ucoincides with the Wiener integral with respect to a Brownian motion W: Indeed, if φis a test function, then (2.3) hξ, φi=Z∞ 0 φ(s)dWs in the sense of equality in law. More precisely, the Wiener integral is defined as the extension to L2(R+) of white noise (see Kuo [6] for a construction of the Wiener integral as extension of white noise). Now, integrating by parts in (2.3), we can write hξ, φi=−Z∞ 0 W(s)˙ φ(s)ds =−hW, ˙ φi, so that ξis the derivative of the Brownian motion Was random distributions. A random distribution is Gaussian if every finite-dimensional projection is a Gaussian random vector. This is the case of white noise and Brownian motion. Further results on random distributions and generalised stochastic proces- 6 AURELI ALABERT AND MARCO FERRANTE ses can be found for instance in the classical papers by Dawson [1] and Fernique [2]. 3. The generalised process solution Consider the equation (3.1) A˙x+Bx =f+ Λξ , where Aand Bare n×nreal matrices, f:R+→Ris continuous, Λ is a n×m constant matrix, and ξis a m-dimensional white noise: ξ= (ξ1, . . . , ξm), with ξiindependent one-dimensional white noises. We first reduce the equation to Kronecker Canonical Form (KCF), see Proposition 2.2: There exist regular matrices Pand Qand integers 0 ≤ d, q ≤n, with d+q=n, such that PAQ =µId0 0N¶and PBQ =µJ0 0Iq¶ where Id, Iqare identities of dimensions dand q,N= blockdiag(N1, . . . , Nr), with Nithe qi×qimatrix Ni=         0 1 0 . . . 0 0 0 1 · · · 0 . . ........ . . . . ....1 0· · · · · · · · · 0         and Jis in Jordan canonical form. We can assume that the blocks of J corresponding to the eigenvalue 0, if any, are located in the lower rows of J. Multiplying equation (3.1) by Pfrom the left, and defining the new variables y=Q−1x, we get (3.2) µId0 0N¶˙y+µJ0 0Iq¶y=f+ Λξ , where for simplicity we use again fand Λ to denote the function Pf and the new “diffusion” matrix PΛ. System (3.2) can be split into two parts. The first one is a classical stochastic differential system of dimension d, and the second one is an “algebraic stochastic system” of dimension q. Denoting by uand vthe variables in the first and the second part respectively, by band cthe related partitioning of the vector function f, and by (Σ, R) the corresponding splitting of Λ into matrices of dimensions d×mand q×m, we can write the two LINEAR STOCHASTIC DIFFERENTIAL-ALGEBRAIC EQUATIONS 7 systems as follows: (3.3)    ˙u1 . . . ˙ud   +J   u1 . . . ud   =   b1 . . . bd   + Σ    ξ1 . . . ξm   , (3.4) N   ˙v1 . . . ˙vq   +   v1 . . . vq   =   c1 . . . cq   +R   ξ1 . . . ξm   . We refer to uas the differential variables and to vas the algebraic variables. For any initial vector u(t0) = u0∈Rd, and t0>0, the solution to (3.3) can be written, in the sense of equality in law, as u(t) = e−J(t−t0)·u0+Zt t0 eJ(s−t0)b(s)ds+Zt t0 eJ(s−t0)ΣdW(s)¸, t ∈]0,∞[, where W(t) is a m-dimensional standard Wiener process. It can also be expressed, if desired, as a generalised process: For φ∈ C∞ c(]0,∞[), hu, φi=Z∞ 0 e−J(t−t0)³u0+Zt t0 eJ(s−t0)b(s)ds´φ(t)dt(3.5) +Z∞ 0·Zt t0 e−J(t−s)ΣdW(s)¸φ(t)dt . On the other hand, system (3.4) consists of a number of decoupled blocks, which are easily solved by backward substitution. For instance, for the first block, (3.6) N1   ˙v1 . . . ˙vq1   +   v1 . . . vq1   =   c1 . . . cq1   +R1   ξ1 . . . ξm   , R1representing the sub-matrix of Rconsisting of the first q1rows, and c= (c1, . . . , cq1) the corresponding part of vector c, a recursive calculation gives the following generalised process solution, for a given φ∈ C∞ c(]0,+∞[): (3.7) hvj, φi= q1 X k=jDck+ m X `=1 (R1)k,`ξ`, φ(k−j)E, j = 1, . . . q1, which can be expressed with the help of a standard Wiener process Was (3.8) hvj, φi= q1 X k=jhZ∞ 0 ck(t)φ(k−j)(t)dt+ m X `=1 (R1)k,`Z∞ 0 φ(k−j)dW`(t)i, j = 1, . . . q1. 8 AURELI ALABERT AND MARCO FERRANTE The remaining blocks are treated in the same way to get the whole vector ysolving the KCF system (3.2). We can thus state the following result: Proposition 3.1. Under assumption (H.1) of Section 2, equation (3.1) admits a unique generalised process solution x=Qy, where Qis the matrix determined by the Kronecker canonical form and y= (u, v), with uthe solution to system (3.3), given by expression (3.5), and vthe solution to system (3.4). 4. The law of the solution In the previous section we have seen that the solution to a linear SDAE with regular pencil and additive white noise can be explicitly given as a functional of the input noise. From the modelling viewpoint, the law of the solution is the important output of the model. Using the explicit form above, one can try to investigate the features of the law in which one might be interested. To illustrate this point, we shall write down the joint law of the solution vector evaluated at a fixed arbitrary test function φand we shall investigate some absolute continuity properties we do not aim at a very general statement, but instead we want to show the sort of arguments that can be used in each particular instance). For notational simplicity, let us assume that the differential system has dimension d= 2 and that the algebraic part consists of a unique block of nilpotency index 2. We will also assume b= 0 and c= 0 in (3.3) and (3.4); the general case can be studied similarly. Thus, we are dealing with the following system (4.1)-(4.2): µ˙u1 ˙u2¶+Jµu1 u2¶=µσ11 . . . σ1m σ21 . . . σ2m¶   ξ1 . . . ξm   ,(4.1) µ0 1 0 0¶µ˙v1 ˙v2¶+µv1 v2¶=µρ11 . . . ρ1m ρ21 . . . ρ2m¶   ξ1 . . . ξm   .(4.2) Given a test function φ∈ C∞ c(]0,+∞[), the solution to the differential part (4.1), with initial condition u0= (u0 1, u0 2) at time t0>0, is hui, φi=Z∞ 0 (e−J(t−t0))iku0 kφ(t)dt+Z∞ 0·Zt t0 (e−J(t−s))ikσk`dW`(s)¸φ(t)dt LINEAR STOCHASTIC DIFFERENTIAL-ALGEBRAIC EQUATIONS 9 for i= 1,2, (see (3.5)), whereas the solution to the algebraic part (4.2) is given by hv1, φi=Z∞ 0 (ρ1`φ(t) + ρ2`˙ φ(t)) dW`(t)(4.3) hv2, φi=Z∞ 0 ρ2`φ(t)dW`(t),(4.4) (see (3.7)), with the convention of summation over repeated indices, that we will keep applying in the sequel without explicit mention. By standard computations, for any given test function φ, the random vector h(u1, u2, v1, v2), φihas a Gaussian distribution, with expectations E[hui, φi] = Z∞ 0 (e−J(t−t0))iku0 kφ(t)dt , i = 1,2 E[hv1, φi] = E[hv2, φi]=0 and covariances Cov[hui, φi,huj, φi] = m X `=1 Z∞ 0hZ∞ s φ(t) 2 X k=1 (e−J(t−s))ikσk` dti (4.5) ×hZ∞ s φ(t) 2 X k=1 (e−J(t−s))jkσk` dtids , i, j = 1,2 16 AURELI ALABERT AND MARCO FERRANTE Let us now reduce the equation to KCF. To simplify the exposition, we shall fix to 1 the values of τi,Rand Li. (A physically meaningful magnitude for Rand Liwould be of order 10−6for the first and of order 104for the latter. Nevertheless the structure of the problem does not change with different constants.) The matrices Pand Q, providing the desired reduction (see Proposition 2.2), are P=     1 2−1 21−1 0−1 1 1 1 1 0 0 −1 0 0 0     , Q =     −1 4−1 2−3 4−1 1 4−1 2−1 40 1 201 20 −1 201 20     . Indeed, multiplying (5.2) by Pfrom the left and defining y=Q−1x, we arrive to (5.3)     1000 0010 0000 0000    ˙y(t) +     1 20 0 0 0 1 0 0 0 0 1 0 0 0 0 1    y(t) =     −τ1τ2−τ3 −τ1−τ2−τ3 0 0 0 0 0 τ3      ξ1 ξ2 ξ3  , We see that the matrix Nof Section 3 has here two blocks: A single zero in the last position ( ˙y4) and a 2-nilpotent block affecting ˙y2and ˙y3. We have therefore an index 2 SDAE. From Propositions 4.1, 4.2 and Theorem 4.3, we can already say that, when applied to any test function φ6= 0, the variables y4,y2and y1, as well as the vectors (y1, y2)and (y1, y4), will be absolutely continuous, whereas y3degenerates to a point. In fact, in this case, we can of course solve completely the system: The differential part is the one-dimensional classical SDE (5.4) ˙y1+1 2y1=−τ1ξ1+τ2ξ2−τ3ξ3, and the algebraic part reads simply (5.5)    ˙y3+y2=−τ1ξ1−τ2ξ2−τ3ξ3 y3= 0 y4=τ3ξ3. The solution to (5.3) can thus be written as y1(t) = e−(t−t0)/2hy(t0) + Zt t0 e−(s−t0)/2(−τ1dW1+τ2dW2−τ3dW3)(s)i y2=−τ1ξ1−τ2ξ2−τ3ξ3 y3= 0 y4=τ3ξ3, LINEAR STOCHASTIC DIFFERENTIAL-ALGEBRAIC EQUATIONS 17 where W1,W2,W3are independent Wiener processes whose generalised derivatives are ξ1,ξ2and ξ3. Multiplying by the matrix Qwe finally obtain the value of the original variables: x1(t) = −x2(t) = −1 2y1(t) u1=−1 4y1−1 2y2−3 4y4 u2=1 4y1−1 2y2, with x1(t0) = −1 4y1(t0)a given intensity at time t0. It is clear that the current intensities, which have almost surely continuous paths, are much more regular than the voltage drops, which are only random distributions. References [1] Donald A. Dawson. Generalized stochastic integrals and equations. Trans. Amer. Math. Soc., 147:473–506, 1970. [2] Xavier Fernique. Processus lin´eaires, processus g´en´eralis´es. Ann. Inst. Fourier (Grenoble), 17(fasc. 1):1–92, 1967. [3] I. M. Gel0fand and N. Ya. Vilenkin. Generalized functions. Vol. 4. Academic Press, New York, 1964. [4] E. Griepentrog and R. M¨arz. Basic properties of some differential-algebraic equations. Z. Anal. Anwendungen, 8(1):25–41, 1989. [5] Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, Cambridge, 1990. [6] Hui Hsiung Kuo. Gaussian measures in Banach spaces. Springer-Verlag, Berlin, 1975. [7] David Nualart. The Malliavin calculus and related topics. Probability and its Applications (New York). Springer-Verlag, New York, 1995. [8] Patrick J. Rabier and Werner C. Rheinboldt. Classical and generalized solutions of time-dependent linear differential-algebraic equations. Linear Algebra Appl., 245:259–293, 1996. [9] Patrick J. Rabier and Werner C. Rheinboldt. Theoretical and numerical analysis of differential-algebraic equations. In Handbook of numerical analysis, Vol. VIII, pages 183–540. North-Holland, Amsterdam, 2002. [10] O. Schein. Stochastic differential-algebraic equations in circuit simulation. PhD thesis, Technische Universit¨at Darmstadt, 1999. [11] O. Schein and G. Denk. Numerical solution of stochastic differential-algebraic equations with applications to transient noise simulation of microelectronic circuits. J. Comput. Appl. Math., 100(1):77–92, 1998. [12] Renate Winkler. Stochastic differential algebraic equations of index 1 and applications in circuit simulation. J. Comput. Appl. Math., 157(2):477–505, 2003. 18 AURELI ALABERT AND MARCO FERRANTE Aureli Alabert Departament de Matem`atiques Universitat Aut`onoma de Barcelona 08193 Bellaterra, Catalonia E-mail:[email protected] Marco Ferrante Dipartimento di Matematica P. e A. Universit`a degli Studi di Padova via Belzoni 7, 35131 Padova, Italy E-mail:[email protected]