"OrthogonalProjection" Method for NDSolve
Introduction
Consider the matrix differential equation
where the initial value 
is given. Assume that 
, that the solution has the property of preserving orthonormality, 
, and that it has full rank for all 
.
From a numerical perspective, a key issue is how to numerically integrate an orthogonal matrix differential system in such a way that the numerical solution remains orthogonal. There are several strategies that are possible. One approach, suggested in [DRV94], is to use an implicit Runge–Kutta method (such as the Gauss scheme). Some alternative strategies are described in [DV99] and [DL01].
The approach taken here is to use any reasonable numerical integration method and then postprocess using a projective procedure at the end of each integration step.
An important feature of this implementation is that the basic integration method can be any built-in numerical method, or even a user-defined procedure. In the following examples an explicit Runge–Kutta method is used for the basic time stepping. However, if greater accuracy is required an extrapolation method could easily be used, for example, by simply setting the appropriate Method option.
Projection Step
At the end of each numerical integration step you need to transform the approximate solution matrix of the differential system to obtain an orthogonal matrix. This can be carried out in several ways (see for example [DRV94] and [H97]):
The Newton and Schulz methods are quadratically convergent, and the number of iterations may vary depending on the error tolerances used in the numerical integration. One or two iterations are usually sufficient for convergence to the orthonormal polar factor (see the following) in IEEE double-precision arithmetic.
QR decomposition is cheaper than singular value decomposition (roughly by a factor of two), but it does not give the closest possible projection.
Definition (Thin singular value decomposition [GVL96]): Given a matrix 
with 
there exist two matrices 
and 
such that 
is the diagonal matrix of singular values of 
, 
, where 
. 
has orthonormal columns and 
is orthogonal.
Definition (Polar decomposition): Given a matrix 
and its singular value decomposition 
, the polar decomposition of 
is given by the product of two matrices 
and 
where 
and 
. 
has orthonormal columns and 
is symmetric positive semidefinite.
The orthonormal polar factor 
of 
is the matrix that solves
for the 2 and Frobenius norms [H96].
Schulz Iteration
The approach chosen is based on the Schulz iteration, which works directly for 
. In contrast, Newton iteration for 
needs to be preceded by QR decomposition.
Comparison with direct computation based on the singular value decomposition is also given.
The Schulz iteration is given by
The Schulz iteration has an arithmetic operation count per iteration of 
floating-point operations, but is rich in matrix multiplication [H97].
In a practical implementation, GEMM-based level 3 BLAS of LAPACK [LAPACK99] can be used in conjunction with architecture-specific optimizations via the Automatically Tuned Linear Algebra Software [ATLAS00]. Such considerations mean that the arithmetic operation count of the Schulz iteration is not necessarily an accurate reflection of the observed computational cost. A useful bound on the departure from orthonormality of 
is in [H89]: 
. Comparison with the Schulz iteration gives the stopping criterion 
for some tolerance 
.
Standard Formulation
Assume that an initial value 
for the current solution of the ODE is given, together with a solution 
from a one-step numerical integration method. Assume that an absolute tolerance 
for controlling the Schulz iteration is also prescribed.
The following algorithm can be used for implementation.
Increment Formulation
NDSolve uses compensated summation to reduce the effect of rounding errors made by repeatedly adding the contribution of small quantities 
to 
at each integration step [H96]. Therefore, the increment 
is returned by the base integrator.
An appropriate orthogonal correction 
for the projective iteration can be determined using the following algorithm.
This modified algorithm is used in "OrthogonalProjection" and shows an advantage of using an iterative process over a direct process, since it is not obvious how an orthogonal correction can be derived for direct methods.
Examples
Orthogonal Error Measurement
of a matrix 
can be defined in terms of the Norm function as follows:
can then be measured using this function [H97]:Square Systems
This example concerns the solution of a matrix differential system on the orthogonal group 
(see [Z98]).
The matrix differential system is given by
are 
, 
, 
. Thus as 
approaches 
, two of the eigenvalues of 
approach -1. The numerical integration is carried out on the interval 
:
Rectangular Systems
In the following example it is shown how the implementation of the orthogonal projection method also works for rectangular matrix differential systems. Formally stated, the interest is in solving ordinary differential equations on the Stiefel manifold, the set of 
×
orthogonal matrices with 
.
Definition The Stiefel manifold of 
×
orthogonal matrices is the set 
, 
, where 
is the 
×
identity matrix.
Solutions that evolve on the Stiefel manifold find numerous applications such as eigenvalue problems in numerical linear algebra, computation of Lyapunov exponents for dynamical systems and signal processing.
Consider an example adapted from [DL01]:
The exact solution is given by
it follows that 
satisfies the following weak skew-symmetric system on 
:
Implementation
The implementation of the method "OrthogonalProjection" has three basic components:
- Initialization. Set up the base method to use in the integration, determining any method coefficients and setting up any workspaces that should be used. This is done once, before any actual integration is carried out, and the resulting MethodData object is validated so that it does not need to be checked at each integration step. At this stage the system dimensions and initial conditions are checked for consistency.
- Perform an orthogonal projection. This performs various tests such as checking that the basic integration proceeded correctly and that the Schulz iteration converges.
Options can be used to modify the stopping criteria for the Schulz iteration. One option provided by the code is "IterationSafetyFactor" which allows control over the tolerance 
of the iteration. The factor is combined with a Unit in the Last Place, determined according to the working precision used in the integration (
for IEEE double precision).
The Frobenius norm used for the stopping criterion can be computed efficiently using the LAPACK LANGE functions [LAPACK99].
The option MaxIterations controls the maximum number of iterations that should be carried out.
Option Summary
option name | default value | |
| Dimensions | {} | specify the dimensions of the matrix differential system |
| "IterationSafetyFactor" |  | specify the safety factor to use in the termination criterion for the Schulz iteration (1) |
| MaxIterations | Automatic | specify the maximum number of iterations to use in the Schulz iteration (2) |
| Method | "StiffnessSwitching" | specify the method to use for the numerical integration |