## `numeric::expMatrix` -- the exponential of a matrix

#### Introduction

`numeric::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.

#### Call(s)

`numeric::expMatrix(A <, method>)`
`numeric::expMatrix(A, x <, method>)`
`numeric::expMatrix(A, X <, method>)`

#### Parameters

`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`

#### Options

`method` | - | specifies the numerical method used for computing the result. The available methods are Diagonalization, Interpolation, Krylov and TaylorExpansion.

#### Returns

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`.

#### Side Effects

The function is sensitive to the environment variable `DIGITS`, which determines the numerical working precision.

#### Related Functions

`exp`

#### Details

- The components of the matrix A must not contain symbolic objects which cannot be converted to numerical values via `float`. Numerical symbolic expressions such as `PI`, `sqrt(2)`, `exp(-1)` etc. are accepted. They are converted to floats.
- The methods Diagonalization and Interpolation do not work for all matrices (see below).
- Special algorithms are implemented for traceless 2 x 2 matrices and skew symmetric 3 x 3 matrices. Specification of a particular method does not have any effect on such matrices.
- If exp(A)*x or exp(A)*X is required, then you should not compute exp(A) first and then multiply the resulting matrix with the vector/matrix x/X. In general the call `numeric::expMatrix``(A, x)`/`numeric::expMatrix``(A, X)` is faster.

#### Option: method

- The method TaylorExpansion is the default algorithm. It produces fast results for matrices with small norms.
- The default method TaylorExpansion computes each individual component of exp(A), exp(A)*x, exp(A)*X to a relative precision of about `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.

- The method Krylov is only available for computing exp(A)*x with a vector x. Also vectors represented by n x 1 matrices are accepted.

  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.

#### Example 1

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:

#### Example 2

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:

#### Example 3

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]

#### Background

- The method TaylorExpansion sums the usual Taylor series

      exp(A) = 1 + A + A^2/2 + ..

  in a suitable numerically stable way.
- The method Diagonalization computes exp(A) = T * diag(exp(l1),exp(l2),..) * T^(-1) by a diagonalization A = T * diag(l1,l2,..) * T^(-1).
- The method Interpolation computes a polynomial P interpolating the function exp at the eigenvalues of A. Evaluation of the matrix polynomial yields exp(A)= P(A).
- The method Krylov reduces A to a Hessenberg matrix H and computes an approximation of exp(A)*x from exp(H). Depending on A and x the dimension of H may be smaller than the dimension of A. Reference: Y. Saad, ``Analysis of some Krylov Subspace Approximations to the Matrix Exponential Operator'', SIAM Journal of Numerical Analysis 29 (1992).
- `numeric::expMatrix` uses polynomial arithmetic to multiply matrices and vectors. Thus sparse matrices are handled efficiently based on MuPADs internal sparse representation of polynomials.

#### Changes

- new function

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000