- Email: [email protected]
Monotonicity of quadratic forms with symplectic Runge—Kutta methods
APPLIED NUMERICAL K~.THEMATICS ELSEVIER
Applied Numerical Mathematics 17 (1995) 293-298
Monotonicity of quadratic forms with symplectic Runge-Kutta methods T i m o Eirola * Institute of Mathematics, Helsinki University of Technology, SF-02150 Espoo, Finland
1. Introduction
It is well known that the Gauss (Gauss-Legendre) methods are algebraically stable, and since this implies B-stability they produce contractive numerical maps when applied to contractive systems [5]. Also it is known that they (or more generally symplectic RKs) preserve any quadratic integrals of the system. In particular, they produce symplectic maps for Hamiltonian systems [6]. Further, in [2] it was noticed that they can be used to preserve positive-definiteness of the numerical solution of Riccati differential equations in the Hamiltonian form. In [7], Humphries and Stuart noticed connections between some of these results and, in [1], Bochev and Scovel showed that symplecticity is a consequence of preservation of quadratic forms and the commutativity of discretization and differentiation with respect to initial condition. Here we demonstrate that all these results follow from a general monotonicity property of these methods for quadratic forms.
2. The general result
Consider an autonomous differential equation
Jc(t) = f ( x ( t ) ) ,
x(t) ~ 12 c ~",
(2.1)
where f ~ ck(12), 12 open, and k >/1. Assume that the solutions can be extended to the whole time axis. For simplicity we assume that a nonautonomous equation is written in this form with one component of f being identically 1. Denote by ~b,, t ~ ~, the flow of f, i.e., O,~b,(x) =f(~b,(x)),
*
~b0(x ) =X.
E-mail: [email protected].
0168-9274/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0 1 6 8 - 9 2 7 4 ( 9 5 ) 0 0 0 3 4 - 8
Z Eirola / Applied Numerical .Mathematics 17 (1995) 293-298
294
Consider an s-stage R u n g e - K u t t a method, applied to (2.1):
dJh(x)=x +h ~-'~bif(pi),
p i = x +h ~ a o f ( p j ) .
t=l
(2.2)
j-I
We say that the RK step is well defined for h > 0, x • S2, if the R u n g e - K u t t a stages Pi (and consequently ~h(x)) are uniquely defined by the implicit function t h e o r e m via continuation from h = 0 , p , = x . We will call ~h the numerical flow, even though it is not a flow. R e m a r k 2.1. For given x let h ( x ) be the supremum of those h such that the R K step is well defined for h, x. Then it is easy to show that h is lower semicontinuous and it follows that h ( C ) is positive for any compact C c f2. If p ~ k is the order of the method, then for 0 < h ~
Vi,j = 1 . . . . . s,
(2.3)
and that the b/s are all nonnegative. The Gauss methods are such. Note that for algebraic stability the requirement on the matrix {m~j} is only that it is positive-semidefinite. Let 8" be a space of functions from a set S to ~" and g~a the subset of functions which have compact ranges in 22. Assume ~ has a topology such that oh,, when considered as a flow in g~a: r/• ~
~ d~, ° r / • ~
1 limt[&
t~0
o~-~]
is continuous,
(2.4)
=fo
for every 77 • ~2. R e m a r k 2.2. We defined the functions of ~ h(~7(S)) for every r / • ~ w
to have compact ranges in order to have positive
Let y be a topological vector space with a closed convex cone K (i.e., a/~ + b K c K for every a,b > 0), which we will call positive and for y,37 e y denote y ~ if y -3~ • K . In some of the applications, K will be the trivial cone containing just the origin. Let now /3 be a continuous bilinear map from ~ × g " to y , and q(~7):=/3(r/, r/) the corresponding quadratic map. We take this somewhat abstract setting in order to get many results as corollaries of one theorem.
Theorem 2.3. Assume that q is nondecreasing along the trajectories of Or in ~a, i.e., q(Ot°~7) ~ g ( ~ ) for any t > 0 and ~7 ~ ~s~. Then also for euery ~7 ~ ~
and h > 0 for which the RK step is well defined.
1~Eirola/Applied NumericalMathematics 17 (1995) 293-298
295
Proof. By continuity o f / 3 and by (2.4) we have for any ~7 ~ g~a O~lim-[q(4~,o~7)-q(rl)] t~0 l
= lira /3
t
t~0
,~b, orl
+/3 rl,
t
= / 3 ( f o r/, r/) +/3(rl, f ° rl).
(2.5)
The R u n g e - K u t t a m e t h o d produces: (2.6)
t~horl = rl + h ~ b i f o,ri-i, ,n-i= rl + h ~ a i j f o Trj. i
I
j=l
Thus we get:
i=l
j=l
i=1
s
,]
j~l
bibi/3(fo'n'i, foT"rj). i,j= 1
+h 2 ~
(2.7)
F r o m the latter of Eqs. (2.6) we get /3 ( f ° rr i , ~7) = fl ( f ° rr ~, 7ri ) - h ~ a J 3 ( f o 7ri , f o .n'j ) , j-1
/3(77, f o z r j ) = / 3 ( ~ i , f o Trj ) - h ~ aj,/3( f o er , f o zrj ). i-I
Putting these in (2.7) gives q(~bhorl):q(rl)+h
ff'_,bi[/3(foTr,,'n,-i)+/3(Tri, i=1
fozri) ] -h 2 ~
mij/3(foeri,
forrj).
i,/-1
(2.8) Now, by (2.3) the h2-terms vanish and the assertion follows from (2.5) and positivity of the bi's. [] R e m a r k 2.4. There is nothing really new in the proof above. T h e essential algebraic formula (2.8) from which the result can be read (and which holds for any R K method) can be seen in the proofs of all the results m e n t i o n e d in the beginning.
3. Applications In the following applications we either have S a singleton set, when ~ can be identified with ~n, or ~ = c r ( o , ~ ) , 0 <~ r <<,k, with the topology of uniform convergence on compact sets of derivatives up to order r. It is easy to show that (2.4) holds in these cases.
72 Eirola /Applied Numerical Mathematics 17 (1995) 293-298
296
R e m a r k 3.1. For B-stability, i.e., q being a negative-semidefinite quadratic form in ~n such that q(cbt(x) - (ot(x')) is nondecreasing for any two initial points x and x', taking ~ ~2n one sees from (2.8) (since the Kronecker product of two symmetric semidefinite matrices is semidefinite) that it suffices for {mij} to be positive-semidefinite. =
3.1. Conservation of quadratic integrals Assume that the flow satisfies
(4~t(x), Q~bt(x )) = (x, Qx),
i.e.,
xT(Q + o T ) f ( x ) = 0
for every t and x, where Q is a constant matrix. Take ~ = I~", /3(x, y) = (x, Qy), and set y = ~, K = {0}. Then the t h e o r e m says that (4'h(X), Q4'h(X)) -- ( X, Qx ) for every h and x for which the RK step is well defined, i.e., the quadratic form is preserved also by the numerical flow. A n application of this is that when one solves /J(t) = S(U(t))U(t) with these methods for an (n X n) matrix U in the situation, where S(U) is skew-symmetric for every U, then the numerical solution will (as well as the true one) always be an orthogonal matrix, provided the initial condition is (see [4]).
3.2. Symplecticity Assume that f is Hamiltonian, i.e., n is even and
f = J ° H, where J = [ 0
-I
~J
It follows that 4', is a symplectic map: O T for every t, x. Take ~ = C ' ( . Q , ~ ) and /3(~, r/)=~x(')TJr/x(-), and set y = C ( O , K = {0}. T h e n the t h e o r e m says that for every r/ in ~s~ we have
g(Ox6 )
=
~"×"),
Snx.
Thus (O~qSh(x))TJO~h(X) = J , i.e., ~h is also a symplectic map.
3.3. Positive-definiteness for Riccati equations Consider solving numerically the symmetric matrix Riccati equation
)((t) = A ( t ) X ( t ) + X ( t ) A ( t )T _ X ( t ) B ( t ) X ( t ) + C(t), (3.1) in ~d×d with a symmetric positive-semidefinite initial condition X(0) = X 0. We assume that the matrices B and C are symmetric and positive-semidefinite. T h e n the solution of (3.1) is symmetric and positive-semidefinite. It was shown in [2] that essentially only the backward Euler m e t h o d preserves this property when applied to (3.1).
T. Eirola /Applied Numerical Mathematics 17 (1995) 293-298
297
On the other hand it is well known that the solution of (3.1) can also be obtained as X ( t ) = Y ( t ) Z ( t ) -1, where Y(t) ]
[A(t)
C(t)
][Y(t) ]
[Y(0)
Denote the matrix above by M. To show the nice properties of Gauss methods for this equation take ~' = Nadxd, y = Ndxd, let K be the cone of positive-semidefinite (d x d) matrices, and set
Then
= (BY-ATz)TY+
Z T ( A y + CZ) = y T B y + z T c z ~ O,
i.e., q is increasing along the trajectories. This is then true also for the numerical trajectories (Yk, Zk). Hence Zkx Y~ and consequently X k := yk(zk) I = (Z kT) 1Z~T YkZk-1 will be positivesemidefinite (definite if the initial condition is). Similarly, taking /£ = 0 and
one easily shows that the numerical solution X~ :-- YkZ~ ~ is also symmetric.
3.4. MonotoniciO~ If in (3.1) one perturbs matrices B, C, X o (pointwise) to /} ~ B , C ~ C 2 0 ~ X 0 without changing A, then it is easy to show that the corresponding solution also satisfies X ( t ) ~ X ( t ) . The same property for the numerical solution is also true if one obtains it by applying a canonical scheme to (3.2). This too can be shown via a (less simple) application of the previous theorem (see [3]). R e m a r k 3.2 (Warning). In the beginning of this work the author expected that conformality would also be preserved by these methods. This means the property that T
(0x¢,(x)) Ox¢bt(x) = a(t, x ) I for some scalar function a. This holds if the Jacobian of f is of the form: skew-symmetric plus scalar times identity. This, however, is not a property that could be written as monotonicity of a quadratic map (unless a - 1) and the Gauss methods do not preserve it.
298
T. Eirola /Applied Numerical Mathematics 17 (1995) 293-298
References [I] P.B. Bochev and C. Scovel, On quadratic invariants and symplectic structure, BIT (1994). [2] L. Dieci and T. Eirola Positive definiteness in the numerical solution of Riccati differential equations, Numer. Math. 67 (1994) 303-313. [3] L. Dieci and T. Eirola, Preserving monotonicity in the numerical solution of Riccati differential equations (to appear). [4] L. Dieci, R.D. Russell and E.S. Van Vleck, Unitary integrators and applications to continuous orthonormalization techniques, SIAM. J. Numer. AnaL 31 (1994) 261-281. [5] E. Hairer and G. Wanner, Solving Ordinary Differential Equations H (Springer-Verlag, Berlin, Heidelberg 1991). [6] J.M. Sanz-Serna, Symplectic integrators for Hamiltonian problems: an overview, Acta Numer. 1 (1992) 243-286. [7] A.M. Stuart and A.R. Humphries, Model problems in numerical stability theory for initial value problems, SIAM ReLy. 36 (1994) 226-257.
Applied Numerical Mathematics 17 (1995) 293-298
Monotonicity of quadratic forms with symplectic Runge-Kutta methods T i m o Eirola * Institute of Mathematics, Helsinki University of Technology, SF-02150 Espoo, Finland
1. Introduction
It is well known that the Gauss (Gauss-Legendre) methods are algebraically stable, and since this implies B-stability they produce contractive numerical maps when applied to contractive systems [5]. Also it is known that they (or more generally symplectic RKs) preserve any quadratic integrals of the system. In particular, they produce symplectic maps for Hamiltonian systems [6]. Further, in [2] it was noticed that they can be used to preserve positive-definiteness of the numerical solution of Riccati differential equations in the Hamiltonian form. In [7], Humphries and Stuart noticed connections between some of these results and, in [1], Bochev and Scovel showed that symplecticity is a consequence of preservation of quadratic forms and the commutativity of discretization and differentiation with respect to initial condition. Here we demonstrate that all these results follow from a general monotonicity property of these methods for quadratic forms.
2. The general result
Consider an autonomous differential equation
Jc(t) = f ( x ( t ) ) ,
x(t) ~ 12 c ~",
(2.1)
where f ~ ck(12), 12 open, and k >/1. Assume that the solutions can be extended to the whole time axis. For simplicity we assume that a nonautonomous equation is written in this form with one component of f being identically 1. Denote by ~b,, t ~ ~, the flow of f, i.e., O,~b,(x) =f(~b,(x)),
*
~b0(x ) =X.
E-mail: [email protected].
0168-9274/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0 1 6 8 - 9 2 7 4 ( 9 5 ) 0 0 0 3 4 - 8
Z Eirola / Applied Numerical .Mathematics 17 (1995) 293-298
294
Consider an s-stage R u n g e - K u t t a method, applied to (2.1):
dJh(x)=x +h ~-'~bif(pi),
p i = x +h ~ a o f ( p j ) .
t=l
(2.2)
j-I
We say that the RK step is well defined for h > 0, x • S2, if the R u n g e - K u t t a stages Pi (and consequently ~h(x)) are uniquely defined by the implicit function t h e o r e m via continuation from h = 0 , p , = x . We will call ~h the numerical flow, even though it is not a flow. R e m a r k 2.1. For given x let h ( x ) be the supremum of those h such that the R K step is well defined for h, x. Then it is easy to show that h is lower semicontinuous and it follows that h ( C ) is positive for any compact C c f2. If p ~ k is the order of the method, then for 0 < h ~
Vi,j = 1 . . . . . s,
(2.3)
and that the b/s are all nonnegative. The Gauss methods are such. Note that for algebraic stability the requirement on the matrix {m~j} is only that it is positive-semidefinite. Let 8" be a space of functions from a set S to ~" and g~a the subset of functions which have compact ranges in 22. Assume ~ has a topology such that oh,, when considered as a flow in g~a: r/• ~
~ d~, ° r / • ~
1 limt[&
t~0
o~-~]
is continuous,
(2.4)
=fo
for every 77 • ~2. R e m a r k 2.2. We defined the functions of ~ h(~7(S)) for every r / • ~ w
to have compact ranges in order to have positive
Let y be a topological vector space with a closed convex cone K (i.e., a/~ + b K c K for every a,b > 0), which we will call positive and for y,37 e y denote y ~ if y -3~ • K . In some of the applications, K will be the trivial cone containing just the origin. Let now /3 be a continuous bilinear map from ~ × g " to y , and q(~7):=/3(r/, r/) the corresponding quadratic map. We take this somewhat abstract setting in order to get many results as corollaries of one theorem.
Theorem 2.3. Assume that q is nondecreasing along the trajectories of Or in ~a, i.e., q(Ot°~7) ~ g ( ~ ) for any t > 0 and ~7 ~ ~s~. Then also for euery ~7 ~ ~
and h > 0 for which the RK step is well defined.
1~Eirola/Applied NumericalMathematics 17 (1995) 293-298
295
Proof. By continuity o f / 3 and by (2.4) we have for any ~7 ~ g~a O~lim-[q(4~,o~7)-q(rl)] t~0 l
= lira /3
t
t~0
,~b, orl
+/3 rl,
t
= / 3 ( f o r/, r/) +/3(rl, f ° rl).
(2.5)
The R u n g e - K u t t a m e t h o d produces: (2.6)
t~horl = rl + h ~ b i f o,ri-i, ,n-i= rl + h ~ a i j f o Trj. i
I
j=l
Thus we get:
i=l
j=l
i=1
s
,]
j~l
bibi/3(fo'n'i, foT"rj). i,j= 1
+h 2 ~
(2.7)
F r o m the latter of Eqs. (2.6) we get /3 ( f ° rr i , ~7) = fl ( f ° rr ~, 7ri ) - h ~ a J 3 ( f o 7ri , f o .n'j ) , j-1
/3(77, f o z r j ) = / 3 ( ~ i , f o Trj ) - h ~ aj,/3( f o er , f o zrj ). i-I
Putting these in (2.7) gives q(~bhorl):q(rl)+h
ff'_,bi[/3(foTr,,'n,-i)+/3(Tri, i=1
fozri) ] -h 2 ~
mij/3(foeri,
forrj).
i,/-1
(2.8) Now, by (2.3) the h2-terms vanish and the assertion follows from (2.5) and positivity of the bi's. [] R e m a r k 2.4. There is nothing really new in the proof above. T h e essential algebraic formula (2.8) from which the result can be read (and which holds for any R K method) can be seen in the proofs of all the results m e n t i o n e d in the beginning.
3. Applications In the following applications we either have S a singleton set, when ~ can be identified with ~n, or ~ = c r ( o , ~ ) , 0 <~ r <<,k, with the topology of uniform convergence on compact sets of derivatives up to order r. It is easy to show that (2.4) holds in these cases.
72 Eirola /Applied Numerical Mathematics 17 (1995) 293-298
296
R e m a r k 3.1. For B-stability, i.e., q being a negative-semidefinite quadratic form in ~n such that q(cbt(x) - (ot(x')) is nondecreasing for any two initial points x and x', taking ~ ~2n one sees from (2.8) (since the Kronecker product of two symmetric semidefinite matrices is semidefinite) that it suffices for {mij} to be positive-semidefinite. =
3.1. Conservation of quadratic integrals Assume that the flow satisfies
(4~t(x), Q~bt(x )) = (x, Qx),
i.e.,
xT(Q + o T ) f ( x ) = 0
for every t and x, where Q is a constant matrix. Take ~ = I~", /3(x, y) = (x, Qy), and set y = ~, K = {0}. Then the t h e o r e m says that (4'h(X), Q4'h(X)) -- ( X, Qx ) for every h and x for which the RK step is well defined, i.e., the quadratic form is preserved also by the numerical flow. A n application of this is that when one solves /J(t) = S(U(t))U(t) with these methods for an (n X n) matrix U in the situation, where S(U) is skew-symmetric for every U, then the numerical solution will (as well as the true one) always be an orthogonal matrix, provided the initial condition is (see [4]).
3.2. Symplecticity Assume that f is Hamiltonian, i.e., n is even and
f = J ° H, where J = [ 0
-I
~J
It follows that 4', is a symplectic map: O T for every t, x. Take ~ = C ' ( . Q , ~ ) and /3(~, r/)=~x(')TJr/x(-), and set y = C ( O , K = {0}. T h e n the t h e o r e m says that for every r/ in ~s~ we have
g(Ox6 )
=
~"×"),
Snx.
Thus (O~qSh(x))TJO~h(X) = J , i.e., ~h is also a symplectic map.
3.3. Positive-definiteness for Riccati equations Consider solving numerically the symmetric matrix Riccati equation
)((t) = A ( t ) X ( t ) + X ( t ) A ( t )T _ X ( t ) B ( t ) X ( t ) + C(t), (3.1) in ~d×d with a symmetric positive-semidefinite initial condition X(0) = X 0. We assume that the matrices B and C are symmetric and positive-semidefinite. T h e n the solution of (3.1) is symmetric and positive-semidefinite. It was shown in [2] that essentially only the backward Euler m e t h o d preserves this property when applied to (3.1).
T. Eirola /Applied Numerical Mathematics 17 (1995) 293-298
297
On the other hand it is well known that the solution of (3.1) can also be obtained as X ( t ) = Y ( t ) Z ( t ) -1, where Y(t) ]
[A(t)
C(t)
][Y(t) ]
[Y(0)
Denote the matrix above by M. To show the nice properties of Gauss methods for this equation take ~' = Nadxd, y = Ndxd, let K be the cone of positive-semidefinite (d x d) matrices, and set
Then
= (BY-ATz)TY+
Z T ( A y + CZ) = y T B y + z T c z ~ O,
i.e., q is increasing along the trajectories. This is then true also for the numerical trajectories (Yk, Zk). Hence Zkx Y~ and consequently X k := yk(zk) I = (Z kT) 1Z~T YkZk-1 will be positivesemidefinite (definite if the initial condition is). Similarly, taking /£ = 0 and
one easily shows that the numerical solution X~ :-- YkZ~ ~ is also symmetric.
3.4. MonotoniciO~ If in (3.1) one perturbs matrices B, C, X o (pointwise) to /} ~ B , C ~ C 2 0 ~ X 0 without changing A, then it is easy to show that the corresponding solution also satisfies X ( t ) ~ X ( t ) . The same property for the numerical solution is also true if one obtains it by applying a canonical scheme to (3.2). This too can be shown via a (less simple) application of the previous theorem (see [3]). R e m a r k 3.2 (Warning). In the beginning of this work the author expected that conformality would also be preserved by these methods. This means the property that T
(0x¢,(x)) Ox¢bt(x) = a(t, x ) I for some scalar function a. This holds if the Jacobian of f is of the form: skew-symmetric plus scalar times identity. This, however, is not a property that could be written as monotonicity of a quadratic map (unless a - 1) and the Gauss methods do not preserve it.
298
T. Eirola /Applied Numerical Mathematics 17 (1995) 293-298
References [I] P.B. Bochev and C. Scovel, On quadratic invariants and symplectic structure, BIT (1994). [2] L. Dieci and T. Eirola Positive definiteness in the numerical solution of Riccati differential equations, Numer. Math. 67 (1994) 303-313. [3] L. Dieci and T. Eirola, Preserving monotonicity in the numerical solution of Riccati differential equations (to appear). [4] L. Dieci, R.D. Russell and E.S. Van Vleck, Unitary integrators and applications to continuous orthonormalization techniques, SIAM. J. Numer. AnaL 31 (1994) 261-281. [5] E. Hairer and G. Wanner, Solving Ordinary Differential Equations H (Springer-Verlag, Berlin, Heidelberg 1991). [6] J.M. Sanz-Serna, Symplectic integrators for Hamiltonian problems: an overview, Acta Numer. 1 (1992) 243-286. [7] A.M. Stuart and A.R. Humphries, Model problems in numerical stability theory for initial value problems, SIAM ReLy. 36 (1994) 226-257.