numeric::matlinsolve
--
solve a linear matrix equationnumeric::matlinsolve
(A, B, ..)
returns the
matrix solution X of the matrix equation
A*X=B.
numeric::matlinsolve(A, B <, Symbolic> <, ShowAssumptions>)
A |
- | an m x n matrix of domain type DOM_ARRAY or of category Cat::Matrix . |
B |
- | an m x p matrix of domain type DOM_ARRAY or of category Cat::Matrix . Column vectors
B may also be represented by a 1-dimensional
array(1..m,[B1,B2,..]) or by a list
[B1,B2,..] . |
Symbolic |
- | prevents conversion of input data to floating point numbers |
ShowAssumptions |
- | returns information on internal assumptions on
symbolic parameters in A and B |
Without the option ShowAssumptions a list
[X,KernelBasis]
is returned. The solution X
is an n x p array. KernelBasis
is an n x
d array. Its columns span the kernel of A.
[FAIL,NIL]
is returned, if the system is not solvable.
With ShowAssumptions a list
[X,KernelBasis,Constraints,Pivots]
is returned. The lists
Constraints
and Pivots
contain equations and
inequalities involving symbolic parameters in A
and
B
. Internally these were assumed to hold true when solving
the system. [FAIL,NIL,[],[]]
is returned, if the system is
not solvable.
Without the option Symbolic the function is
sensitive to the environment variable DIGITS
, which determines the
numerical working precision.
linalg::matlinsolve
, linsolve
, numeric::inverse
, numeric::linsolve
, solve
numeric::matlinsolve
is a fast numerical linear
solver. It is also a recommended solver for linear systems with exact or
symbolic coefficients (use option Symbolic).Without Symbolic the input
data are converted to floating point numbers. The matrix A
must not contain non-convertible parameters, unless Symbolic is used! If such objects are found, then
numeric::matlinsolve
automatically switches to its
symbolic mode, issuing a warning. Symbolic parameters in B
are accepted without warning.
DIGITS
.X
is a special solution of
A*X=B.KernelBasis
is the most
general solution of A*X=0. Its columns span the
d-dimensional kernel of A.KernelBasis
is the null vector represented by
array(1..n,1..1,[0,..,0])
.KernelBasis=array(1..n,1..1,[0,..,0])
, then the
dimension d of the kernel of A is
d:=0
. Otherwise it is
d:=op(KernelBasis,[0,3,2]])
.Due to roundoff errors some or all basis vectors in
the kernel of A
may be missed in the numeric mode.
X
in conjunction with
KernelBasis
provides the general solution of
A*X=B. Solutions generated without the option ShowAssumptions are valid for arbitrary complex values of
the symbolic parameters which may be present in A
and
B
. If no such solution exists, then
[FAIL,NIL]
is returned. Solutions that are valid only for
special values of the symbolic parameters may be obtained with ShowAssumptions. Cf. examples 3, 4, 7.numeric::matlinsolve
internally uses a sparse
representation of the matrices. It is suitable for solving large sparse
systems. Cf. example 5. Note, however, that the
input requires specification of entire matrices. For this reason you
should consider using numeric::linsolve
which allows to
input the linear system as a sparse list of symbolic equations.numeric::matlinsolve
does not react to
any assumptions on symbolic parameters in A,B
that are set
via assume
.
Gaussian elimination with partial pivoting is used. Without option Symbolic the pivoting strategy takes care of numerical stabilization. With Symbolic exact data are assumed. The symbolic pivoting strategy tries do maximize speed and does not take care of numerical stabilization! Cf. example 6.
Cat::Matrix
objects A
from
matrix domains such as Dom::Matrix
(..)
or Dom::SquareMatrix
(..)
are internally converted to arrays over expressions via
A::dom::expr(A)
. Note that linalg::matlinsolve
must be
used, when the solution is to be computed over the component domain.
Cf. example 8. Note that the option Symbolic should be used, if the entries cannot be
converted to numerical expressions.
This option should not be used for matrices with floating point entries! Numerical instabilities may occur in floating point operations. Cf. example 6.
This option changes the format of the return value
to [X, KernelBasis, Constraints, Pivots]
.
X
and KernelBasis
represent the general
solution subject to Constraints
and
Pivots
.Constraints
is a list of equations for symbolic
parameters in B which are necessary and sufficient for
AX=B to be solvable.
Such constraints arise, if Gaussian elimination leads to equations
of the form 0=c, where c is some expression
involving symbolic parameters contained in B. All such
equations are collected in Constraints
.
numeric::matlinsolve
assumes that these equations are
satisfied and returns a special solution X
.
If no such constraints arise, then the return value of
Constraints
is the empty list.
Constraints
usually is a
non-linear list of equation for the symbolic parameters. It is
not investigated by numeric::matlinsolve
, i.e., solutions
may be returned, even if the Constraints
cannot be
satisfied. Cf. example 3.
This option changes the return strategy for
``unsolvable'' systems. Without the option ShowAssumptions the result [FAIL,NIL]
is
returned, whenever Gaussian elimination produces an equation
0=c with non-zero c. With ShowAssumptions such equations are returned via
Constraints
, provided c involves symbolic
parameters. If c is a purely numerical value, then
[FAIL,NIL,[],[]]
is returned.
Pivots
is a list of inequalities involving symbolic
parameters in A. Internally, division by pivot elements
occurs in the Gaussian elimination. The expressions collected in
Pivots
are the numerators of those pivot elements that
involve symbolic parameters contained in A. If only
numerical pivot elements are used, then the return value of
Pivots
is the empty list.We use equivalent input formats
B1
,...,B4
to represent a vector with
components [a,PI]. First, this vector is defined as a
2-dimensional array:
>> A := array(1..2, 1..3, [[1, 2, 3],[1, 1, 2]]):
>> B1 := array(1..2, 1..1, [[a], [PI]]):
>> numeric::matlinsolve(A, B1)
-- +- -+ +- -+ -- | | 6.283185307 - 1.0 a | | -1.0 | | | | | | | | | | 1.0 a - 3.141592654 |, | -1.0 | | | | | | | | | | 0 | | 1 | | -- +- -+ +- -+ --
Next, we use a 1-dimensional array and compute an exact solution:
>> B2 := array(1..2, [a, PI]):
>> numeric::matlinsolve(A, B2, Symbolic)
-- +- -+ +- -+ -- | | 2 PI - a | | -1 | | | | | | | | | | a - PI |, | -1 | | | | | | | | | | 0 | | 1 | | -- +- -+ +- -+ --
Now a list is used to specify the vector. No internal
assumptions were used by numeric::matlinsolve
to obtain
the solution:
>> B3 := [a, PI]:
>> numeric::matlinsolve(A, B3, ShowAssumptions)
-- +- -+ +- -+ -- | | 6.283185307 - 1.0 a | | -1.0 | | | | | | | | | | 1.0 a - 3.141592654 |, | -1.0 |, [], [] | | | | | | | | | 0 | | 1 | | -- +- -+ +- -+ --
Finally, we use Dom::Matrix
objects to specify the
system. Note that the result are still arrays:
>> A := Dom::Matrix()([[1, 2, 3],[1, 1, 2]]):
>> B4 := Dom::Matrix()([a, PI]):
>> numeric::matlinsolve(A, B4)
-- +- -+ +- -+ -- | | 6.283185307 - 1.0 a | | -1.0 | | | | | | | | | | 1.0 a - 3.141592654 |, | -1.0 | | | | | | | | | | 0 | | 1 | | -- +- -+ +- -+ --
>> delete A, B1, B2, B3, B4:
We invert a matrix by solving A*X=1:
>> A := array(1..3, 1..3, [[1, 1, 0], [0, 1, 1], [0, 0, 1]]):
>> B := array(1..3, 1..3, [[1, 0, 0], [0, 1, 0], [0, 0, 1]]):
>> InverseOfA := numeric::matlinsolve(A, B, Symbolic)[1]
+- -+ | 1, -1, 1 | | | | 0, 1, -1 | | | | 0, 0, 1 | +- -+
>> delete A, B, InverseOfA:
We solve an equation with a symbolic parameter x:
>> M := Dom::Matrix():
>> A := M(3, 3, [[2, 2, 3], [1, 1, 2], [3, 3, 5]]):
>> B := M(3, 1, [sin(x)^2, cos(x)^2, 0]):
>> [X, Kernel, Constraints, Pivots] := numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)
-- +- -+ -- | | 2 | +- -+ | | | 5 sin(x) | | -1 | | | | | | | 2 2 | | | 0 |, | 1 |, [cos(x) + sin(x) = 0], [] | | | | | | | | | 2 | | 0 | | | | - 3 sin(x) | +- -+ | -- +- -+ --
This solution holds subject to the constraint
sin(x)^2+cos(x)^2=0 on the parameter x.
numeric::matlinsolve
does not investigate the
Constraints
and does not realize that they cannot be
satisfied. We check the consistency of the ``result'' by inserting the
solution into the original system. For convenience we convert the
arrays X
and Kernel
to Dom::Matrix
objects and use the
overloaded operators *
and -
for matrix
multiplication and subtraction:
>> A*M(X) - B, A*M(Kernel)
+- -+ | 0 | +- -+ | | | 0 | | 2 2 | | | | - cos(x) - sin(x) |, | 0 | | | | | | 0 | | 0 | +- -+ +- -+
>> delete M, A, B, X, Kernel, Constraints, Pivots:
We give a further demonstration of the option ShowAssumptions. The following system does not have a
solution for all values of the parameter a
:
>> A := array(1..2, 1..2, [[1, 1], [1, 1]]):
>> B := array(1..2, 1..1, [[1], [a]]):
>> numeric::matlinsolve(A, B, Symbolic)
[FAIL, NIL]
With ShowAssumptions
numeric::matlinsolve
investigates, under which conditions
(on a
) there is a solution:
>> numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)
-- +- -+ +- -+ -- | | 1 | | -1 | | | | |, | |, [a - 1 = 0], [] | | | 0 | | 1 | | -- +- -+ +- -+ --
We conclude that there is a 1-dimensional solution space for a=1.
>> delete A, B:
We solve a large sparse system with 3 diagonal bands:
>> n := 100: A := Dom::Matrix()(n, n, [1, -2, 1], Banded):
>> B := array(1..n, [1 $ n]): numeric::matlinsolve(A, B)
-- +- -+ +- -+ -- | | -50.0 | | 0 | | | | | | | | | | -99.0 | | 0 | | | | | | | | ... ... | | | | | | | | -50.0 | | 0 | | -- +- -+ +- -+ --
>> delete n, A, B:
The option Symbolic should not be used for equations with floating point coefficients, because the symbolic pivoting strategy favors efficiency instead of numerical stability.
>> A := array(1..2, 1..2, [[1, 10^20], [1, 1]]):
>> B := array(1..2, 1..1, [[10^20], [0]]):
The float approximation of the exact solution is:
>> map(numeric::matlinsolve(A, B, Symbolic)[1], float)
+- -+ | -1.0 | | | | 1.0 | +- -+
We now convert the exact input data to floating point approximations:
>> A := map(A, float): B := map(B, float):
The default pivoting strategy stabilizes floating point operations. Consequently, one gets a correct result:
>> numeric::matlinsolve(A, B)[1]
+- -+ | -1.0 | | | | 1.0 | +- -+
With the option Symbolic the pivoting strategy optimizes speed, assuming exact arithmetic. Numerical instabilities may occur, if floating point coefficients are involved. The following result is caused by internal round-off effects (``cancellation''):
>> numeric::matlinsolve(A, B, Symbolic)[1]
+- -+ | 0 | | | | 1.0 | +- -+
>> delete A, B:
We demonstrate, how a complete solution of the following linear system with symbolic parameters may be found:
>> A := array(1..3, 1..2, [[1, 1], [a, b], [1, c]]):
>> B := array(1..3, 1..1, [[1], [1], [1]]):
>> numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)
-- +- -+ -- | | b - 1 | | | | ----- | +- -+ | | | b - a | | 0 | | | | |, | |, [a c - c - a + 1 = 0], [b - a <> 0] | | | 1 - a | | 0 | | | | ----- | +- -+ | | | b - a | | -- +- -+ --
This is the general solution, assuming a<>b. We now set b=a to investigate further solution branches:
>> A := subs(A, b = a):
>> numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)
-- +- -+ +- -+ -- | | 1 | | 0 | | | | |, | |, [1 - a = 0], [c - 1 <> 0] | | | 0 | | 0 | | -- +- -+ +- -+ --
This is the general solution for a=b, assuming c<>1. We finally set c=1 to obtain the last solution branch:
>> A := subs(A, c = 1):
>> numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)
-- +- -+ +- -+ -- | | 1 | | -1 | | | | |, | |, [1 - a = 0], [] | | | 0 | | 1 | | -- +- -+ +- -+ --
From the constraint on a
we conclude that
there is a 1-dimensional solution space for the special values
a=b=c=1 of the symbolic parameters.
>> delete A, B:
Matrices from a domain such as Dom::Matrix
(..)
are
converted to arrays with numbers or expressions. Hence
numeric::matlinsolve
finds no solution for the following
system:
>> M := Dom::Matrix(Dom::IntegerMod(7)):
>> A := M([[1, 4], [6, 3], [3, 2]]): B := M([[9], [5], [0]]):
>> numeric::matlinsolve(A, B)
[FAIL, NIL]
Use linalg::matlinsolve
to solve the
system over the coefficient field of the matrices. A solution does
exist over the field Dom::IntegerMod
(7)
:
>> linalg::matlinsolve(A, B)
+- -+ | 1 mod 7 | | | | 2 mod 7 | +- -+
>> delete M, A, B: