Previous Page Next Page Contents

numeric::matlinsolve -- solve a linear matrix equation

Introduction

numeric::matlinsolve(A, B, ..) returns the matrix solution X of the matrix equation A*X=B.

Call(s)

numeric::matlinsolve(A, B <, Symbolic> <, ShowAssumptions>)

Parameters

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,..].

Options

Symbolic - prevents conversion of input data to floating point numbers
ShowAssumptions - returns information on internal assumptions on symbolic parameters in A and B

Returns

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.

Side Effects

Without the option Symbolic the function is sensitive to the environment variable DIGITS, which determines the numerical working precision.

Related Functions

linalg::matlinsolve, linsolve, numeric::inverse, numeric::linsolve, solve

Details

Option: Symbolic

Option: ShowAssumptions

Example 1

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:

Example 2

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:

Example 3

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:

Example 4

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:

Example 5

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:

Example 6

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:

Example 7

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:

Example 8

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:

Changes




Do you have questions or comments?


Copyright © SciFace Software GmbH & Co. KG 2000