Exponential of a matrix
This functionality does not run in MATLAB.
numeric::expMatrix(A) returns the exponential
a square matrix A.
numeric::expMatrix(A, x) with a vector x returns
numeric::expMatrix(A, X) with a matrix X returns
If no return type is specified via the option
= d, the domain type of the result depends on the type of
the input matrix
For a dense matrixA of
Dom::DenseMatrix(Ring), the result is again
a matrix of type
Dom::DenseMatrix() over the ring
For all other matrices A of
the result is returned as a
Dom::Matrix() over the ring of
expressions. This includes input matrices
The components of A must
not contain symbolic objects which cannot be converted to numerical
Numerical symbolic expressions such as π,
are accepted. They are converted to floats.
The specification of a method such as
SoftwareFloats, i.e., the result is computed
via the software arithmetic of the MuPAD® kernel.
not work for all matrices (see below).
SoftwareFloats, special algorithms are
implemented for traceless 2×2 matrices
and skew symmetric 3×3 matrices.
Specification of a particular method does not have any effect for
required, one should not compute
and then multiply the resulting matrix with the vector/matrix x/X.
In general, the call
numeric::expMatrix(A, x) or
X), respectively, is faster.
The function is sensitive to the environment variable
which determines the numerical working precision.
We consider a lower triangular matrix given by an array:
A := array(1..2, 1..2, [[1, 0] , [1, PI]]): expA := numeric::expMatrix(A)
We consider a vector given by a list
by an equivalent 1-dimensional
x1 := [1, 1]: x2 := array(1..2, [1, 1]):
Further, an equivalent input vector X of
X := matrix(x1):
The following three calls all yield a vector represented by an 2×1 array corresponding to the type of the input matrix A:
numeric::expMatrix(A, x1), numeric::expMatrix(A, x2, Krylov), numeric::expMatrix(A, X, Diagonalization)
For further processing, the array
converted to an element of the matrix domain
expA := matrix(expA):
Now, the overloaded arithmetical operators
can be used for further computations:
delete A, expA, x1, x2, X:
We demonstrate the different precision goals of the methods. Note that software arithmetic is used when a method is specified:
A := array(1..3, 1..3, [[ 1000, 1, 0 ], [ 0, 1, 1 ], [1/10^100, 0, -1000]]):
The default method
TaylorExpansion computes each
Diagonalization produces a result,
which is accurate in the sense that
Indeed, the largest components of
Diagonalization does not even
get the right order of magnitude of the smaller components:
Note that is very sensitive to small changes in A. After elimination of the small lower triangular element, both methods yield the same result with correct digits for all entries:
B := array(1..3, 1..3, [[ 1000, 1, 0 ], [ 0 , 1, 1 ], [ 0 , 0, -1000]]): numeric::expMatrix(B, SoftwareFloats)
delete A, B:
real positive eigenvalues. For large dimension, most of these eigenvalues
are small and may be regarded as a single cluster. Consequently, the
Krylov is useful:
numeric::expMatrix(linalg::hilbert(100), [1 $ 100], Krylov)
One of the flags
One of the flags
Compared to the
If the result cannot be computed with hardware floats, software arithmetic by the MuPAD kernel is tried.
If the current value of
There may be several reasons for hardware arithmetic to fail:
The trailing digits in floating-point results computed with
The specification of a method implies
The default method
This method is fast when x is spanned by few eigenvectors of A. Further, if A has only few clusters of similar eigenvalues, then this method can be much faster than the other methods. Cf. Example 3.
Option, specified as
All results are float matrices/vectors. For an n×n matrix A:
numeric::expMatrix(A, method) returns
an n×n matrix,
numeric::expMatrix(A, x, method) returns
an n×1 matrix,
numeric::expMatrix(A, X, method) returns
an n×m matrix.
The domain type of the result depends on the domain type of
the input matrix A unless
a return type is requested explicitly via
TaylorExpansion sums the usual
in a suitable numerically stable way.
a diagonalization A = T diag(λ1, λ2,
…) T- 1.
Interpolation computes a polynomial P interpolating
the function exp at
the eigenvalues of A.
Evaluation of the matrix polynomial yields
Krylov reduces A to
a Hessenberg matrix H and
computes an approximation of
Depending on A and x,
the dimension of H may
be smaller than the dimension of A.
numeric::expMatrix uses polynomial arithmetic
to multiply matrices and vectors. Thus, sparse matrices are handled
efficiently based on the MuPAD internal sparse representation
Y. Saad, "Analysis of some Krylov Subspace Approximations to the Matrix Exponential Operator", SIAM Journal of Numerical Analysis 29 (1992).