## `numeric::fMatrix` -- functional calculus for numerical square matrices

#### Introduction

`numeric::fMatrix``(f, A, ..)` computes the matrix f(A,..) with a function f and a square matrix A.

#### Call(s)

`numeric::fMatrix(f, A, p1, p2, ..)`

#### Parameters

`f` | - | a procedure representing a scalar function f: C -> C or f: C x P x P x .. -> C, where P is a set of parameters
`A` | - | a square matrix of domain type `DOM_ARRAY` or of category `Cat::Matrix`
`p1, p2, ..` | - | arbitrary MuPAD objects accepted by f as parameters

#### Returns

The matrix f(A, ..) is returned as an array.

#### Side Effects

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

#### Related Functions

`numeric::expMatrix`, `numeric::inverse`

#### Details

- The components of 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 matrix A must be diagonalizable. `numeric::fMatrix` aborts with an error message, if it detects numerically that A is not diagonalizable. For most non-diagonalizable matrices, however, the numerical algorithm fails to detect this fact and the returned matrix is dominated by roundoff effects. It is the user's responsibility to make sure that the diagonalization is feasible and well conditioned.

- Symmetric/Hermitean and skew/skew Hermitean matrices can always be diagonalized in a numerically stable way. `numeric::fMatrix` produces reliable numerical results for such matrices.
- The procedure f must accept complex floating point numbers as first argument. It may return arbitrary MuPAD expressions, provided these can be multiplied with floating point numbers.
- The parameters `p1, p2, ..` may be numerical or symbolic objects. They must be accepted by `f` as 2nd argument, 3rd argument etc.
- In contrast to the components of A, numerical symbolic objects such as `PI`, `sqrt(2)` etc. passed as parameters `p1, p2,..` are not converted to floats.
- Inversion or exponentiation of a matrix may be realized with the functions f: a -> 1/a and f = exp, respectively. However, it is recommended to use the specialized algorithms `numeric::inverse` and `numeric::expMatrix` instead. Also matrix evaluation of low degree polynomials should be done with standard matrix arithmetic rather than with `numeric::fMatrix`.

#### Example 1

We compute the matrix power A^100:

    >> A := array(1..2, 1..2, [[2, PI], [exp(-10), 0]]):

    >> numeric::fMatrix(x -> x^100, A)

                       +-                                -+
                       |  1.272133133e30, 1.998190806e30  |
                       |                                  |
                       |  2.887634784e25, 4.535724387e25  |
                       +-                                -+

Alternatively you may use the function `_power` which takes the exponent as a second parameter.

    >> numeric::fMatrix(_power, A, 100)

    >> delete A:

#### Example 2

We compute the square root of a matrix:

    >> A := array(1..2, 1..2, [[0, 1], [-1, 0]]):

    >> B := numeric::fMatrix(sqrt, A)

                  array(1..2, 1..2,
                    (1, 1) = 0.7071067812 - 1.084202173e-19 I,
                    (1, 2) = 0.7071067812 + 2.710505431e-20 I,
                    (2, 1) = - 0.7071067812 - 1.084202173e-19 I,
                    (2, 2) = 0.7071067812 - 5.421010863e-20 I
                  )

The small imaginary parts are caused by numerical roundoff. We eliminate them by extracting the real parts of the components:

    >> B := map(B, Re)

                         +-                             -+
                         |   0.7071067812, 0.7071067812  |
                         |                               |
                         |  -0.7071067812, 0.7071067812  |
                         +-                             -+

We verify that `B^2` matrix is `A`. For convenience we convert `B` to an element of a matrix domain and compute the square by the overloaded operator `^`:

    >> B := Dom::Matrix(Dom::Complex)(B): B^2

                      +-                                   -+
                      |  5.421010863e-20,        1.0        |
                      |                                     |
                      |        -1.0,      -1.084202173e-19  |
                      +-                                   -+

This coincides with `A` up to numerical roundoff.

    >> delete A, B:

#### Example 3

We compute exp(t*PI*A) with a symbolic parameter t:

    >> A := array(1..2,1..2,[[0,1],[-1,0]]):

    >> numeric::fMatrix(exp@_mult, A, t*PI)

            array(1..2, 1..2,
              (1, 1) = 0.5 exp(-1.0 I t PI) + 0.5 exp(1.0 I t PI),
              (1, 2) = 0.5 I exp(-1.0 I t PI) - 0.5 I exp(1.0 I t PI),
              (2, 1) = 0.5 I exp(1.0 I t PI) - 0.5 I exp(-1.0 I t PI),
              (2, 2) = 0.5 exp(-1.0 I t PI) + 0.5 exp(1.0 I t PI)
            )

    >> delete A:

#### Background

- A numerical diagonalization A = X * diag(l1,l2,..) * X^(-1) is computed. The columns of X are the (right) eigenvectors of A, the diagonal entries l1,l2,.. are the corresponding eigenvalues. The function f is mapped to the eigenvalues, the matrix result is computed by

              f(A, p1, p2, ..) = X * diag(f(l1, p1, p2, ..),f(l2, p1, p2, ..), ..) * X^(-1).

  The eigenvector matrix X may be obtained via `numeric::eigenvectors(A)[2].`
- The condition number norm(X)*norm(X^(-1)) of the eigenvector matrix is a measure indicating how well conditioned the diagonalization of the matrix A is. If this number is larger than 10^DIGITS, then not a single digit of the diagonalization data is trustworthy.
- The call `numeric::fMatrix(exp, A)` is equivalent to

  `numeric::expMatrix(A, Diagonalization)`

#### Changes

- new function

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