Solve a system of linear equations
This functionality does not run in MATLAB.
numeric::linsolve(eqs
, <vars
>, options
)
numeric::linsolve(eqs, vars)
solves a system
of linear equations eqs
for the unknowns vars
.
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]
.
Note:
Without the option 
The numerical working precision is set by the environment variable DIGITS
.
The solutions are returned as a list of solved equations of the form
,
where x_{1}, x_{2}, … 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; 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, 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.
If no unknowns are specified by vars
, numeric::linsolve
solves
for all symbolic objects in eqs
. The unknowns are
determined internally by indets(eqs, PolyExpr)
.
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, FAIL
is returned. Solutions that are
valid only for special values of the symbolic parameters may be obtained
with the option ShowAssumptions
. See Example 2, Example 3, Example 4, and Example 11.
The solved equations representing the solution are suitable
as input for assign
and subs
. See Example 8.
numeric::linsolve
is suitable for solving
large sparse systems. See 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 specify both the equations
as well as the unknowns by lists to guarantee the desired form of
the coefficient matrix. When using sets, the data may be reordered
internally leading to a loss of band structure and, consequently,
of efficiency. See Example 6.
Note:

Note:

Note:
Gaussian elimination with partial pivoting is used. Without
the option 
Without the option Symbolic
, the function
is sensitive to the environment variable DIGITS
, which determines
the numerical working precision.
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])
With the option Symbolic
, exact arithmetic
is used. The following system has a 1parameter 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)
The unknowns may be expressions:
numeric::linsolve([f(0)  sin(x + 1) = 2, f(0) = 1  sin(x + 1)], [f(0), sin(x + 1)])
The following system does not have a solution:
numeric::linsolve([x + y = 1, x + y = 2], [x, y])
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)
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)
delete eqs:
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)
With ShowAssumptions
, numeric::linsolve
investigates
under which conditions (on the parameter a
) there
is a solution:
numeric::linsolve([x + y = 1, x + y = a], [x, y], Symbolic, ShowAssumptions)
We conclude that there is a 1parameter 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, the constraint becomes part of
the solution:
numeric::linsolve([x + y = 1, x + y = a], [x, y, a], Symbolic, ShowAssumptions)
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)
With floatingpoint arithmetic, PI is
converted to 3.1415...
. The system has no solution:
numeric::linsolve([x = x  y + 1, y = PI], [x, y], ShowAssumptions)
Since numeric::linsolve
does not do a systematic
internal check for nonlinearities, 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 nonlinear
systems:
a := sin(x): numeric::linsolve([y = 1  a, x = y], [x, y], Symbolic)
numeric::linsolve([a*x + y = 1, x = y], [x, y], Symbolic)
Polynomial nonlinearities are usually detected. Regarding x,
y, c
as unknowns, the following quadratic system yields
an error:
numeric::linsolve([x*c + y = 1, x = y], Symbolic)
Error: This system does not seem to be linear. [numeric::linsolve]
Error: This system does not seem to be linear. [numeric::linsolve]
This system is linear in x, y
if c
is
regarded as a parameter:
numeric::linsolve([x*c + y = 1, x = y], [x, y], Symbolic)
delete a:
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 := 500: x[0] := 0: x[n + 1] := 0: eqs := [x[i1]  2*x[i] + x[i+1] = 1 $ i = 1..n]: vars := [x[i] $ i = 1..n]: numeric::linsolve(eqs, vars)
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})
delete n, x, eqs, vars:
The option Symbolic
should not be used for
equations with floatingpoint 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)
We now convert the exact coefficients to floatingpoint numbers:
feqs := map(eqs, map, float)
The default pivoting strategy stabilizes floatingpoint operations. Consequently, one gets a correct result:
numeric::linsolve(feqs, [x, y])
With Symbolic
, the pivoting strategy optimizes
speed, assuming exact arithmetic. Numerical instabilities may occur
if floatingpoint coefficients are involved. The following incorrect
result is caused by internal roundoff effects ("cancellation"):
numeric::linsolve(feqs, [x, y], Symbolic)
delete eqs, feqs:
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)
subs(eqs, Solution)
The solution can be assigned to the unknowns via assign
:
assign(Solution): x, y, eqs
delete eqs, Solution, Constraints, Pivots, x:
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 unknowns z, w
are
such parameters. They do 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 := [x, y, z, w]: Solution := numeric::linsolve(eqs, vars, Symbolic)
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")
delete eqs, vars, Solution, NewSolutionList:
We demonstrate the difference between hardware and software
arithmetic. The following problem is very illconditioned. The results,
both with HardwareFloats
as well as with SoftwareFloats
,
are marred by numerical roundoff:
n:= 10: eqs:= [(_plus(x[j]/(i + j 1) $ j = 1..n) = 1) $ i = 1..n]: vars:= [x[i] $ i = 1..n]: numeric::linsolve(eqs, vars, SoftwareFloats); numeric::linsolve(eqs, vars, HardwareFloats)
This is the exact solution:
numeric::linsolve(eqs, vars, Symbolic);
delete eqs, vars:
We demonstrate how a complete solution of the following linear
system in x
, y
with symbolic
parameters a
, b
, c
, d
may
be found:
eqs := [x + y = d, a*x + b*y = 1, x + c*y = 1]: numeric::linsolve(eqs, [x, y], Symbolic, ShowAssumptions)
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)
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)
From the constraints on the symbolic parameters a
and d
,
we conclude that there is a special 1parameter solution x =
1  y for a = b = c = d =
1.
delete eqs:

A list, set, 

A list or set of unknowns to solve for. Unknowns may be identifiers or indexed identifiers or arithmetical expressions. 

With With Compared to the If no If the result cannot be computed with hardware floats, software arithmetic by the MuPAD kernel is tried. If the current value of There may be several reasons for hardware arithmetic to fail:
If neither If Note that The trailing digits in floatingpoint results computed with
 

Prevents conversion of input data to floatingpoint numbers.
This option overrides This option must be used if the coefficients of the equations contain symbolic parameters that cannot be converted to floatingpoint numbers.
 

Returns information on internal assumptions on symbolic parameters
in This option is only useful if the equations contain symbolic
parameters. Consequently, it should only be used in conjunction with
the option
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 If no such constraints arise, the return value of
See Example 2, Example 3, Example 4, and Example 11.  

Suppresses warnings If symbolic coefficients are found, 
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.