Invariant Galerkin Ansatz Spaces and Davison-Maki Methods for the Numerical Solution of Differential Riccati Equations
IInvariant Galerkin Ansatz Spaces and
Davison-Maki Methods for the NumericalSolution of Differential Riccati Equations
Maximilian Behr Peter Benner Jan HeilandOctober 30, 2019
Abstract
The differential Riccati equation appears in different fields of applied mathematics like controland system theory. Recently Galerkin methods based on Krylov subspaces were developed for theautonomous differential Riccati equation. These methods overcome the prohibitively large storagerequirements and computational costs of the numerical solution. In view of memory efficientapproximation, we review and extend known solution formulas and identify invariant subspacesfor a possibly low-dimensional solution representation. Based on these theoretical findings, wepropose a Galerkin projection onto a space related to a low-rank approximation of the algebraicRiccati equation. For the numerical implementation, we provide an alternative interpretation ofthe modified
Davison-Maki method via the transformed flow of the differential Riccati equation,which enables us to rule out known stability issues of the method in combination with theproposed projection scheme. We present numerical experiments for large-scale autonomousdifferential Riccati equations and compare our approach with high-order splitting schemes.
Contents
1. Introduction 22. Preliminaries 33. Algebraic and Differential Riccati Equations 34. Radon’s Lemma 5
5. Galerkin Approach for Large-Scale Differential Riccati Equations 14
6. Numerical Experiments 22 a r X i v : . [ m a t h . NA ] O c t . Conclusion 26A. Numerical Results for Galerkin Approach 30B. Numerical Results for Splitting Schemes 32
1. Introduction
In this paper we consider the autonomous differential Riccati equation˙ X ( t ) = A T X ( t ) + X ( t ) A − X ( t ) BB T X ( t ) + C T C,X (0) = X . The equation plays an important role in model order reduction, optimal control, differentialgames and stability analysis [1, 7, 15, 17, 32, 38, 40]. We focus in this work on the large-scale case.In this setting, the numerical approximation of X comes with high memory requirements andhigh computational costs. Just the storage of the solution at the relevant time instances wouldscale with N t n , where n is the dimension of the problem and N t is the number of time steps.The approach of first discretizing in time and then focusing on efficient approximation of theresulting algebraic equations has been the main course of research on this problem setup, see,e.g., [13–15, 27, 33, 35, 36, 41, 42, 52, 54]. In all these approaches, the approximation of at least onelarge-scale algebraic equation has to be solved and stored for every time step so that the memorydemands still scale with N t n . Conceptually, it seems more beneficial for the autonomous differentialRiccati equation to first reduce the problem dimensions to, say, k (cid:28) n and then approach thereduced equation as this leads to storage requirements in the order of N t k . In this respect, Krylovsubspace methods have been proposed [4, 23–26, 31, 34] that generate a trial space for the numericalsolution using an Arnoldi method. The resulting Galerkin projected system is of lower order andcan be solved with low memory demand and with various methods that exist for differential Riccatiequations of small or moderate size.In this work, we develop a Galerkin approach, where the trial space is based on the numericalsolution of the algebraic Riccati equation. This extends the concepts of our previous work on anumerical scheme for differential Lyapunov equations [10].The paper is organized as follows. In Section 3 we introduce the algebraic and differential Riccatiequations and review the relevant fundamental properties about their solutions. In Section 4we review Radon’s Lemma and work out its implication that the differential Riccati equation isconnected to a flow on a Grassmanian manifold. Moreover, in Section 4.2, we apply
Radon’s Lemma to obtain solution formulas for the differential Riccati equation based on the solution of the algebraicRiccati equation that we will use to explain and illustrate the major source of numerical instabilitiesof the
Davison-Maki method for the numerical solution of the differential Riccati equation; seeSection 4.3 Then we will use the connection to the Grassmanian manifold to derive the modifiedDavison-Maki method in a way that overcomes these instabilities. In Section 5, we develop aGalerkin approach for the solution of the differential Riccati equation in the matrix exponentialrepresentation that results from
Radon’s Lemma . We combine the monotonicity of the solution ofthe differential and relevant properties of the solution of the algebraic Riccati equation to definea suitable and numerically computable trial space for the approximation of the solution of thedifferential Riccati equation. We propose to solve the resulting Galerkin system with the modifiedDavison-Maki method . Numerical results are presented in Section 6 and Appendices A and B.2 . Preliminaries
In this section we set the notation and review some basic results from linear algebra. The identitymatrix and zero matrix of size n are written by I n and 0. The image or column space of a matrix A ∈ R n × m is denoted by range ( A ), and its kernel or null space by ker( A ). The 1–norm, 2–norm,Frobenius norm and Frobenius inner product are denoted by k·k , k·k , k·k F and h· , ·i F , respectively.The spectrum of a quadratic matrix A is denoted by Λ( A ). Generally, the spectrum is a subsetof C . A matrix is called stable, if its spectrum is contained in the left open complex half plane C − , i.e. Λ( A ) ⊆ C − . If A is real and symmetric, all eigenvalues are real and λ ↓ k ( A ) represents the k –largest eigenvalue. Therefore, λ ↓ ( A ) ≥ λ ↓ ( A ) ≥ · · · ≥ λ ↓ n ( A ) are the eigenvalues of A orderedin a non-decreasing fashion. The Loewner partial ordering on the set of real symmetric matricesis defined by A (cid:52) B , which means B − A is positive semidefinite, [28, Ch. 7.7]. The orthogonalcomplement of a linear subspace U ⊆ R n is denoted by U ⊥ ⊆ R n . For A ∈ R n × n and B ∈ R n × b ,the image of the Krylov matrix generated by A and B is denoted by K ( A, B ) := range (cid:16)h
B, AB, . . . , A n − B i(cid:17) ⊆ R n . The linear space K ( A, B ) is A –invariant.
3. Algebraic and Differential Riccati Equations
In this section we introduce the algebraic and differential Riccati equation (ARE/DRE) and thealgebraic Lyapunov equation (ALE).Consider
A, X ∈ R n × n and C ∈ R c × n and B ∈ R n × b . Throughout this paper, we assume that X is a symmetric positive semidefinite matrix and consider the DRE˙ X ( t ) = R ( X ( t )) := A T X ( t ) + X ( t ) A − X ( t ) BB T X ( t ) + C T C, (1a) X (0) = X . (1b)Stationary points of (1a) are solutions of the corresponding ARE0 = R ( X ) = A T X + XA − XBB T X + C T C. (2)The linear version ( B = 0) of the ARE is the ALE0 = A T X + XA + C T C. (3)We review some fundamental results about existence, uniqueness and properties of the solution ofthe DRE (1), ARE (2) and the ALE (3). Theorem 3.1 (Existence and Uniqueness of Solutions to the ALE (3), [1, Thm. 1.1.3, 1.1.7]) . If Λ( A ) ∩ Λ( − A ) = ∅ , then the ALE (3) admits a unique solution X L ∈ R n × n . The solution X L issymmetric. If A is stable, then X L is symmetric positive semidefinite and given by X L = ∞ Z e tA T C T Ce tA d t. (4) Theorem 3.2 (Existence and Uniqueness of Solutions for the ARE (2), [1, Lem. 2.4.1, Cor.2.4.3], [32, Ch. 10]) . Let ( A, B ) be stabilizable and ( A, C ) be detectable, then the ARE (2) has a unique stabilizing solution ∞ ∈ R n × n . This means R ( X ∞ ) = 0 and Λ( A − BB T X ∞ ) ⊆ C − . Moreover X ∞ is symmetricpositive semidefinite and there is no other symmetric positive semidefinite solution of the ARE (2) . Theorem 3.3 (Range of the Solution of the ARE (2), [9, Thm. 3.2]) . Let ( A, B ) be stabilizable, ( A, C ) be detectable and X ∞ ∈ R n × n be the unique stabilizing solution ofthe ARE (2) . Then the following relation holds: range ( X ∞ ) = K (cid:16) A T , C T (cid:17) . The inclusion K (cid:16) A T , C T (cid:17) ⊆ range ( X ∞ ) in Theorem 3.3 is actually true for each symmetric solutionof the ARE (2), cf. [1, Lemma 2.4.9]. In [2, Ch. 3.3] a Kalman decomposition is used to showthat rank ( X ∞ ) = dim (cid:16) K (cid:16) A T , C T (cid:17)(cid:17) . A connection between the space K (cid:16) A T , C T (cid:17) and a certainKrylov subspace generated by the associated Hamiltonian matrix, which can be used for numericalapproximation of the solution of the ARE (2), was presented in [11, Thm. 10].Typically, solutions of quadratic differential equations like the DRE (1) exhibit a finite-time escapephenomena. By means of comparison arguments and the fact that − BB T is negative semidefiniteone can show that the solution exists for all t ≥
0. With additional assumptions, the solutionconverges monotonically to the unique solution of ARE (2) and is, thus, bounded.
Theorem 3.4 (Existence and Uniqueness of Solutions of the DRE (1), [1, Thm. 4.1.6, 4.1.8], [32, Ch.10]) . The
DRE (1) has a unique solution X : [0 , ∞ ) → R n × n . The solution X has the following properties:• X ( t ) is symmetric positive semidefinite for all t ≥ .• If ˙ X (0) = R ( X ) (cid:60) , then t X ( t ) is monotonically non-decreasing on [0 , ∞ ) , i.e. X ( t ) (cid:52) X ( t ) for all t , t such that ≤ t ≤ t . Theorem 3.5 (Invariant Subspace of the Solution of the DRE (1), cp. [9, Thm. 3.1]) . Let the columns of Q ∈ R n × p span an orthonormal basis of K (cid:16) A T , C T (cid:17) and define the linear space Q := n QY Q T | Y ∈ R p × p o ⊆ R n × n or Q := { } ⊆ R n × n , if C = 0 . Then the following holds: X ( t ) ∈ Q for all t ≥ , where X is the unique solution of the DRE (1) with X = 0 . With this relation, one can readily confirm that the solution of the DRE (1) evolves in an invariantsubspace of R n × n .For numerical approximations of the solutions of large-scale ALEs, AREs and DREs, one typicallyseeks for low-rank approximations, i.e. a q (cid:28) n so that the relation in Theorem 3.5 is still valid up toa given tolerance, to avoid overly demanding memory requirements. Therefore, the relevant literaturefeatures numerous contributions which study the decay rate of λ ↓ k ( X ) or λ ↓ k ( X ) λ ↓ ( X ) for increasing k ; see,e.g., [5, 6, 8, 21, 22, 44, 45, 51] on the eigenvalue decay of the solution of the ALE and [11, 44, 54] forresults on the ARE and DRE.For the autonomous DRE (1), one can derive estimates based on the monotonicity. Assume that R ( X ) (cid:60)
0, then by Theorem 3.4 the function t X ( t ) is monotonically non-decreasing on [0 , ∞ ),where X is the unique solution of the DRE (1). A direct consequence of the Courant-Fischer-Weylmin-max principle [28, Cor. 7.7.4] implies that t λ ↓ k ( X ( t )) is also monotonically non-decreasingon [0 , ∞ ). Therefore the number of eigenvalues of X ( t ) greater than or equal to a given threshold ε > xample 3.1 (Eigenvalue Decay) . We illustrate this by an example in
Figure 1 . We have chosen X = 0 , C = h , . . . , i = B T and A to be tridiagonal with entries , − , − on the subdiagonal, diagonal and superdiagonal, respectively.The matrices are of size n = 100 and the DRE was solved numerically to a high precision on the timeinterval [0 , . For this we have used the variable-precision arithmetic vpa of MATLAB ® with significant digits and Algorithm 2 with step size h = 2 − . The eigenvalues of X ( t ) arearranged in a non-increasing order and plotted for t ∈ { . , , . . . , } . The functions t λ ↓ k ( X ( t )) are highlighted in red for k ∈ { , , , , } . All eigenvalues below − were truncated from Figure 1 . The shadowed red plane is drawn at the level · − , which is approximately machineprecision in double arithmetic. − − − − − k t ∈ [0 , λ ↓ k ( X ( t )) Fig. 1. Eigenvalues λ ↓ k ( X ( t )) of the numerical solution of DRE (1).
4. Radon’s Lemma
In this section, we consider the non-symmetric differential Riccati equation abbreviated by NDREas a generalization of the DRE. We will make heavy use of
Radon’s Lemma that shows that theNDRE is locally equivalent to a linear differential equation of twice the size. Vice versa, the solutionof the NDRE defines the solution of an associated linear system.
Radon’s Lemma (Thm. 4.1) has several consequences. In Section 4.1, we review the fact that thesolution of the NDRE induces a flow on the Grassmanian manifold. This flow has a simpler structureas it is based on a matrix exponential. In Section 4.2 we show how solution formulas can beobtained by applying suitable linear transformations, which decouple the linear differential equation.Then, in view of numerical approximation, we review the
Davison-Maki method and the modifiedDavison-Maki method in Section 4.3. We use the solution formula from Section 4.2 to explain, whythe
Davison-Maki method applied to the DRE usually suffers from numerical instabilities and showthat an exploitation of the structure of the transformed flow on the Grassmanian manifold leads toa suitable modification of the
Davison-Maki method . Theorem 4.1 (Radon’s Lemma, [1, Thm. 3.1.1]) . Let M ∈ R n × n , M ∈ R n × m , M , M ∈ R m × n , M ∈ R m × m and I ⊆ R be an open intervalsuch that ∈ I . We consider the NDRE˙ W ( t ) = M W ( t ) − W ( t ) M − W ( t ) M ( t ) W ( t ) + M , (5a)5 (0) = M . (5b) The following holds:1. Let W : I → R m × n be the solution of (5) and U : I → R n × n be the solution of the linear initialvalue problem ˙ U ( t ) = ( M + M W ( t )) U ( t ) , U (0) = I n . (6) Moreover let V ( t ) := W ( t ) U ( t ) . Then U : I → R n × n and V : I → R m × n define the solution of " ˙ U ( t )˙ V ( t ) = M " U ( t ) V ( t ) := " M M M M U ( t ) V ( t ) , " U (0) V (0) = " I n M . (7)
2. If " UV : I → R ( n + m ) × n is a solution of (7) and the matrix U ( t ) is nonsingular for all t ∈ I ,then W : I → R m × n , W ( t ) = V ( t ) U ( t ) − is a solution of (5) . Radon’s Lemma (Thm. 4.1) also holds for time-dependent continuous matrix valued functionsas coefficients. Note that, usually, the solution of the NDRE (5) has finite time escape, whilethe solution of system (7) exists for all t ∈ R . However, one can consider the solution W of theNDRE (5) on the interval of existence. As the function U is a solution of the linear initial valueproblem (6) and U (0) = I n is nonsingular, the determinant of U ( t ) can not vanish on the interval I .It follows that the matrix U ( t ) is nonsingular for all t ∈ I , c.f. [57, §15]. Therefore as long as thesolution of the NDRE (5) exists, it can be recovered from the solution of system (7). In this section we review the fact that the solution of the NDRE (5) is locally equivalent to a flowon the Grassmanian manifold. This connection was first observed in [49] and the corresponding flowwas further studied in [39, 50]. The content of this subsection is a summary of [50, §2]. One mainobservation from Radon’s Lemma (Thm. 4.1) is that the solution W of the NDRE (5) depends onlyon the linear space spanned by U ( t ) and V ( t ). This can be seen by the following arguments. Let " UV : I → R ( n + m ) × n be a solution of (7) and t ∈ I . Moreover assume that ˜ U ∈ R n × n , ˜ V ∈ R m × n are such that range " ˜ U ˜ V = range " U ( t ) V ( t ) . The linear spaces are equal, if and only if there is a nonsingular matrix T ∈ R n × n such that " ˜ U ˜ V = " U ( t ) V ( t ) T. Since U ( t ) is nonsingular, we have˜ V ˜ U − = V ( t ) T T − U ( t ) − = V ( t ) U ( t ) − = W ( t ) . " U ( t ) V ( t ) ⊆ R n + m that defines the solution W ( t ),rather than the chosen basis " U ( t ) V ( t ) to represent the space. Sincerange " U ( t ) V ( t ) = range e tM " I n M = e tM range " I n M , and the (nonsingular) matrix exponential is applied to an n -dimensional subspace range " I n M ,we obtain a time-dependent family of n -dimensional subspaces of R n + m . The Grassmanian manifold G n (cid:0) R n + m (cid:1) consists of all n -dimensional subspaces of R n + m . Therefore the flow associated to theNDRE (5) on G n (cid:0) R n + m (cid:1) is given by ϕ : R × G n (cid:16) R n + m (cid:17) → G n (cid:16) R n + m (cid:17) , ( t, S ) e tM S. The flow exists for all t ∈ R and has the flow properties ϕ (0 , S ) = S and ϕ ( t , ϕ ( S, t )) = ϕ ( t + t , S )for all S ∈ G n (cid:0) R n + m (cid:1) and t , t ∈ R .In addition it holds that U ( t ) is nonsingular as long as W exists. This motivates us to consider theset of all graph subspaces of G n (cid:0) R n + m (cid:1) G n (cid:16) R n + m (cid:17) := ( range " UV U ∈ R n × n , V ∈ R m × n , det U = 0 ) ⊆ G n (cid:16) R n + m (cid:17) , together with the function ψ : G n (cid:16) R n + m (cid:17) → R m × n , range " UV V U − . The function ψ is well defined, as it does not depend on the basis of the graph subspace. Thus, wehave that W ( t ) = ψ range " U ( t ) V ( t ) = ψ ϕ t, range " I n M , and ψ − ( W ( t )) = range " I n W ( t ) = range " I n V ( t ) U ( t ) − = range " U ( t ) V ( t ) = ϕ t, range " I n M , as long as the solution W exists. Therefore the solution of the NDRE (5) induces a flow on theGrassmanian manifold. The solution W can be recovered from the flow by using ψ , and, vice versa,the flow can be obtained from the solution of the NDRE (5) using ψ − . Radon’s Lemma (Thm. 4.1) enables a certain solution representations for the DRE (1): Theorem 3.4ensures that the DRE (1) has a unique solution for t ≥
0. By Radon’s Lemma (Thm. 4.1) we have7hat U ( t ) is nonsingular for all t ≥ H := " A − BB T − C T C − A T ∈ R n × n be the Hamiltonian matrix corresponding to the DRE (1).The matrices U ( t ) and V ( t ) are determined by the linear initial value problem " ˙ U ( t )˙ V ( t ) = − H " U ( t ) V ( t ) , " U (0) V (0) = " I n X . (8)We obtain " U ( t ) V ( t ) = e − tH " I n X . The strategy is to decompose the Hamiltonian matrix H , such that (8) decouples. Theorem 4.2 (Solution representation I for DRE (1), [48]) . Let X ∈ R n × n be any solution of the ARE (2) . Then the solution of the
DRE (1) for t ≥ is givenby X ( t ) = X − e t ( A − BB T X T ) T ˜ X I n − t Z e s ( A − BB T X ) BB T e s ( A − BB T X T ) T d s ˜ X − e t ( A − BB T X ) , ˜ X := X − X . Proof.
We use T := " I n X I n and apply a similarity transformation to H , T − HT = " I n − X I n A − BB T − C T C − A T I n X I n = " A − BB T X − BB T − ( A − BB T X T ) T =: ˜ H. This gives " U ( t ) V ( t ) = e − tH " I n X = e − tT ˜ HT − " I n X = T e − t ˜ H T − " I n X = T e − t ˜ H " I n X − X =: T " ˜ U ( t )˜ V ( t ) . Clearly ˜ U and ˜ V are determined by the solution of the initial value problem " ˙˜ U ( t )˙˜ V ( t ) = − ˜ H " ˜ U ( t )˜ V ( t ) = " − ( A − BB T X ) BB T A − BB T X T ) T ˜ U ( t )˜ V ( t ) , " ˜ U (0)˜ V (0) = " I n X − X . By using the variation of constants formula [57, §18] we obtain that ˜ U and ˜ V are given by˜ V ( t ) = − e t ( A − BB T X T ) T ( X − X ) , ˜ U ( t ) = e − t ( A − BB T X ) + t Z e − ( t − s )( A − BB T X ) BB T ˜ V ( s )d s e − t ( A − BB T X ) I n − t Z e s ( A − BB T X ) BB T e s ( A − BB T X T ) T d s ( X − X ) . Since ˜ U ( t ) = U ( t ) is nonsingular for all t ≥ t ≥
0. Finally we obtain V ( t ) = X ˜ U ( t ) + ˜ V ( t ) ,X ( t ) = V ( t ) U ( t ) − = X + ˜ V ( t ) ˜ U ( t ) − . The formula was presented in [48] without proof. Since the existence of the involved inverse is nottrivially established, we provide a proof.
Theorem 4.3 (Solution Representation II for DRE (1), [18, Thm. 1], [47]) . Let ( A, B ) be stabilizable and ( A, C ) be detectable and X ∞ ∈ R n × n be the unique symmetric positivedefinite stabilizing solution of the ARE (2) . Moreover let ˆ A := A − BB T X ∞ and X L ∈ R n × n be theunique symmetric positive semidefinite solution of the Lyapunov equation ˆ AX L + X L ˆ A T + BB T = 0 . (9) Then the solution of the
DRE (1) for t ≥ is given by X ( t ) = X ∞ − e t ˆ A T ( X ∞ − X ) (cid:16) I n − ( X L − e t ˆ A X L e t ˆ A T )( X ∞ − X ) (cid:17) − e t ˆ A . Proof.
Similar to the proof of Theorem 4.2 we use similarity transformations to decompose theHamiltonian matrix H . This is also known as a Riccati-Lyapunov transformation [1, Ch. 3.1.1.].We obtain T : = " I n X ∞ I n , T − HT = " ˆ A − BB T − ˆ A T =: ˜ H, ˜ T : = " I n − X L I n , ˜ T − ˜ H ˜ T = " ˆ A − ˆ A T =: ˆ H. We thus get " U ( t ) V ( t ) = e − tH " I n X = e − t ( T ˜ T ) ˆ H ( T ˜ T ) − " I n X = ( T ˜ T ) e − t ˆ H ( T ˜ T ) − " I n X = " I n − X L X ∞ I n − X ∞ X L e − t ˆ A e t ˆ A T I n − X L X ∞ X L − X ∞ I n I n X = " e − t ˆ A ( I n − X L X ∞ ) + X L e t ˆ A X ∞ e − t ˆ A X L − X L e t ˆ A T X ∞ e − t ˆ A ( I n − X L X ∞ ) − ( I n − X ∞ X L ) e t ˆ A X ∞ X ∞ e − t ˆ A X L + ( I n − X ∞ X L ) e t ˆ A I n X = " e − t ˆ A ( I n − X L ( X ∞ − X )) + X L e t ˆ A T ( X ∞ − X ) X ∞ e − t ˆ A ( I n − X L ( X ∞ − X )) − ( X ∞ X L + I n ) e t ˆ A T ( X ∞ − X ) . (10)Now observe that U ( t ) = e − t ˆ A (cid:16) I n − ( X L − e t ˆ A X L e t ˆ A T )( X ∞ − X ) (cid:17) , ( t ) = X ∞ e − t ˆ A (cid:16) I n − ( X L − e t ˆ A X L e t ˆ A T )( X ∞ − X ) (cid:17) − e t ˆ A T ( X ∞ − X )= X ∞ U ( t ) − e t ˆ A T ( X ∞ − X ) , therefore X ( t ) = V ( t ) U ( t ) − = X ∞ − e t ˆ A T ( X ∞ − X ) (cid:16) I n − (cid:16) X L − e t ˆ A X L e t ˆ A T (cid:17) ( X ∞ − X ) (cid:17) − e t ˆ A . In [3, Ch. 15.4] one can find another solution formula, which holds under more restrictive assumptions.A solution formula based on the Jordan canonical form is given in [1, Thm. 3.2.1].
The
Davison-Maki method for the NDRE (5) was proposed in [20]. The method is based onfirst computing the matrix exponential e hM for a given step size h >
0. According to Radon’sLemma (Thm. 4.1) we have that " U ( h ) V ( h ) = e hM " I n M , W ( h ) = V ( h ) U ( h ) − . The next step is then to make use of the semigroup property of the matrix exponential " U (2 h ) V (2 h ) = e hM " I n M = (cid:16) e hM (cid:17) " I n M , W (2 h ) = V (2 h ) U (2 h ) − . For the further steps we obtain " U ( kh ) V ( kh ) = (cid:16) e hM (cid:17) k " I n M , W ( kh ) = V ( kh ) U ( kh ) − . (11)Another variant of the Davison-Maki method updates U and V instead of the matrix exponential.The variant follows from " U ( kh ) V ( kh ) = e khM " I n M = e hM e ( k − hM " I n M = e hM " U (( k − h ) V (( k − h ) . (12)Both variants of the method are given in Algorithm 1.10 lgorithm 1 Davison-Maki method for the NDRE (5) [20, 30]
Assumption:
The NDRE (5) has a solution W : [0 , t f ) → R m × n . Input:
Real matrices M and M ij as in Theorem 4.1, step size h > t f > Output:
Matrices W k , such that W ( kh ) = W k for k ∈ N and kh < t f . W = M ; k = 1;% Compute matrix exponential e.g. by a scaling and squaring method: Θ h = exp h " M M M M ; Variant with matrix exponential update: Θ = Θ h ; while kh < t f do Partition n m (cid:20) (cid:21) n Θ Θ m Θ Θ = Θ; U dm = Θ + Θ M ; V dm = Θ + Θ M ; W k = V dm U dm − ; Θ = ΘΘ h ; k = k + 1; end while Variant with updating U and V : U dm = I n ; V dm = M ; Partition n m (cid:20) (cid:21) n Θ Θ m Θ Θ = Θ; while kh < t f do U dm = Θ U dm + Θ V dm ; V dm = Θ U dm + Θ V dm ; W k = V dm U − dm ; k = k + 1; end while When the
Davison-Maki method (Alg. 1) is applied to the DRE (1), usually numerical instabilitiesoccur which are due to the fact that each block of e − tH as well as U ( t ) and V ( t ) contains thematrix e − t ˆ A , cp. equation (10). Since ˆ A = A − BB T X ∞ is stable, the matrix exponential of − t ˆ A exhibits exponential growth which becomes problematic for large t . The occurrence of thesenumerical problems with the Davison-Maki method (Alg. 1) was also pointed out in [19, 30, 37, 56].Another reason is that the spectrum of a real Hamiltonian matrix comes in quadruples, thatis Λ( H ) = { λ , . . . , λ n , − λ , . . . , − λ n } with Re( λ i ) ≤
0. Therefore, usually, the spectrum of theHamiltonian contains eigenvalues with positive real part and, thus, also it’s matrix exponentialgrows [43, Prop. 2.3.1].A suitable modification of the
Davison-Maki method (Alg. 1) was proposed in [30], but the modifiedmethod originates back to [29, p. 9]. By Radon’s Lemma (Thm. 4.1), as laid out in Section 4.1, we11ave the identity W ( kh ) = ψ range " U ( kh ) V ( kh ) = ψ e khM range " I n M = ψ e hM range " U (( k − h ) V (( k − h ) = ψ e hM range " I n W (( k − h ) = ψ range e hM " I n W (( k − h ) . Therefore the iteration for the modified Davison-Maki method is given by " ˜ U ˜ V := e hM " I n W (( k − h ) , W ( kh ) = ˜ V ˜ U − . (13)The modified Davison-Maki method is given in Algorithm 2. Algorithm 2
Modified Davison-Maki method for the NDRE (5) [29, 30]
Assumption:
The NDRE (5) has a solution W : [0 , t f ) → R m × n . Input:
Real matrices M and M ij as in Theorem 4.1, step size h >
0, final time t f > tol exp > Output:
Matrices W k , such that W ( kh ) = W k for k ∈ N and kh < t f . W = M ; k = 1;% Compute matrix exponential e.g. by a scaling and squaring method: Θ = exp h " M M M M ;% Check the norm of the matrix exponential: if k Θ k > tol exp then return Error(”1-Norm of the matrix exponential is too large, decrease the step size h “). end if Partition n m (cid:20) (cid:21) n Θ Θ m Θ Θ = Θ; while kh < t f do U mod dm = Θ + Θ W k − ; V mod dm = Θ + Θ W k − ; W k = V mod dm U − mod dm ; k = k + 1; end while A decrease of the step size h >
0, does not improve the accuracy in general, because the iteration isexact. The accuracy is determined by the accuracy of the matrix exponential computation and thematrix inversion. The step size cannot be chosen arbitrary large as the matrix exponential maybecome too large in norm. In practice we suggest to compute the norm of the matrix exponentialbefore the iteration starts. If the norm is too large, then the step size has to be decreased. In the k -th iteration of Algorithm 2 we have " U mod dm V mod dm = e hM " I n W (( k − h ) = " Θ Θ Θ Θ I n W (( k − h ) , k U mod dm k ≤ k Θ k + k Θ k k W (( k − h ) k , k V mod dm k ≤ k Θ k + k Θ k k W (( k − h ) k . For small step sizes of h > e hM ≈ I n + m + hM and Θ ≈ I n + hM , Θ ≈ hM ,Θ ≈ hM and Θ ≈ I m + hM . Therefore for small enough step size and moderate norm ofthe solution k W ( t ) k , the norm of the iterates cannot grow heavily in contrast to Algorithm 1. Ifthe norm of the iterates becomes too large during iteration, the step size should be decreased.Assume that the matrix exponential in line 3 of Algorithm 2 was approximated by using the scalingand squaring method, then the intermediates of the squaring phase can be used and the matrixexponential needs not be recomputed from scratch. Example 4.1 (Exponential Growth Davison-Maki method) . We applied the Davison-Maki method (Alg. 1) with step size h = 2 − to a DRE with the samematrices
A, B, C and X as for Example 3.1 . We plot the -norm of the iterates U dm and V dm aswell as the -norm condition number of U dm on the interval [0 , . The plot shows that all quantitiesgrow exponentially over time. Therefore, eventually, either a floating point overflow will occur orthe matrix inversion ceases to be executed accurately. Figure 3 shows the same quantities for theiterates U mod dm and V mod dm of the modified Davison-Maki method (Alg. 2) . . . . . t ∈ [0 , c o nd ( U dm ) (a) Condition number of U dm . . . . . t ∈ [0 , k U dm k (b) 2-norm U dm . . . . . t ∈ [0 , k V dm k (c) 2-norm V dm . Fig. 2. Davison-Maki method Algorithm 1 and the growth of U dm and V dm . . . . . − t ∈ [0 , c o nd ( U moddm ) (a) Condition number of U mod dm . . . . . − t ∈ [0 , k U moddm k (b) 2-norm U mod dm . . . . . − t ∈ [0 , k V moddm k (c) 2-norm V mod dm . Fig. 3. Modified Davison-Maki method Algorithm 2 and the growth of U mod dm and V mod dm .If a symmetric solution is expected, then line 11 in Algorithm 2 should be altered with W k = (cid:16) W k + W Tk (cid:17) , because due to numerical errors the symmetry will be lost after some iterations.Any computational efficient norm can also be used for the matrix exponential in Algorithm 2line 4. The modified Davison-Maki method is also more efficient than the Davison-Maki method in13oth variants, because less matrix-matrix products are needed by time step, compare Algorithm 2line 8-13 with Algorithm 1 line 5-12 and line 16-21.The computational cost apart from matrix exponential computation grows linearly with the timestep size h , compare Algorithm 2 line 8-13.The intermediates U dm , V dm from Algorithm 2 and U mod dm , V mod dm from Algorithm 1 are usuallydifferent. The next lemma shows the connection. Lemma 4.1.
In the k -th iteration of Algorithm 1 and
Algorithm 2 , the iterates U dm , V dm and U mod dm , V mod dm fulfill range " U dm V dm = range " U mod dm V mod dm . Proof.
From equation (11) it followsrange " U dm V dm = range e hkM " I n M . Equation (13) givesrange " U mod dm V mod dm = range e hM " I n W (( k − h ) = range e hM " I n V (( k − h ) U (( k − h ) − = range e hM " U (( k − h ) V (( k − h ) = range e khM " I n M .
5. Galerkin Approach for Large-Scale Differential Riccati Equations
In this section we develop a feasible numerical approach for large-scale differential Riccati equations.We consider the DRE (1) and assume that X = 0. We develop the Galerkin approach based ontwo theoretical considerations. First we use the solution formula of Theorem 4.3. We show thatthe range of the solution X ∞ of the ARE is invariant under the action of the closed-loop matrix A − BB T X ∞ . It follows then that the action of the matrix exponential of the closed-loop matrix on X ∞ has the same property. This makes the approach consistent in the sense that the evolution doesnot leave the ansatz space and provides reasoning that the consistency error made by a numericalapproximation to these subspaces can be made arbitrarily small. Moreover, this invariance propertyallows for a straight-forward low-dimensional approximation of the matrix exponential. After thatwe show that, for our proposed choice of a Galerkin basis, a quick decay of the eigenvalues of thesolution of the ARE implies a decent approximation of the solution X ( t ) of the DRE.The result is a low-dimensional solution space with an accessible formula for the relevant matrixexponential so that we can use the modified Davison-Maki (Algorithm 2) for an efficient solution ofthe projected Galerkin system. First we prove that the range space of the solution X ∞ of the ARE is invariant under the action ofthe transposed closed-loop matrix (cid:16) A − BB T X ∞ (cid:17) T .14 emma 5.1. Let ( A, B ) be stabilizable, ( A, C ) be detectable and X ∞ ∈ R n × n be the unique stabilizing solution ofthe ARE (2) . Then range( X ∞ ) is ( A − BB T X ∞ ) T –invariant.Proof. We can assume that X ∞ = 0. Let the columns of Q ∞ ∈ R n × p be an orthonormal basis forrange( X ∞ ). Then Q ∞ Q T ∞ is the orthogonal projection onto range( X ∞ ). We obtain Q ∞ Q T ∞ X ∞ = X ∞ . By Theorem 3.3, the columns of Q ∞ are also an orthonormal basis for K (cid:16) A T , C T (cid:17) . The space K (cid:16) A T , C T (cid:17) is A T –invariant. We obtain A T Q ∞ = Q ∞ Q T ∞ A T Q ∞ . Finally, we have ( A − BB T X ∞ ) T Q ∞ = Q ∞ Q T ∞ A T Q ∞ − Q ∞ Q T ∞ X ∞ BB T Q ∞ = Q ∞ (cid:16) Q T ∞ A T Q ∞ − Q T ∞ X ∞ BB T Q ∞ (cid:17) . This means range( X ∞ ) is ( A − BB T X ∞ ) T –invariant.According to Theorem 4.3 the solution of the DRE (1) is for t ≥ X ( t ) = X ∞ − e t ˆ A T X ∞ (cid:16) I n − (cid:16) X L − e t ˆ A X L e t ˆ A T (cid:17) X ∞ (cid:17) − e t ˆ A , where ˆ A = A − BB T X ∞ . The identity ( I n − P ( t )) − = I n + ( I n − P ( t )) − P ( t ) leads to X ( t ) = X ∞ − e t ˆ A T X ∞ e t ˆ A − e t ˆ A T X ∞ (cid:16) I n − (cid:16) X L − e t ˆ A X L e t ˆ A T (cid:17) X ∞ (cid:17) − (cid:16) X L − e t ˆ A X L e t ˆ A T (cid:17) X ∞ e t ˆ A . (14) Derivation by using the exact solution X ∞ of the ARE By Lemma 5.1 it holds that range ( X ∞ ) is invariant under ˆ A T . Assume now that X ∞ is given infactorized form, this means that X ∞ = Z ∞ Z T ∞ and Z ∞ ∈ R n × p and 1 ≤ p = rank ( X ∞ ) ≤ n . Ifrank( X ∞ ) = 0, then also X ∞ = 0 as well as the solution X ( t ). Now it holds that range( X ∞ ) =range( Z ∞ ) and consequently range( Z ∞ ) is invariant under ˆ A T .By means of the compact singular value decomposition of Z ∞ , we obtain matrices Q ∞ ∈ R n × p , S ∞ ∈ R p × p and V ∞ ∈ R p × p , such that Z ∞ = Q ∞ S ∞ V T ∞ , range ( Q ∞ ) = range ( Z ∞ ) and Z ∞ Z T ∞ = Q ∞ S ∞ Q T ∞ . Because of the invariance we get e t ˆ A T Q ∞ = Q ∞ e tQ T ∞ ˆ A T Q ∞ . Now observe that e t ˆ A T X ∞ = e t ˆ A T Z ∞ Z T ∞ = e t ˆ A T Q ∞ S ∞ Q T ∞ = Q ∞ e tQ T ∞ ˆ A T Q ∞ S ∞ Q T ∞ . (15)15herefore the solution X ( t ) can be written in the form X ( t ) = X ∞ − Q ∞ ˜ X ( t ) Q T ∞ . (16)We use the DRE (1) and equation (16) and get a differential equation for ˜ X ( t )˙˜ X ( t ) = Q T ∞ ˆ A T Q ∞ ˜ X ( t ) + ˜ X ( t ) Q T ∞ ˆ AQ ∞ + ˜ XQ T ∞ BB T Q ∞ ˜ X ( t ) , (17a)˜ X (0) = Q T ∞ X ∞ Q ∞ . (17b) Derivation by using a low-rank approximation X N of the exact solution X ∞ of theARE Let now Z N Z TN = X N ≈ X ∞ be a low-rank approximation obtained by a numerical method. Wereplace X ∞ by X N in formula 14 and obtain X ( t ) ≈ X N − e t ( A − BB T X N ) T X N e t ( A − BB T X N ) − e t ( A − BB T X N ) T X N (cid:16) I n − (cid:16) X L − e t ˆ A X L e t ˆ A T (cid:17) X ∞ (cid:17) − (cid:16) X L − e t ˆ A X L e t ˆ A T (cid:17) X N e t ( A − BB T X N ) . Let Z N = Q N S N V TN be the compact singular value decomposition of the low-rank factor. Accordingto formula 15, we propose to approximate the action of the matrix exponential by e t ( A − BB T X N ) T X N = e t ( A − BB T X N ) T Z N Z TN = e t ( A − BB T X N ) T Q N S N Q TN ≈ Q N e tQ TN ( A − BB T X N ) T Q N S N Q TN . Therefore we obtain the Galerkin ansatz X ( t ) ≈ X N − Q N ˜ X N ( t ) Q TN for the numerical approximation.Again we use the DRE (1) and get a differential equation for ˜ X N ( t )˙˜ X N ( t ) = Q TN (cid:16) A − BB T X N (cid:17) T Q N ˜ X N ( t ) + ˜ X N ( t ) Q TN (cid:16) A − BB T X N (cid:17) Q N + ˜ X N ( t ) Q TN BB T Q N ˜ X N ( t ) + Q TN R ( X N ) Q N . ˜ X N (0) = Q TN X N Q N . We assume that the numerical low-rank approximation is accurate enough such that R ( X N ) ≈ (cid:13)(cid:13)(cid:13) Q TN R ( X N ) Q N (cid:13)(cid:13)(cid:13) ≤ kR ( X N ) k ≈ . This means that the projected residual Q TN R ( X N ) Q N is even smaller than the residual of the ARE R ( X N ) and, therefore, we can neglect the residual. Let X ∞ = Z ∞ Z T ∞ be the exact solution of the ARE (2). Moreover let Z ∞ = Q ∞ S ∞ V T ∞ be itscompact singular value decomposition, such that Q ∞ ∈ R n × p , S ∞ ∈ R p × p and V ∞ ∈ R p × p and Z ∞ = Q ∞ S ∞ V T ∞ . The compact singular value decomposition of Z ∞ gives a spectral decompositionof X ∞ that is X ∞ = Z ∞ Z T ∞ = Q ∞ S ∞ Q T ∞ , and S ∞ = diag (cid:16) λ ↓ ( X ∞ ) , . . . , λ ↓ p ( X ∞ ) (cid:17) . S ∞ contains all non-zero eigenvalues of X ∞ in a non-increasingfashion. We have that range ( X ∞ ) = range ( Z ∞ ) = range ( Q ∞ ). Because of Theorem 3.3 it holdsthat range ( Q ∞ ) = K (cid:16) A T , C T (cid:17) . According to Theorem 3.5 we can represent the solution in thefollowing form X ( t ) = Q ∞ Q T ∞ X ( t ) Q ∞ Q T ∞ . This representation has the advantage that the entries of Q T ∞ X ( t ) Q ∞ can be bounded by theeigenvalues of X ∞ . Theorem 5.1.
Let ( A, B ) be stabilizable and ( A, C ) be detectable. Moreover let X ∞ ∈ R n × n be the unique symmetricpositive semidefinite solution of the ARE (2) and q , . . . , q n ∈ R n be a system of orthonormaleigenvectors of X ∞ corresponding to the eigenvalues λ ↓ ( X ∞ ) , . . . , λ ↓ n ( X ∞ ) ∈ R . Then for all i, j = 1 , . . . , n and t ≥ the following holds: (cid:12)(cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12)(cid:12) ≤ q λ ↓ i ( X ∞ ) λ ↓ j ( X ∞ ) , (18) where X is the unique solution of the DRE (1) with X = 0 .Proof. According to Theorem 3.1 the inequality 0 (cid:52) X ( t ) (cid:52) X ∞ holds for all t ≥
0. By multiplyingthe inequality with q Ti from the left and q i from the right we obtain0 ≤ q Ti X ( t ) q i ≤ q Ti X ∞ q i = λ ↓ i ( X ∞ ) q Ti q i = λ ↓ i ( X ∞ ) . Now let i = j and α, β ∈ R . Again by multiplying the inequality with αq i + βq j we obtain0 ≤ ( αq i + βq j ) T X ( t ) ( αq i + βq j ) ≤ ( αq i + βq j ) T X ∞ ( αq i + βq j ) . Since X ( t ) and X ∞ are symmetric, it applies that α q Ti X ( t ) q i + 2 αβq Ti X ( t ) q j + β q Tj X ( t ) q j ≤ α q Ti X ∞ q i + 2 αβq Ti X ∞ q j + β q Tj X ∞ q j . As q i and q j are different orthonormal eigenvectors of X ∞ , we obtain for the right hand side α q Ti X ∞ q i + 2 αβq Ti X ∞ q j + β q Tj X ∞ q j = α λ ↓ i ( X ∞ ) + 2 αβλ ↓ j ( X ∞ ) q Ti q j + β λ ↓ i ( X ∞ )= α λ ↓ i ( X ∞ ) + β λ ↓ i ( X ∞ ) . As X ( t ) is symmetric positive semidefinite, the following inequality holds for the left hand side. α q Ti X ( t ) q i + 2 αβq Ti X ( t ) q j + β q Tj X ( t ) q j ≥ αβq Ti X ( t ) q j . Now we have0 ≤ α λ ↓ i ( X ∞ ) − αβq Ti X ( t ) q j + β λ ↓ j ( X ∞ ) = h α β i " λ ↓ i ( X ∞ ) − q iT X ( t ) q j − q iT X ( t ) q j λ ↓ j ( X ∞ ) αβ . Since this holds for all α, β ∈ R the matrix " λ ↓ i ( X ∞ ) − q iT X ( t ) q j − q iT X ( t ) q j λ ↓ j ( X ∞ )
17s symmetric positive semidefinite. Therefore its determinant must be non-negative,0 ≤ λ ↓ i ( X ∞ ) λ ↓ j ( X ∞ ) − (cid:16) q iT X ( t ) q j (cid:17) . Finally this leads to (cid:12)(cid:12)(cid:12) q iT X ( t ) q j (cid:12)(cid:12)(cid:12) ≤ q λ ↓ i ( X ∞ ) λ ↓ j ( X ∞ ) . Let the columns of Q ∞ be q , . . . , q p . Due to the decay of the eigenvalues λ ↓ k ( X ∞ ) of the solutionof the ARE (2) and the inequality (18) from Theorem 5.1, the values (cid:12)(cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12)(cid:12) also decay for i + j increasing. We have that X ( t ) = Q ∞ Q T ∞ X ( t ) Q ∞ Q T ∞ = p X i,j =1 (cid:16) q Ti X ( t ) q j (cid:17) q i q Tj . For quick enough eigenvalue decay, we expect that (cid:12)(cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12)(cid:12) ≤ q λ ↓ i ( X ∞ ) λ ↓ j ( X ∞ ) ≈ i + j large enough. We truncate the series and obtain X ( t ) ≈ k X i,j =1 (cid:16) q Ti X ( t ) q j (cid:17) q i q Tj = Q ∞ ,k Q T ∞ ,k X ( t ) Q ∞ ,k Q T ∞ ,k , where Q ∞ ,k = h q , . . . , q k i ∈ R n × k . We also consider the appropriate real linear space Q ∞ ,k := n Q ∞ ,k Y Q T ∞ ,k | Y ∈ R q × q o ⊆ R n × n together with the orthogonal projection P k : R n × n → Q ∞ ,k , P ∞ ,k ( X ) = Q ∞ ,k Q T ∞ ,k XQ ∞ ,k Q T ∞ ,k As the columns of Q ∞ ,k are orthonormal, it holds that P ∞ ,k ( X ) = P ∞ ,k ( X ). Moreover the projection P ∞ ,k is orthogonal, because h X − P ∞ ,k ( X ) , Q ∞ ,k Y Q ∞ ,k i F = h X − Q ∞ ,k Q T ∞ ,k XQ ∞ ,k Q T ∞ ,k , Q ∞ ,k Y Q T ∞ ,k i F = h X, Q ∞ ,k Y Q T ∞ ,k i F − h Q ∞ ,k Q T ∞ ,k XQ ∞ ,k Q T ∞ ,k , Q ∞ ,k Y Q T ∞ ,k i F = h X, Q ∞ ,k Y Q T ∞ ,k i F − h X, Q ∞ ,k Y Q T ∞ ,k i F = 0for all Y ∈ R k × k . Therefore the best approximation of X ( t ) in Q ∞ ,k is given by k X i,j =1 (cid:16) q Ti X ( t ) q j (cid:17) q i q Tj = P ∞ ,k ( X ( t )) = argmin X ∈Q ∞ ,k k X − X ( t ) k F and for the approximation error we obtain k X ( t ) − P ∞ ,k ( X ( t )) k F = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) p X i,j =1 i>k ∨ j>k (cid:16) q Ti X ( t ) q j (cid:17) q i q Tj (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) F = vuuuut p X i,j =1 i>k ∨ j>k (cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12) vuuuut p X i,j =1 i>k ∨ j>k λ ↓ i ( X ∞ ) λ ↓ j ( X ∞ ) . Since the eigenvalues λ ↓ p +1 ( X ∞ ) , . . . , λ ↓ n ( X ∞ ) are 0 we obtain k X ( t ) − P ∞ ,k ( X ( t )) k F ≤ vuuuut n X i,j =1 i>k ∨ j>k λ ↓ i ( X ∞ ) λ ↓ j ( X ∞ ) . (19)We propose therefore to setup a trial space for the Galerkin approach using a system of eigenvectorscorresponding to the largest eigenvalues. This can be obtained by using a low-rank method to obtaina numerical approximation of the solution of the ARE. Then a compact singular value decompositionof the numerical low-rank approximation of X ∞ can be used to obtain an approximation of theeigenvectors corresponding to the largest eigenvalues. The small singular values can be safelytruncated from the singular value decomposition by virtue of Thm. 5.1. This reduces also thedimension of the trial space. Let Z N = Q N S N V TN . be the truncated reduced singular value decomposition of the low-rank approximation. With that,the trial space for the Galerkin approach is given by n Q N ˜ XQ TN | ˜ X ∈ R p × p o , and, as X ( t ) converges to X ∞ and X ∞ ≈ Z N Z TN , we propose the Galerkin ansatz X ( t ) ≈ Z N Z TN − Q N ˜ X ( t ) Q TN . Example 5.1 (Decay of Absolute Values of Entries) . We illustrate the decay of (cid:12)(cid:12)(cid:12) q i X ( t ) q Tj (cid:12)(cid:12)(cid:12) in Figures 4-8 . We have chosen the same matrices as for the
Example 3.1 . To improve the visualization all values below machine precision were set to machineprecision. The eigenvalue decay of the solution X ∞ of the corresponding ARE is shown in
Figure 9 . − − − − i j (cid:12)(cid:12) q T i X ( t ) q j (cid:12)(cid:12) − − − − − − − − Fig. 4. Decay of (cid:12)(cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12)(cid:12) for t = 1. − − − − i j (cid:12)(cid:12) q T i X ( t ) q j (cid:12)(cid:12) − − − − − − − − Fig. 5. Decay of (cid:12)(cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12)(cid:12) for t = 3. 19
25 50 75 125507510 − − − − i j (cid:12)(cid:12) q T i X ( t ) q j (cid:12)(cid:12) − − − − − − − − Fig. 6. Decay of (cid:12)(cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12)(cid:12) for t = 5. − − − − i j (cid:12)(cid:12) q T i X ( t ) q j (cid:12)(cid:12) − − − − − − − − Fig. 7. Decay of (cid:12)(cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12)(cid:12) for t = 7. − − − − i j (cid:12)(cid:12) q T i X ( t ) q j (cid:12)(cid:12) − − − − − − − − Fig. 8. Decay of (cid:12)(cid:12)(cid:12) q Ti X ( t ) q j (cid:12)(cid:12)(cid:12) for t = 9. − − − − − k λ ↓ k ( X ∞ ) Fig. 9. The eigenvalue decay of X ∞ . Remark 5.1.
With minor adjustments, all arguments also hold for the generalized
DRE M T ˙ X ( t ) M = A T X ( t ) M + M T X ( t ) A − M T X ( t ) BB T X ( t ) M + C T C, (20a) X (0) = 0 , (20b) with M ∈ R n × n nonsingular that can accommodate, e.g., a mass matrix from a finite elementdiscretization. In summary, the proposed approach reads as written down in Algorithm 3.20 lgorithm 3
Galerkin approach for the generalized DRE (20) (
ARE-Galerkin ) Assumption: (cid:0) AM − , B (cid:1) is stabilizable and ( AM − , CM − ) is detectable. Input:
M, A ∈ R n × n , B ∈ R n × b , C ∈ R c × n . Output: X ( t ) ≈ Z ∞ Z T ∞ − Q ∞ ˜ X ( t ) Q T ∞ that approximates the solution to M T ˙ X ( t ) M = A T X ( t ) M + M T X ( t ) A − M T X ( t ) BB T X ( t ) M + C T C, X (0) = 0.% Solve the ARE: A T X ∞ M + M T X ∞ A − M T X ∞ BB T X ∞ M + C T C = 0 for X ∞ ≈ Z ∞ Z T ∞ and Z ∞ ∈ R n × q ;% Compute compact singular value decomposition: [ Q ∞ , S ∞ , ∼ ] = svd ( Z ∞ , tol = ε machine · S ∞ (1 , idx = diag ( S ∞ ) ≥ tol ; S ∞ = S ∞ ( idx, idx ); Q ∞ = Q ∞ (: , idx ); Z ∞ = Q ∞ S ∞ ;% Compute matrices: A F = Q T ∞ (cid:16) AM − − BB T Z ∞ Z T ∞ (cid:17) Q ∞ ; B F = Q T ∞ B ;% Solve the differential equation using Algorithm 2: ˙˜ X ( t ) = A TF ˜ X ( t ) + ˜ X ( t ) A F + ˜ X ( t ) B F B TF ˜ X ( t ) , ˜ X (0) = S ∞ ; 21 . Numerical Experiments To quantify the performance of Algorithm 3, we consider a number of differential Riccati equationsthat are used to define optimal controls. Concretely, we consider the generalized differential Riccatiequation M T ˙ X ( t ) M = A T X ( t ) M + M T X ( t ) A − M T X ( t ) BB T M + C T C, (21a) X (0) = 0 . (21b)and their realizations. First, we consider the RAIL benchmark example, that is a finite elementdiscretization of a heat equation; see [16] for the model description. The second example,
CONV DIFF ,derives from a finite-differences discretized heat equation with convection on the unit square withhomogenous Dirichlet boundary conditions, ∂∂t x ( ξ, t ) − ∆ x ( ξ, t ) − v · ∇ x ( ξ, t ) = f ( ξ ) u ( t ) in Ω × (0 , T )where Ω = (0 , and v = [10 , T ; see [46].On both examples, we compare the proposed method with the splitting methods developed in [52,53].The splitting methods are based on a splitting of the DRE into an affine and nonlinear subproblem.The advantages of that approach lie in the fact that the nonlinear subproblem can be solved by anexplicit solution formula. The numerical solution of the linear subproblem is based on approximatingthe action of a matrix exponential by means of Krylov subspace methods. We used the MATLABimplementation DREsplit [55] of the splitting methods for our experiments. In the tests, we employedthe
Lie and
Strang splitting of order 1 and 2 respectively, as well as the symmetric splitting of order4 , LIE , STRANG , SYMMETRIC2 , SYMMETRIC4 , SYMMETRIC6 and
SYMMETRIC8 .To evaluate the error, we computed a reference solution X ref ( t ) using SYMMETRIC8 with constanttime step size h . The basic information about the setup of the benchmark problems are given inTable 1.problem n matrices interval reference solution M symmetric positive definite, RAIL A symmetric, M − A stable, [0 , SYMMETRIC8 , h = 2 − B ∈ R n × , C ∈ R × n M = I n , CONV DIFF A nonsymmetric and stable, [0 , . SYMMETRIC8 , h = 2 − B ∈ R n × , C ∈ R × n Table 1: Information about benchmark problems.All computations are carried out on a machine with 2 × Xeon ® Skylake Silver 4110 @ 2.10GHzCPU with 8 cores, 192 GB Ram and MATLAB 2018a. We have used the low-rank Newton ADIiteration implemented in
MEX-M.E.S.S. [12] to solve the algebraic Riccati equations; as required forour approach as laid out in Algorithm 3.We report the absolute and relative errors k X ( t ) − X ref ( t ) k and k X ( t ) − X ref ( t ) kk X ref ( t ) k , X ( t ) is the numerical approximation and X ref ( t ) is the reference solution in 2-norm andFrobenius norm. We also report the norm of the reference solution k X ref ( t ) k as well as theconvergence to the stationary point k X ref ( t ) − X ∞ k .Numerical results for the Galerkin approximation from Algorithm 3 and for the splitting schemebased solvers and be found in Appendices A and B. The computational costs for both methods aregiven in Section 6.2. Also, we evaluate the best approximation in the trial space of the referencesolution, which is given by X best ( t ) := Q ∞ Q T ∞ X ref ( t ) Q ∞ Q T ∞ = argmin X ∈ { Q ∞ ˜ XQ T ∞ | ˜ X ∈ R p × p } k X − X ref ( t ) k F , where Q ∞ is the matrix from Algorithm 3 Line 8.The code of the implementation and the precomputed reference solution are available as mentionedin Figure 10. Code and Data Availability
The source code of the implementations used to compute the presented results is available from: doi:10.5281/zenodo.2629737https://gitlab.mpi-magdeburg.mpg.de/behr/behbh19 dre are galerkin code under the GPLv2+ license and is authored by Maximilian Behr.Fig. 10. Link to code and data.
The initial step of Algorithm 3 requires the solution to the associated ARE. For this task we call
MEX-M.E.S.S. that iteratively computes the numerical solution to the following absolute and relativeresiduals (cid:13)(cid:13)(cid:13) A T Z ∞ Z T ∞ M + M T Z ∞ Z T ∞ A − M T Z ∞ Z T ∞ BB T Z ∞ Z T ∞ M + C T C (cid:13)(cid:13)(cid:13) and (cid:13)(cid:13)(cid:13) A T Z ∞ Z T ∞ M + M T Z ∞ Z T ∞ A − M T Z ∞ Z T ∞ BB T Z ∞ Z T ∞ M + C T C (cid:13)(cid:13)(cid:13) k C T C k . The achieved values for the different test setups as well as the number of columns of the corresponding Z ∞ after truncation (see Step 5 of Algorithm 3), that define the dimension of the reduced model,are listed in Table 2. instance n size of Galerkin system absolute residual relative residual RAIL ,
177 319 5 . · − . · − CONV DIFF ,
400 56 1 . · − . · − Table 2: Residuals for the ARE 0 = A T XM + M XA − M T XBB T XM + C T C .The 1-norm bound for the matrix exponential tol exp from Algorithm 2 was set to 1 · . The23esulting step sizes are given in Table 3. Instance n Step sizes h RAIL , (cid:8) , − , . . . , − (cid:9) CONV DIFF , (cid:8) − , − , . . . , − (cid:9) Table 3: Step sizes h for modified Davison-Maki method Algorithm 2.We plot the numerical errors in Figures 15–18 and 21–24. Figures 19, 20, 25 and 26 show the normof the reference solution and the convergence to the stationary point.In view of the performance, we can interpret the presented numbers and plots as follows: Firstly,the accuracy of the modified Davison-Maki method ; cf. Figure 16 and 18 is independent of the stepsize, as discussed in Section 4.3. Still we compute the solution on different time grids, since forcontrol applications the values of the solution might be needed at many time instances.The computational times for ARE-Galerkin include the solve of the corresponding ARE and thesubsequent integration of the projected dense DRE. Since the efforts for the time integration exactlydoubles with a bisection of the step size, from the timings for the
RAIL problem, with, e.g., 42s( h = 2 − ) and 77s ( h = 2 − ) (see Figure 11), one infers that most of the time is spent to solvethe dense DRE. Conversely, for the CONV DIFF benchmark problem, most of the time (45s) wasused to solve the ARE. As the resulting Galerkin projected DRE system is of size 56 only, thecomputational costs for the time integration are vanishingly small. Accordingly, the differences inthe effort caused by finer time grids are hardly visible; see Figure 13.The reference solution for the
RAIL problem is large in norm what makes the absolute errorcomparatively large; see the Figure 19 in Appendix A.In both examples, in terms of accuracy, the
ARE-Galerkin approximation is nearly at the samelevel as the high order splitting schemes, cf. Figures 16, 32 and Figures 22, 38. We note, however,that the
ARE-Galerkin method does not give the best possible approximation in the trial space;compare the error levels for X best .In any case, the ARE-Galerkin method clearly outperforms the splitting methods in terms ofcomputational time versus accuracy in all test examples. − − − − − · · · · · RAIL , n = 5177, Step Size h C o m pu t a t i o n a l T i m e i nS ec o nd s Fig. 11. Timing for
ARE-Galerkin . − − · · · · · RAIL , n = 5177, Step Size h C o m pu t a t i o n a l T i m e i nS ec o nd s Fig. 12. Timing for splitting schemes.24 − − − − − · · · · · CONV DIFF , n = 6400, Step Size h C o m pu t a t i o n a l T i m e i nS ec o nd s Fig. 13. Timing for
ARE-Galerkin . − − − − · · · · · CONV DIFF , n = 6400, Step Size h C o m pu t a t i o n a l T i m e i nS ec o nd s Fig. 14. Timing for splitting schemes.
LIE STRANG SYMMETRIC2 SYMMETRIC4 SYMMETRIC6 SYMMETRIC8ARE-Galerkin . Conclusion We have reviewed and extended fundamental properties of the solution to the differential and algebraicRiccati equation and heavily relied on the solution representation provided by Radon’s Lemma toanalyze variants of
Davison-Maki methods and to derive an efficient Galerkin projection scheme.Numerical tests confirmed that the resulting projected scheme outperforms existing methods interms of computation time, memory requirements, and approximation quality. In particular, storagerequirements have been the bottleneck in the numerical considerations of large-scale differentialRiccati equations.Our proposed Galerkin ansatz bases on a low-rank approximation of the associated algebraicRiccati equation (ARE) for which there are efficient solvers. Moreover, the information on theresidual and on eigenvalue decay, that come with the low-rank iteration for the ARE can be directlytransferred into estimates for the approximation quality of our approach the more that the use ofthe
Davison-Maki methods leads to an exact time discretization.Future work will deal with the treatment of nonzero initial conditions. While the formulas areeasily extended to this case, the invariance properties and the eigenvalue comparisons, that were thebackbone of our numerical approach, are no longer given in general. For the (not so) special casethat the initial condition X writes as X = C T W C with a symmetric positive definite weightingmatrix W , the flow invariance as established in Section 5.1 still holds so that the presented algorithmcan be applied without modification. For a general low-rank initial condition X = Z Z T the flowinvariance can be achieved by taking the columns of the solution to the ARE (2) with C T C replacedby h C T , Z i h C T , Z i T as the Galerkin ansatz space. In any case, the inequality 0 (cid:52) X ( t ) (cid:52) X ∞ ,where X ∞ is the solution of the ARE, does not hold anymore and, thus, approximation qualitycannot be assured as in (19). However, if one can find a matrix ˜ X for which the comparison0 (cid:52) X ( t ) (cid:52) ˜ X holds, all arguments of Section 5.2 apply accordingly. References [1] H. Abou-Kandil, G. Freiling, V. Ionescu, and G. Jank.
Matrix Riccati Equations in Controland Systems Theory . Birkh¨auser, Basel, Switzerland, 2003.[2] L. Amodei and J.-M. Buchot. An invariant subspace method for large-scale algebraic Riccatiequation.
Appl. Numer. Math. , 60(11):1067–1082, 2010.[3] B. D. O. Anderson and J. B. Moore.
Linear Optimal Control . Prentice-Hall, Englewood Cliffs,NJ, 1971.[4] V. Angelova, M. Hached, and K. Jbilou. Approximate solutions to large nonsymmet-ric differential Riccati problems with applications to transport theory. Technical ReportarXiv:1801.01291v2, arXiv, 2019. math.NA.[5] A. C. Antoulas, D. C. Sorensen, and Y. Zhou. On the decay rate of Hankel singular values andrelated issues.
Syst. Cont. Lett. , 46(5):323–342, 2002.[6] J. Baker, M. Embree, and J. Sabino. Fast singular value decay for Lyapunov solutions withnonnormal coefficients.
SIAM J. Matrix Anal. Appl. , 36(2):656–668, 2015.[7] M. Beck and S. J. A. Malham. Computing the Maslov index for large systems.
Proc. Amer.Math. Soc. , 143(5):2159–2173, 2015.[8] B. Beckermann and A. Townsend. On the singular values of matrices with displacementstructure.
SIAM J. Matrix Anal. Appl. , 38(4):1227–1248, 2017.269] M. Behr, P. Benner, and J. Heiland. On an Invariance Principle for the Solution Space of theDifferential Riccati Equation.
Proc. Appl. Math. Mech. , 18(1), 2018.[10] M. Behr, P. Benner, and J. Heiland. Solution formulas for differential Sylvester and Lyapunovequations. e-print 1811.08327, arXiv, November 2018. math.NA.[11] P. Benner and Z. Bujanovi´c. On the solution of large-scale algebraic Riccati equations by usinglow-dimensional invariant subspaces.
Linear Algebra Appl. , 488:430–459, 2016.[12] P. Benner, M. K¨ohler, and J. Saak. M.E.S.S. – Matrix Equations Sparse Solver. .[13] P. Benner and H. Mena. BDF methods for large-scale differential Riccati equations. InB. De Moor, B. Motmans, J. Willems, P. Van Dooren, and V. Blondel, editors,
Proc. 16th Intl.Symp. Mathematical Theory of Network and Systems, MTNS 2004 , 2004.[14] P. Benner and H. Mena. Rosenbrock methods for solving Riccati differential equations.
IEEETrans. Autom. Control , 58(11):2950–2957, 2013.[15] P. Benner and H. Mena. Numerical solution of the infinite-dimensional LQR-problem and theassociated differential Riccati equations.
J. Numer. Math. , 26(1):1–20, 2018.[16] P. Benner and J. Saak. A semi-discretized heat transfer model for optimal cooling of steelprofiles. In P. Benner, V. Mehrmann, and D. Sorensen, editors,
Dimension Reduction of Large-Scale Systems , volume 45 of
Lect. Notes Comput. Sci. Eng. , pages 353–356. Springer-Verlag,Berlin/Heidelberg, Germany, 2005.[17] R. A. Brocket.
Finite Dimensional Linear Systems . Wiley, New York, 1970.[18] F. M. Callier, J. Winkin, and J. L. Willems. Convergence of the time-invariant Riccati differentialequation and LQ-problem: mechanisms of attraction.
Internat. J. Control , 59(4):983–1000,1994.[19] C. H. Choi. A survey of numerical methods for solving matrix Riccati differential equations. In
IEEE Proceedings on Southeastcon , pages 696–700 vol.2, 1990.[20] E. J. Davison and M. C. Maki. The numerical solution of the matrix Riccati differentialequation.
IEEE Trans. Autom. Control , 18:71–73, 1973.[21] L. Grasedyck. Existence of a low rank or H -matrix approximant to the solution of a Sylvesterequation. Numer. Lin. Alg. Appl. , 11(4):371–389, 2004.[22] L. Grubiˇsi´c and D. Kressner. On the eigenvalue decay of solutions to operator Lyapunovequations.
Syst. Cont. Lett. , 73:42–47, 2014.[23] Y. G¨uldoˇgan, M. Hached, K. Jbilou, and M. Kurulay. Low-rank approximate solutions tolarge-scale differential matrix Riccati equations.
Applicationes Mathematicae , 45(2):233–254,2018.[24] M. Hached and K. Jbilou. Computational Krylov-based methods for large-scale differentialSylvester matrix problems.
Numer. Lin. Alg. Appl. , 25(5):e2187, 14, 2018.[25] M. Hached and K. Jbilou. Numerical methods for differential linear matrix equations via Krylovsubspace methods. Technical Report arXiv:1805.10192v1, arXiv, 2018. math.NA. 2726] M. Hached and K. Jbilou. Numerical solutions to large-scale differential Lyapunov matrixequations.
Numer. Algorithms , 2018.[27] J. Heiland.
Decoupling and Optimization of Differential-Algebraic Equations with Applicationin Flow Control . Dissertation, TU Berlin, 2014.[28] R. A. Horn and C. R. Johnson.
Matrix Analysis . Cambridge University Press, Cambridge,1985.[29] R. E. Kalman and T. S. Englar. A user’s manual for the automatic synthesis program. RIASReport CR-475, NASA, 1966.[30] C. Kenney and R. B. Leipnik. Numerical integration of the differential matrix Riccati equation.
IEEE Trans. Autom. Control , 30:962–970, 1985.[31] G. Kirsten and V. Simoncini. Order reduction methods for solving large-scale differential matrixRiccati equations. e-print 1905.12119, arXiv, 2019. math.NA.[32] H. W. Knobloch and H. Kwakernaak.
Lineare Kontrolltheorie . Springer-Verlag, Berlin, 1985.In German.[33] M. K¨ohler, N. Lang, and J. Saak. Solving differential matrix equations using Parareal.
Proc.Appl. Math. Mech. , 16(1):847–848, 2016.[34] A. Koskela and H. Mena. A structure preserving Krylov subspace method for large scaledifferential Riccati equations. e-print arXiv:1705.07507, arXiv, 2017. math.NA.[35] N. Lang.
Numerical Methods for Large-Scale Linear Time-Varying Control Systems and relatedDifferential Matrix Equations . Dissertation, Technische Universit¨at Chemnitz, Germany, 2017.[36] N. Lang, H. Mena, and J. Saak. On the benefits of the
LDL T factorization for large-scaledifferential matrix equation solvers. Linear Algebra Appl. , 480:44–71, 2015.[37] A. J. Laub. Schur techniques for Riccati differential equations. In D. Hinrichsen and A. Isidori,editors,
Feedback Control of Linear and Nonlinear Systems , pages 165–174. Springer-Verlag,New York, 1982.[38] A. Locatelli.
Optimal Control: An Introduction . Birkh¨auser, Basel, Switzerland, 2001.[39] C. Martin. Grassmannian manifolds, Riccati Equations and Feedback Invariants of LinearSystems. In
Geometrical Methods for the Theory of Linear Systems , pages 195–211. Springer,1980.[40] T. McCauley. Computing the Maslov index from singularities of a matrix Riccati equation.
J.Dyn. Diff. Equat. , 29(4):1487–1502, 2017.[41] H. Mena.
Numerical Solution of Differential Riccati Equations Arising in Optimal Control ofPartial Differential Equations . Dissertation, Escuela Polit´ecnica Nacional, Ecuador, 2007.[42] H. Mena, L.-M. Pfurtscheller, and T. Stillfjord. GPU acceleration of splitting schemes appliedto differential matrix equations.
Numer. Algorithms , 2019.[43] K. R. Meyer and D. C. Offin.
Introduction to Hamiltonian Dynamical Systems and the N-BodyProblem , volume 90 of
Applied Mathematical Sciences . Springer, Cham, third edition, 2017.2844] M. Opmeer. Decay of singular values of the Gramians of infinite-dimensional systems. In
Proceedings 2015 European Control Conference (ECC) , pages 1183–1188, Linz, Austria, 2015.IEEE.[45] T. Penzl. Eigenvalue decay bounds for solutions of Lyapunov equations: the symmetric case.
Syst. Cont. Lett. , 40:139–144, 2000.[46] T. Penzl.
Lyapack
Users Guide. Technical Report SFB393/00-33, Sonderforschungsbereich393
Numerische Simulation auf massiv parallelen Rechnern , TU Chemnitz, 09107 Chemnitz,Germany, 2000. Available from .[47] V. Radisavljevic. Improved Potter-Anderson-Moore algorithm for the differential Riccatiequation.
Appl. Math. Comput. , 218(8):4641–4646, 2011.[48] I. Rusnak. Almost analytic representation for the solution of the differential matrix Riccatiequation.
IEEE Trans. Autom. Control , 33(2):191–193, 1988.[49] C. R. Schneider. Global aspects of the matrix Riccati equation.
Math. Systems Theory ,7(3):281–286, 1973.[50] M. A. Shayman. Phase portrait of the matrix Riccati equation.
SIAM J. Control Optim. ,24(1):1–65, 1986.[51] D. C. Sorensen and Y. Zhou. Bounds on eigenvalue decay rates and sensitivity of solutions toLyapunov equations. Technical Report TR02-07, Dept. of Comp. Appl. Math., Rice University,Houston, TX, 2002.[52] T. Stillfjord. Low-rank second-order splitting of large-scale differential Riccati equations.
IEEETrans. Autom. Control , 60(10):2791–2796, 2015.[53] T. Stillfjord. Adaptive high-order splitting schemes for large-scale differential Riccati equations.
Numer. Algorithms , 78:1129–1151, 2018.[54] T. Stillfjord. Singular value decay of operator-valued differential Lyapunov and Riccati equations.
SIAM J. Control Optim. , 56:3598–3618, 2018.[55] T. Stilljford. DREsplit. .[56] D. R. Vaughan. A negative exponential solution for the matrix Riccati equation.
IEEE Trans.Autom. Control , 14:72–75, 1969.[57] W. Walter.
Ordinary differential equations , volume 182 of
Graduate Texts in Mathematics .Springer-Verlag, New York, 1998. 29 . Numerical Results for Galerkin Approach
RAIL , n = 5177 and M T ˙ X ( t ) M = A T X ( t ) M + M T X ( t ) A − M T X ( t ) BB T X ( t ) M + C T C, X (0) = 0. − − − − − RAIL , n = 5177, t ∈ [0 , h = 2 a b s o l u t ee rr o r Fig. 15. Absolute error of the Galerkin andBest approximation. − − − − − − − − − − − RAIL , n = 5177, t ∈ [0 , h = 2 r e l a t i v ee rr o r Fig. 16. Relative error of the Galerkin andBest approximation. − − − − − RAIL , n = 5177, t ∈ [0 , h = 2 − a b s o l u t ee rr o r Fig. 17. Absolute error of the Galerkin andBest approximation. − − − − − − − − − − − RAIL , n = 5177, t ∈ [0 , h = 2 − r e l a t i v ee rr o r Fig. 18. Relative error of the Galerkin andBest approximation.
ARE-Galerkin X best ( t ) RAIL , n = 5177, t ∈ [0 , n o r m o f r e f e r e n ce s o l u t i o n k X ref ( t ) k k X ref ( t ) k F Fig. 19. Norm of the reference solution. RAIL , n = 5177, t ∈ [0 , c o n v e r g e n ce t o s t a t i o n a r y p o i n t k X ref ( t ) − X ∞ k Fig. 20. Convergence to the stationary point.30
ONV DIFF , n = 6400 and ˙ X ( t ) = A T X ( t ) + X ( t ) A − X ( t ) BB T X ( t ) + C T C, X (0) = 0. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − a b s o l u t ee rr o r Fig. 21. Absolute error of the Galerkin andBest approximation. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − r e l a t i v ee rr o r Fig. 22. Relative error of the Galerkin andBest approximation. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − a b s o l u t ee rr o r Fig. 23. Absolute error of the Galerkin andBest approximation. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − r e l a t i v ee rr o r Fig. 24. Relative error of the Galerkin andBest approximation.
ARE-Galerkin X best ( t ) .
02 0 .
04 0 .
06 0 .
08 0 . . − CONV DIFF , n = 6400, t ∈ [0 , . n o r m o f r e f e r e n ce s o l u t i o n k X ref ( t ) k k X ref ( t ) k F Fig. 25. Norm of the reference solution. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . c o n v e r g e n ce t o s t a t i o n a r y p o i n t k X ref ( t ) − X ∞ k Fig. 26. Convergence to the stationary point.31 . Numerical Results for Splitting Schemes
RAIL , n = 5177 and M T ˙ X ( t ) M = A T X ( t ) M + M T X ( t ) A − M T X ( t ) BB T X ( t ) M + C T C, X (0) = 0. − − − RAIL , n = 5177, t ∈ [0 , h = 2 a b s o l u t ee rr o r Fig. 27. Absolute error of the splittingscheme approximation. − − − − − − − − − − − − − RAIL , n = 5177, t ∈ [0 , h = 2 r e l a t i v ee rr o r Fig. 28. Relative error of the splitting schemeapproximation. − − − RAIL , n = 5177, t ∈ [0 , h = 2 a b s o l u t ee rr o r Fig. 29. Absolute error of the splittingscheme approximation. − − − − − − − − − − − − − RAIL , n = 5177, t ∈ [0 , h = 2 r e l a t i v ee rr o r Fig. 30. Relative error of the splitting schemeapproximation. − − − RAIL , n = 5177, t ∈ [0 , h = 2 − a b s o l u t ee rr o r Fig. 31. Absolute error of the splittingscheme approximation. − − − − − − − − − − − − − RAIL , n = 5177, t ∈ [0 , h = 2 − r e l a t i v ee rr o r Fig. 32. Relative error of the splitting schemeapproximation.
LIE STRANG SYMMETRIC2 SYMMETRIC4 SYMMETRIC6 SYMMETRIC8 ONV DIFF , n = 6400 and ˙ X ( t ) = A T X ( t ) + X ( t ) A − X ( t ) BB T X ( t ) + C T C, X (0) = 0. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − a b s o l u t ee rr o r Fig. 33. Absolute error of the splittingscheme approximation. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − r e l a t i v ee rr o r Fig. 34. Relative error of the splitting schemeapproximation. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − a b s o l u t ee rr o r Fig. 35. Absolute error of the splittingscheme approximation. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − r e l a t i v ee rr o r Fig. 36. Relative error of the splitting schemeapproximation. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − a b s o l u t ee rr o r Fig. 37. Absolute error of the splittingscheme approximation. .
02 0 .
04 0 .
06 0 .
08 0 . . − − − − − − − − − − − − − CONV DIFF , n = 6400, t ∈ [0 , . h = 2 − r e l a t i v ee rr o r Fig. 38. Relative error of the splitting schemeapproximation.