numeric::expMatrix
-- the
exponential of a matrixnumeric::expMatrix
(A, ..)
returns the
exponential exp(A) of a square matrix A.
numeric::expMatrix
(A, x, ..)
with a vector
x returns the vector exp(A)*x.
numeric::expMatrix
(A, X, ..)
with a matrix
X returns the matrix exp(A)*X.
numeric::expMatrix(A <, method>)
numeric::expMatrix(A, x <, method>)
numeric::expMatrix(A, X <, method>)
A |
- | a numerical square matrix of domain type DOM_ARRAY or of category Cat::Matrix |
x |
- | a vector represented by a list
[x[1],..,x[n]] or a 1-dimensional array
array(1..n,[x[1],..,x[n]]) |
X |
- | an n x m matrix of domain type DOM_ARRAY or Dom::Matrix(Ring) with a suitable
coefficient ring Ring |
method |
- | specifies the numerical method used for computing the result. The available methods are Diagonalization, Interpolation, Krylov and TaylorExpansion. |
All results are float matrices/vectors.
The call numeric::expMatrix
(A <,
method>)
returns exp(A) as a matrix of domain type
DOM_ARRAY
.
The call numeric::expMatrix
(A, x, <,
method>)
returns exp(A) x as a vector of
the same domain type as the input vector x
, i.e., either
as a list or as a 1-dimensional array
array(1..n,[..])
.
The call numeric::expMatrix
(A, X, <,
method>)
returns exp(A) X as an n x
m matrix of domain type DOM_ARRAY
.
The function is sensitive to the environment variable DIGITS
, which determines the
numerical working precision.
float
. Numerical symbolic
expressions such as PI
, sqrt(2)
,
exp(-1)
etc. are accepted. They are converted to
floats.numeric::expMatrix
(A,
x)
/numeric::expMatrix
(A, X)
is
faster.10^(-DIGITS)
, unless numerical roundoff prevents
reaching this precision goal. Roughly speaking: all digits of all
components of the result are reliable up to roundoff effects.The methods Diagonalization, Interpolation and Krylov compute the result to a relative precision w.r.t. the norm:
norm(error) <= 10^(-DIGITS)*norm(result) .Consequently, if the result has components of different orders of magnitude, then the smaller components have larger relative errors than the large components. Not all digits of the small components are reliable! Cf. example 2.
The method Diagonalization
only works for diagonalizable matrices. For matrices without a basis of
eigenvectors numeric::expMatrix
may either produce an
error or the returned result is dominated by roundoff effects. For
symmetric/Hermitean or skew/skew-Hermitean matrices this method
produces reliable results.
The method Interpolation may become numerically unstable for certain matrices. The algorithm tries to detect such unstabilities and stops with an error message.
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.
We consider the matrix
>> A := array(1..2, 1..2, [[1, 0] , [1, PI]]):
>> expA := numeric::expMatrix(A)
+- -+ | 2.718281829, 0 | | | | 9.536085572, 23.14069263 | +- -+
We consider a vector given by a list x1
and
by an equivalent 1-dimensional array x2
,
respectively:
>> x1 := [1, 1]: x2 := array(1..2, [1, 1]):
The constructor M
of the matrix domain
Dom::Matrix()
converts the
list x1
to an equivalent 2 x 1 matrix
X
(a column):
>> M := Dom::Matrix(): X := M(x1):
The first call yields a list, the second a 1-dimensional array, the third a 2 x 1 array:
>> numeric::expMatrix(A, x1), numeric::expMatrix(A, x2, Krylov), numeric::expMatrix(A, X, Diagonalization)
+- -+ [2.718281829, 32.6767782], | 2.718281829, 32.6767782 |, +- -+ +- -+ | 2.718281829 | | | | 32.6767782 | +- -+
For further processing the array expA
is
converted to an element of the matrix domain:
>> expA := M(expA):
Now the overloaded arithmetical operators
+
, *
, ^
etc. can be used for
further computations:
>> expA*X
+- -+ | 2.718281829 | | | | 32.6767782 | +- -+
>> delete A, expA, x1, x2, M, X:
We demonstrate the different precision goals of the methods.
>> A := array(1..3, 1..3, [[ 1000, 1, 0 ], [ 0, 1, 1 ], [1/10^100, 0, -1000]]):
The default method TaylorExpansion computes each component of exp(A) correctly:
>> numeric::expMatrix(A)
+- -+ | 1.970071114e434, 1.972043157e431, 9.860215786e427 | | | | 9.860215786e327, 9.870085872e324, 4.935042936e321 | | | | 9.85035557e330, 9.860215786e327, 4.930107893e324 | +- -+
The method Diagonalization produces a result, which is accurate in the sense that norm(error) <= 10^(-DIGITS)* norm(exp(A)) holds. Indeed, the largest components of exp(A) are correct. However, Diagonalization does not even get the right order of magnitude of the smaller components:
>> numeric::expMatrix(A, Diagonalization)
+- -+ | 1.970071114e434, 1.972043157e431, 0 | | | | 0, 2.718281829, 0 | | | | 0, 0, 5.075958898e-435 | +- -+
Note that exp(A) 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)
+- -+ | 1.970071114e434, 1.972043157e431, 9.860215786e427 | | | | 0, 2.718281829, 0.002715566262 | | | | 0, 0, 5.075958897e-435 | +- -+
>> numeric::expMatrix(B, Diagonalization)
+- -+ | 1.970071114e434, 1.972043157e431, 9.860215786e427 | | | | 0.0, 2.718281829, 0.002715566262 | | | | 0.0, 0.0, 5.075958898e-435 | +- -+
>> delete A, B:
Hilbert matrices H[i,j]=1/(i+j-1) have real positive eigenvalues. For large dimension most of these eigenvalues are small and may be regarded as a single cluster. Consequently, the option Krylov is useful:
>> numeric::expMatrix(linalg::hilbert(100), [1 $ 100], Krylov)
[25.47080919, 18.59337041, ... , 2.863083064, 2.848538965]
exp(A) = 1 + A + A^2/2 + ..in a suitable numerically stable way.
numeric::expMatrix
uses polynomial arithmetic to
multiply matrices and vectors. Thus sparse matrices are handled
efficiently based on MuPADs internal sparse representation of
polynomials.