numeric::linsolve
-- solve a
system of linear equationsnumeric::linsolve
(eqs, ..)
solves a system
of linear equations.
numeric::linsolve(eqs <, vars> <, Symbolic> <, ShowAssumptions>)
eqs |
- | a list or set of linear equations or expressions |
|
- | 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 . |
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.
Without the option Symbolic the function is
sensitive to the environment variable DIGITS
, which determines the
numerical working precision.
linalg::matlinsolve
, linsolve
, numeric::fsolve
, numeric::inverse
, numeric::matlinsolve
, numeric::polyroots
, numeric::polysysroots
,
numeric::realroots
,
polylib::realroots
,
solve
numeric::linsolve
is a fast numerical linear solver.
It is also a recommended solver for
linear systems with exact or symbolic coefficients (using Symbolic).[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.
DIGITS
.[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.
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.assign
and
subs
. Cf. example 8.numeric::linsolve
is suitable for solving large sparse
systems. Cf. example 6.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.
vars
vars
, then
numeric::linsolve
solves for all symbolic objects in
eqs
. The unknowns are determined internally by
indets(eqs,PolyExpr)
.This option should not be used for equations with floating point coefficients! Numerical instabilities may occur in floating point operations. Cf. example 7.
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.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
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:
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], [], []]
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, [], []]
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 --
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:
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:
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:
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:
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: