numeric::fMatrix
-- functional
calculus for numerical square matricesnumeric::fMatrix
(f, A, ..)
computes the
matrix f(A,..) with a function f and a square
matrix A.
numeric::fMatrix(f, A, p1, p2, ..)
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 |
The matrix f(A, ..) is returned as an array.
The function is sensitive to the environment variable DIGITS
, which determines the
numerical working precision.
numeric::expMatrix
, numeric::inverse
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.
numeric::fMatrix
produces reliable numerical results for such matrices.p1, p2, ..
may be numerical or symbolic
objects. They must be accepted by f
as 2nd argument, 3rd
argument etc.PI
, sqrt(2)
etc. passed as
parameters p1, p2,..
are not converted to floats.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
.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:
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:
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:
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].
numeric::fMatrix(exp, A)
is equivalent to
numeric::expMatrix(A,
Diagonalization)