Previous Page Next Page Contents

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

Option: method

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

Changes




Do you have questions or comments?


Copyright © SciFace Software GmbH & Co. KG 2000