Time-stepping and preserving orthonormality (Q678209)

Certain applications (for example, the computation of Lyapunov exponents or the construction of smooth singular value decompositions of parametrized matrices) lead to initial value problems for ordinary differential equations whose (matrix) solution preserves orthonormality. While for some special problem classes there exist time-stepping methods that have the same property, one will often have to resort to using a standard integrator and regularly replace the resulting approximate solution by an orthonormal matrix. The optimal replacement (closest with respect to the Frobenius norm) is given by the orthonormal polar factor \(U\) in \(A= UH\) (where \(U^TU= I\) and \(H\) is symmetric positive definite). The author analyzes this alternative to the more standard approach based on using the orthonormal factor \(Q\) in the QR decomposition of the solution [see, e.g., \textit{L. Dieci}, \textit{R. D. Russell} and \textit{E. S. van Vleck}, SIAM J. Numer. Anal. 31, No. 1, 261-281 (1994; Zbl 0815.65096)], with the focus being on finite-time error bounds and adaptive time-stepping. The numerical implementation and, in particular, the effect of the projection on the global error and the departure from orthogonality are discussed in detail.