## `numeric::linsolve` -- solve a system of linear equations

#### Introduction

`numeric::linsolve``(eqs, ..)` solves a system of linear equations.

#### Call(s)

`numeric::linsolve(eqs <, vars> <, Symbolic> <, ShowAssumptions>)`

#### Parameters

`eqs` | - | a list or set of linear equations or expressions

#### Options

``vars`` | - | a list or set of unknowns to solve for. Unknowns may be identifiers or indexed identifiers or expressions.
`Symbolic` | - | prevents conversion of input data to floating point numbers.
`ShowAssumptions` | - | returns information on internal assumptions on symbolic parameters in `eqs`.

#### Returns

Without the option ShowAssumptions a list of simplified equations is returned. It represents the general solution of the system `eqs`. `FAIL` is returned, if the system is not solvable.

With ShowAssumptions a list `[Solution, Constraints, Pivots]` is returned. `Solution` is a list of simplified equations representing the general solution of `eqs`. The lists `Constraints` and `Pivots` contain equations and inequalities involving symbolic parameters in `eqs`. Internally these were assumed to hold true when solving the system. `[FAIL,[],[]]` 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::fsolve`, `numeric::inverse`, `numeric::matlinsolve`, `numeric::polyroots`, `numeric::polysysroots`, `numeric::realroots`, `polylib::realroots`, `solve`

#### Details

- `numeric::linsolve` is a fast numerical linear solver. It is also a recommended solver for linear systems with exact or symbolic coefficients (using Symbolic).
- Expressions are interpreted as homogeneous equations. E.g., the input `[x=y-1,x-y]` is interpreted as the system of equations `[x=y-1,x-y=0]`.

- Without the option Symbolic the input data are converted to floating point numbers. The coefficient matrix A of the system A*x=b represented by `eqs` must not contain non-convertible parameters, unless the option Symbolic is used! If such objects are found, then `numeric::linsolve` automatically switches to its symbolic mode, issuing a warning. Symbolic parameters in the ``right hand side'' b are accepted without warning.

- The numerical working precision is set by the environment variable `DIGITS`.
- The solutions are returned as a list of solved equations of the form

      [x1=value1,x2=value2,..] ,

  where x1,x2,.. are the unknowns. These simplified equations should be regarded as constraints on the unknowns. E.g., if an unknown x[1], say, does not turn up in the form [x[1]=.., ..] in the solution, then there is no constraint on this unknown and it is an arbitrary parameter. Generally, all unknowns that do not turn up on the left hand side of the solved equations are arbitrary parameters spanning the solution space. Cf. example 9.

  In particular, if the empty list is returned as the solution, then there are no constraints whatsoever on the unknowns, i.e., the system is trivial.

- The ordering of the solved equations corresponds to the ordering of the unknowns `vars`. It is recommended that the user specifies `vars` by a a list of unknowns. This guarantees that the solved equations are returned in the expected order. If `vars` are specified by a set, or if no `vars` are specified at all, then an internal ordering is used.
- `numeric::linsolve` returns the general solution of the system `eqs`. It is valid for arbitrary complex values of the symbolic parameters which may be present in `eqs`. If no such solution exists, then `FAIL` is returned. Solutions that are valid only for special values of the symbolic parameters may be obtained with the option ShowAssumptions. Cf. examples 2, 3, 4, 10.
- The solved equations representing the solution are suitable as input for `assign` and `subs`. Cf. example 8.
- `numeric::linsolve` is suitable for solving large sparse systems. Cf. example 6.
- If `eqs` represents a system with a banded coefficient matrix, then this is detected and used by `numeric::linsolve`. Note that in this case it is important to enter `eqs` as a list and to specify the unknowns as a list in order to guarantee the desired form of the coefficient matrix. The elements of sets may be reordered internally, leading to a loss of band structure and, consequently, of efficiency. Cf. example 6.

- `numeric::linsolve` is tuned for speed. For this reason it does not check systematically that the equations `eqs` are indeed linear in the unknowns! For non-linear equations strange things may happen, `numeric::linsolve` might even return wrong results! Cf. example 5.

- `numeric::linsolve` does not react to any assumptions on the unknowns or on symbolic parameters that are set via `assume`.

- Gaussian elimination with partial pivoting is used. Without the option Symbolic, floating point arithmetic is used and the pivoting strategy takes care of numerical stabilization. With Symbolic exact data are assumed and the pivoting strategy tries do maximize speed, not taking care of numerical stabilization! Cf. example 7.

#### Option: `vars`

- If no unknowns are specified by `vars`, then `numeric::linsolve` solves for all symbolic objects in `eqs`. The unknowns are determined internally by `indets(eqs,PolyExpr)`.

#### Option: Symbolic

- This option prevents conversion of the input data to floats.
- This option must be used, if the coefficients of the equations contain symbolic parameters that cannot be converted to floating point numbers.

- This option should not be used for equations with floating point coefficients! Numerical instabilities may occur in floating point operations. Cf. example 7.

#### Option: ShowAssumptions

- This option is only useful, if the equations contain symbolic parameters. Consequently, it should only be used in conjunction with the option Symbolic.

- The format of the return value is changed to `[Solution, Constraints, Pivots]`.

- `Solution` is a set of simplified equations representing the general solution subject to `Constraints` and `Pivots`.
- `Constraints` is a list of equations for symbolic parameters in `eqs`, which are necessary and sufficient to make the system solvable.

  Such constraints arise, if Gaussian elimination of the original equations leads to equations of the form 0=c, where c is some expression involving symbolic parameters in the ``right hand side'' of the system. All such equations are collected in `Constraints`. `numeric::linsolve` assumes that these equations are satisfied and returns a solution.

  If no such constraints arise, then the return value of `Constraints` is the empty list.

- This option changes the return strategy for ``unsolvable'' systems. Without the option Symbolic `FAIL` 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,[],[]]` is returned.

- `Pivots` is a list of inequalities involving symbolic parameters in the coefficient matrix A of the linear system A*x=b represented by `eqs`. Internally, division by pivot elements occurs in the Gaussian elimination. The expressions collected in `Pivots` are the numerators of those pivot elements that contain symbolic parameters. If only numerical pivot elements were used, then the return value of `Pivots` is the empty list.
- Cf. examples 2, 3, 4, 10.

#### Example 1

Equations and variables may be entered as sets or lists:

    >> numeric::linsolve({x = y - 1, x + y = z}, {x, y}),
       numeric::linsolve([x = y - 1, x + y = z], {x, y}),
       numeric::linsolve({x = y - 1, x + y = z}, [x, y]),
       numeric::linsolve([x = y - 1, x + y = z], [x, y])

          [y = 0.5 z + 0.5, x = 0.5 z - 0.5],

             [y = 0.5 z + 0.5, x = 0.5 z - 0.5],

             [x = 0.5 z - 0.5, y = 0.5 z + 0.5],

             [x = 0.5 z - 0.5, y = 0.5 z + 0.5]

With the option Symbolic exact arithmetic is used. The following system has a 1-parameter set of solution, the unknown `x[3]` is arbitrary:

    >> numeric::linsolve([x[1] + x[2] = 2, x[1] - x[2] = 2*x[3]],
                         [x[1], x[2], x[3]], Symbolic)

                        [x[1] = x[3] + 1, x[2] = 1 - x[3]]

The unknowns may be expressions:

    >> numeric::linsolve([f(0) - sin(x + 1) = 2, f(0) = 1 - sin(x + 1)],
                         [f(0), sin(x + 1)])

                          [f(0) = 1.5, sin(x + 1) = -0.5]

The following system does not have a solution:

    >> numeric::linsolve([x + y = 1, x + y = 2], [x, y])

                                       FAIL

#### Example 2

We demonstrate some examples with symbolic coefficients. Note that the option Symbolic has to be used:

    >> eqs := [x + a*y = b, x + A*y = b]:

    >> numeric::linsolve(eqs, [x, y], Symbolic)

                                  [x = b, y = 0]

Note that for a=A this is not the general solution. Using the option ShowAssumptions it turns out, that the above result is the general solution subject to the assumption a<>A:

    >> numeric::linsolve(eqs, [x, y], Symbolic, ShowAssumptions)

                        [[x = b, y = 0], [], [A - a <> 0]]

    >> delete eqs:

#### Example 3

We give a further demonstration of the option ShowAssumptions. The following system does not have a solution for all values of the parameter `a`:

    >> numeric::linsolve([x + y = 1, x + y = a], [x, y], Symbolic)

                                       FAIL

With ShowAssumptions `numeric::linsolve` investigates, under which conditions (on `a`) there is a solution:

    >> numeric::linsolve([x + y = 1, x + y = a], [x, y], Symbolic,
                         ShowAssumptions)

                          [[x = 1 - y], [a - 1 = 0], []]

We conclude that there is a 1-parameter set of solutions for a=1. The constraint in `a` is a linear equation, since the parameter `a` enters the equations linearly. If `a` is regarded as an unknown rather than as a parameter, then the constraint becomes part of the solution:

    >> numeric::linsolve([x + y = 1, x + y = a], [x, y, a], Symbolic,
                         ShowAssumptions)

                           [[x = 1 - y, a = 1], [], []]

#### Example 4

With exact arithmetic `PI` is regarded as a symbolic parameter. The following system has a solution subject to the constraint `PI=1`:

    >> numeric::linsolve([x = x - y + 1, y = PI], [x, y],
                         Symbolic, ShowAssumptions)

                           [[y = PI], [1 - PI = 0], []]

With floating point arithmetic `PI` is converted to `3.1415...`. The system has no solution:

    >> numeric::linsolve([x = x - y + 1, y = PI], [x, y],
                         ShowAssumptions)

                                  [FAIL, [], []]

#### Example 5

Since `numeric::linsolve` does not do a systematic internal check for non-linearities, the user should make sure that the equations to be solved are indeed linear in the unknowns. Otherwise strange things may happen. Garbage is produced for the following non-linear systems:

    >> a := sin(x):

    >> numeric::linsolve([y = 1 - a, x = y], [x, y], Symbolic)

                         [x = 1 - sin(0), y = 1 - sin(0)]

    >> numeric::linsolve([a*x + y = 1, x = y], [x, y], Symbolic)

                    --           1                  1       --
                    |  x = -------------, y = -------------  |
                    |             3                  3       |
                    --     sin(x16 ) + 1      sin(x16 ) + 1 --

Polynomial non-linearities are usually detected. Regarding `x,y,a` as unknowns the following quadratic system yields an error:

    >> delete a: numeric::linsolve([x*a + y = 1, x = y], Symbolic)

          Error: this system does not seem to be linear
          [numeric::linsolve]

This system is linear in `x,y`, if `a` is regarded as a parameter:

    >> numeric::linsolve([x*a + y = 1, x = y], [x, y], Symbolic)

                            --       1          1   --
                            |  x = -----, y = -----  |
                            --     a + 1      a + 1 --

#### Example 6

We solve a large sparse system. The coefficient matrix has only 3 diagonal bands. Note that both the equations as well as the variables are passed as lists. This guarantees that the band structure is not lost internally:

    >> n := 100: x[0] := 0: x[n + 1] := 0:
       eqs := [x[i-1] - 2*x[i] + x[i+1] = 1 $ i = 1..n]:
       vars := [x[i] $ i = 1..n]:
       numeric::linsolve(eqs, vars)

          [x[1] = -50.0, x[2] = -99.0, x[3] = -147.0, x[4] = -194.0,

             ..., x[98] = -147.0, x[99] = -99.0, x[100] = -50.0]

The band structure is lost, if the equations or the unknowns are specified by sets. The following call takes more time than the previous call:

    >> numeric::linsolve({op(eqs)}, {x[i] $ i = 1..n})

          [x[86] = -645.0, x[49] = -1274.0, x[12] = -534.0,

             ... , x[87] = -609.0, x[50] = -1275.0, x[13] = -572.0]

    >> delete n, x:

#### Example 7

The option Symbolic should not be used for equations with floating point coefficients, because the symbolic pivoting strategy favors efficiency instead of numerical stability.

    >> eqs := [x + 10^20*y = 10^20, x + y = 0]:

The float approximation of the exact solution is:

    >> map(numeric::linsolve(eqs, [x, y], Symbolic), map, float)

                                [x = -1.0, y = 1.0]

We now convert the exact coefficients to floating point numbers:

    >> feqs := map(eqs, map, float)

                       [x + 1.0e20 y = 1.0e20, x + y = 0.0]

The default pivoting strategy stabilizes floating point operations. Consequently, one gets a correct result:

    >> numeric::linsolve(feqs, [x, y])

                                [x = -1.0, y = 1.0]

With Symbolic the pivoting strategy optimizes speed, assuming exact arithmetic. Numerical instabilities may occur, if floating point coefficients are involved. The following incorrect result is caused by internal round-off effects (``cancellation''):

    >> numeric::linsolve(feqs, [x, y], Symbolic)

                                 [x = 0, y = 1.0]

    >> delete eqs, feqs:

#### Example 8

We demonstrate that the simplified equations representing the solution can be used for further processing with `subs`:

    >> eqs := [x + y = 1, x + y = a]:

    >> [Solution, Constraints, Pivots] :=
       numeric::linsolve(eqs, [x, y], ShowAssumptions)

                      [[x = 1.0 - 1.0 y], [a - 1.0 = 0], []]

    >> subs(eqs, Solution)

                                [1.0 = 1, 1.0 = a]

The solution can be assigned to the unknowns via `assign`:

    >> assign(Solution): x, y, eqs

                        1.0 - 1.0 y, y, [1.0 = 1, 1.0 = a]

    >> delete eqs, Solution, Constraints, Pivots, x, y:

#### Example 9

If the solution of the linear system is not unique, then some of the unknowns are used as ``free parameters'' spanning the solution space. In the following example the unknown `z` is such a parameter. It does not turn up on the left hand side of the solved equations:

    >> eqs := [x + y = z, x + 2*y = 0, 2*x - z = -3*y, y + z = 0]:

    >> vars := [w, x, y, z]:

    >> Solution := numeric::linsolve(eqs, vars, Symbolic)

                                 [x = 2 z, y = -z]

You may define a function such as the following `NewSolutionList` to rename your free parameters to ``myName1'', ``myName2'' etc. and fill up your list of solved equations accordingly:

    >> NewSolutionList :=
       proc(Solution : DOM_LIST, vars : DOM_LIST, myName : DOM_STRING)
       local i, solvedVars, newEquation;
       begin
         solvedVars := map(Solution, op, 1);
         for i from 1 to nops(vars) do
            if not has(solvedVars, vars[i]) then
               newEquation := vars[i] = genident(myName);
               Solution := listlib::insertAt(
                   subs(Solution, newEquation), newEquation, i)
            end_if
         end_for:
         Solution
       end_proc:

    >> NewSolutionList(Solution, vars, "FreeParameter")

          [w = FreeParameter1, x = 2 FreeParameter2,

             y = -FreeParameter2, z = FreeParameter2]

    >> delete eqs, vars, Solution, NewSolutionList:

#### Example 10

We demonstrate, how a complete solution of the following linear system in `x,y` with symbolic parameters may be found:

    >> eqs := [x + y = A, a*x + b*y = 1, x + c*y = 1]:

    >> numeric::linsolve(eqs, [x, y], Symbolic, ShowAssumptions)

          -- --     A b - 1      1 - A a --
          |  |  x = -------, y = -------  |,
          -- --      b - a        b - a  --

                                                             --
             [b - a - c - A b + A a c + 1 = 0], [b - a <> 0]  |
                                                             --

This is the general solution, assuming a<>b. We now set b=a to investigate further solution branches:

    >> eqs := subs(eqs, b = a):

    >> numeric::linsolve(eqs, [x, y], Symbolic, ShowAssumptions)

          -- --     A c - 1      1 - A --                              --
          |  |  x = -------, y = -----  |, [1 - A a = 0], [c - 1 <> 0]  |
          -- --      c - 1       c - 1 --                              --

This is the general solution for a=b, assuming c<>1. We finally set c=1 to obtain the last solution branch:

    >> eqs := subs(eqs, c = 1):

    >> numeric::linsolve(eqs, [x, y], Symbolic, ShowAssumptions)

                    [[x = A - y], [1 - A = 0, 1 - A a = 0], []]

From the constraints on the symbolic parameters `a` and `A` we conclude that there is a special 1-parameter solution x=1-y for a=b=c=A=1.

    >> delete eqs:

#### Changes

- new function

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