JOURNAL OF COMPUTATIONAL AND APPLIED MATHEMATICS
ELSEVIER
Journal of Computational and Applied Mathematics 82 (1997) 261-278
Variable step-size techniques in continuous Runge-Kutta methods for isospectral dynamical systems I L. Lopez a, C. Mastroserio a, T. Politi b'* a Dipartimento di Matematica, Universitd di Bari, Via E. Orabona 4, 1-70125 Bari, Italy b Dipartimento di Matematica, Politecnico di Bari, Via E. Orabona 4, 1-70125 Bari, Italy Received 8 August 1996; received in revised form 20 January 1997
Abstract
In this paper we consider numerical methods for the dynamical system
L'=[B(L),L],
L(O) =Lo,
(*)
where L0 is a n x n symmetric matrix, [B(L),L] is the commutator of B(L) and L, and B(L) is a skew-symmetric matrix for each symmetric matrix L. The differential system is isospectral, i.e., L(t) preserves the eigenvalues of L0, for t~>0. The matrix B(L) characterizes the flow, and for special B(.), the solution matrix L(t) tends, as t increases, to a diagonal matrix with the same eigenvalues of L0. In [11] a modification of the MGLRK methods, introduced in [2], has been proposed. These procedures are based on a numerical approximation of the Flaschka formulation of ( . ) by Runge-Kutta (RK) methods. Our numerical schemes (denoted by EdGLRKs) consist in solving the system ( . ) by a continuous explicit Runge-Kutta method (CERK) and then performing a single step of a Gauss-Legendre RK method, for the Flaschka formulation of (*), in order to convert the approximation of L(t) to an isospectral solution. The problems of choosing a constant time step or a variable time step strategy are both of great importance in the application of these methods. In this paper, we introduce a definition of stability for the isospectral numerical methods. This definition involves a potential function associated to the isospectral flow. For the class EdGLRKs we propose a variable step-size strategy, based on this potential function, and an optimal constant time step h in the stability interval. The variable time step strategy will be compared with a known variable step-size strategy for RK methods applied to these dynamical systems. Numerical tests will be given and a comparison with the QR algorithm will be shown.
Keywords." Isospectral flows; Variable step-size; Continuous Runge-Kutta methods A M S classification." 65L06; 65L15
*
Corresponding author. E-mail:
[email protected]. Work supported, in part, by the Italian C.N.R. contract no. 96.00240.CT01.
0377-0427/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved PH S 03 7 7 - 0 4 2 7 ( 9 7 ) 0 0 0 4 8 - 4
262
L. Lopez et al./ Journal of Computational and Applied Mathematics 82 (1997) 261-278
1. Introduction to the problem
In the last few years a certain effort has been devoted to the numerical solution of isospectral flows (see [2, 3, 5, 9, 11, 12]). The usual form of these dynamical systems is given by L'(t) = [B(L(t)),L(t)],
L(O) =Lo,
t>~O,
(1)
where L0 is a given n x n symmetric real matrix, B(.) is Lipschitz matrix function mapping symmetric matrices into skew-symmetric matrices and [., .] is the standard Lie-bracket operator [B(L ), L] = B(L )L - LB(L ).
(2)
The theoretical solution L(t) of (1) and (2) is a symmetric real matrix preserving the eigenvalues, that is a ( L ( t ) ) = a(L0) for t~>0, where a(.) is the spectrum (see [3, 7, 10, 14]). The flow (1), (2) is said to be isospectral. Special cases are the Toda flow and the double-bracket flow. In the former case, L(t) is tridiagonal and B(L) = L+ - (L+)T,
(3)
where L+ is the strictly upper triangular part of the matrix L(t). When t --+ + cxz, the flow L(t) tends to a diagonal matrix L(oo), consisting of the eigenvalues of L0 and the dynamics of (1) is strictly related to the QR algorithm. The double-bracket flow, introduced by Brockett in [1], may be derived by (1) assuming
(4)
B ( L ) = [N,L],
where N is a fixed n x n target symmetric matrix. When N is a diagonal matrix with distinct eigenvalues, then L(t) tends, as t increases, to a diagonal matrix L(oo). If we set N=diag(n,n-
1,..., 1),
(5)
the double-bracket flow becomes equivalent to the Toda flow (see [2]). Then, in the following, we will only consider double-bracket flows, since Toda flows may be included in this class using (5). Recently, numerical isospectral methods, based on the Flaschka formulation of (1), have been proposed, see, e.g., the MGLRK methods [2] and the EdGLRKs methods [11]. In the application of these methods both the problems of choosing a constant time step and a variable time-step strategy are of great importance. For instance, an unsuitable constant time step could provide numerical instability, while a step-size variable version of these methods will be effective either to compare these methods to the existing numerical methods for finding eigenvalues or to follow the numerical solution, on the whole integration interval, accurately. In this paper, we introduce a definition of numerical stability for isospectral methods based on a potential function associated to the isospectral double-bracket flows. For the class EdGLRKs, we propose a variable step-size strategy, derived by this potential function, and an optimal constant time-step h for which the numerical stability is satisfied. The variable time-step strategy will be compared with a known variable step-size strategy for RK methods. Finally, numerical tests will be given and a comparison with the QR algorithm will be shown.
L. Lopez et al./ Journal of Computationaland Applied Mathematics 82 (1997) 261-278
263
2. Introduction to the numerical methods
In the Flaschka formulation of (1), the solution L(t) is written as
L(t) = U(t)LoUT(t),
t~>O,
(6)
with U(t) solution of the differential system
U'=B(L)U,
U(0) = I,
t~>O
(7)
where, since B(L) is skew-symmetric, U(t) is an unitary matrix for t~>0 (see [3, 4, 7, 14]). In particular, if the sequence tk+l = tk +hk, for k = 0,...,K, discretizes the time interval, then a computational form of the Flaschka formulation of (1) consists of a numerical approximation of (6) by T Lk+l = Vk+lLkVk+1,
(8)
where Lk+l ~--L(tk+~) and Vk+~ are unitary matrices, providing an approximation of U(t) at tk+l, obtained by Gauss-Legendre RK methods applied to the unitary differential system
U'(t)=B(L(t))U(t),
U(tk)=L
tE[tk,tk+~], k>~O
(see [6]). In [11], Lopez and Politi have proposed a modification of the isospectral MGLRK methods derived in [2]. It consists in solving
L'(t) = [B(L(t)),L(t)],
L(tk) =L~,
tE[tk, tk+,],
(9)
by a continuous v-stage explicit RK method (CERK) defined by C*
A*
(10)
b'T(0)
with c* =(c*,... ,cv)**, b*(O)=(b*(O),.. . ,by*(0))V, for 0~
