Functional calculus for numerical square matrices
This functionality does not run in MATLAB.
numeric::fMatrix(f, A, p1, p2, …, options)
numeric::fMatrix(f, A) computes the matrix f(A) with a function f and a square matrix A.
If no return type is specified via the option ReturnType = d, the domain type of the result depends on the type of the input matrix A:
For a dense matrix A of type Dom::DenseMatrix() the result is a dense matrix of type Dom::DenseMatrix() over the ring of MuPAD® expressions.
For all other matrices A of category Cat::Matrix, the result is returned as a matrix of type Dom::Matrix() over the ring of MuPAD expressions. This includes input matrices A of type Dom::Matrix(...), Dom::SquareMatrix(...), Dom::MatrixGroup(...) etc.
The components of A must not contain symbolic objects which cannot be converted to numerical values via float. Numerical symbolic expressions such as π, , etc. are accepted. They are converted to floats.
Note: 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 round-off effects. It is the user's responsibility to make sure that the diagonalization is feasible and well conditioned.
Symmetric/Hermitian and skew/skew Hermitian 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 π, etc. passed as parameters p1, p2, … are not converted to floats.
Inversion or exponentiation of a matrix may be realized with the functions and 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.
The function is sensitive to the environment variable DIGITS, which determines the numerical working precision.
We compute the matrix power A100:
A := array(1..2, 1..2, [[2, PI], [exp(-10), 0]]): numeric::fMatrix(x -> x^100, A)
Alternatively, you may use the function _power which takes the exponent as a second parameter:
numeric::fMatrix(_power, A, 100)
We compute the square root of a matrix:
A := matrix([[0, 1], [-1, 1]]): B := numeric::fMatrix(sqrt, A)
The small imaginary parts are caused by numerical round-off. We eliminate them by extracting the real parts of the components:
B := map(B, Re)
We verify that B^2 is A. Since A was passed as a matrix of type Dom::Matrix(), the matrix B is also of this type. We may compute the square by the overloaded standard arithmetic using the operator ^:
delete A, B:
We compute with a symbolic parameter t:
A := array(1..2, 1..2, [[0, 1], [-1, 0]]): numeric::fMatrix(exp@_mult, A, t*PI)
We demonstrate the difference between HardwareFloats and SoftwareFloats. The diagonalization of the following matrix is ill-conditioned. The result is dominated by round-off effects:
A := array(1..3, 1..3, [[10, 1, 1 ], [ 0, 1, 1 ], [ 1, 0, 10^(-14)]]): numeric::fMatrix(ln, A, SoftwareFloats)
numeric::fMatrix(ln, A, HardwareFloats)
In the following case, the round-off effects of SoftwareFloats makes the algorithm think that the matrix cannot be diagonalized. Consequently, FAIL is returned. With HardwareFloats, however, a result is computed:
A := array(1..3, 1..3, [[ 1 , 1 , 1 ], [ 0 , 1 , 1 ], [ 10^(-30), 0 , 10^(-30)]]): numeric::fMatrix(ln, A, SoftwareFloats)
numeric::fMatrix(ln, A, HardwareFloats)
A procedure representing a scalar function or , where P is a set of parameters
p1, p2, …
Arbitrary MuPAD objects accepted by f as additional input parameters
Hard, HardwareFloats, Soft, SoftwareFloats
With Hard (or HardwareFloats), computations are done using fast hardware float arithmetic from within a MuPAD session. Hard and HardwareFloats are equivalent. With this option, the input data are converted to hardware floats and processed by compiled C code. The result is reconverted to MuPAD floats and returned to the MuPAD session.
With Soft (or SoftwareFloats) computations are dome using software float arithmetic provided by the MuPAD kernel. Soft and SoftwareFloats are equivalent. SoftwareFloats is used by default if the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY.
Compared to the SoftwareFloats used by the MuPAD kernel, the computation with HardwareFloats may be many times faster. Note, however, that the precision of hardware arithmetic is limited to about 15 digits. Further, the size of floating-point numbers may not be larger than approximately 10308 and not smaller than approximately 10- 308.
If no HardwareFloats or SoftwareFloats are requested explicitly, the following strategy is used: If the current value of DIGITS is smaller than 16 or if the matrix A is a hardware float array of domain type DOM_HFARRAY, then hardware arithmetic is tried. If this is successful, the result is returned.
If the result cannot be computed with hardware floats, software arithmetic by the MuPAD kernel is tried.
If the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY, or if one of the options Soft, SoftwareFloats or Symbolic is specified, MuPAD computes the result with its software arithmetic without trying to use hardware floats first.
There may be several reasons for hardware arithmetic to fail:
If neither HardwareFloats nor SoftwareFloats is specified, the user is not informed whether hardware floats or software floats are used.
If HardwareFloats are specified but fail due to one of the reasons above, a warning is issued that the (much slower) software floating-point arithmetic of the MuPAD kernel is used.
Option, specified as ReturnType = d
Depending on the type of the input matrix A, the matrix f(A) is returned as a matrix of type DOM_ARRAY, DOM_HFARRAY, Dom::Matrix() or Dom::DenseMatrix(). If the algorithm thinks that A is not diagonalizable, then FAIL is returned.
A numerical diagonalization A = X diag(λ1, λ2, …) X- 1 is computed. The columns of X are the (right) eigenvectors of A, the diagonal entries λ1, λ2, … are the corresponding eigenvalues. The function f is mapped to the eigenvalues, the matrix result is computed by
The eigenvector matrix X may be obtained via numeric::eigenvectors(A).
The condition number of the eigenvector matrix is a measure indicating how well conditioned the diagonalization of the matrix A is. If this number is larger than 10DIGITS, then not a single digit of the diagonalization data is trustworthy.
The call numeric::fMatrix(exp, A) corresponds to numeric::expMatrix(A, Diagonalization).