Documentation

numeric::linsolve

Solve a system of linear equations

MuPAD® notebooks will be removed in a future release. Use MATLAB® live scripts instead.

MATLAB live scripts support most MuPAD functionality, though there are some differences. For more information, see Convert MuPAD Notebooks to MATLAB Live Scripts.

Syntax

numeric::linsolve(eqs, <vars>, options)

Description

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 Symbolic, the input data are converted to floating-point numbers. The coefficient matrix A of the system Ax = 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. This warning may be suppressed via NoWarning. 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 ,

where x1, x2, … are the unknowns. These simplified equations should be regarded as constraints on the unknowns. E.g., if an unknown x1, say, does not turn up in the form [x1 = …, …] 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

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! See Example 5.

Note

numeric::linsolve does not react to any properties of the unknowns or of symbolic parameters that are set via assume.

Note

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! See Example 7.

Environment Interactions

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

Examples

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])    With the option Symbolic, exact arithmetic is used. The following system has a 1-parameter set of solution; the unknown x3 is arbitrary:

numeric::linsolve([x + x = 2, x - x = 2*x],
[x, x, x], 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]) 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) 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 aA:

numeric::linsolve(eqs, [x, y], Symbolic, ShowAssumptions) 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) 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 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, the constraint becomes part of the solution:

numeric::linsolve([x + y = 1, x + y = a], [x, y, a], Symbolic,
ShowAssumptions) 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) 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) 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) numeric::linsolve([a*x + y = 1, x = y], [x, y], Symbolic) Polynomial non-linearities 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:

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 := 500: x := 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) 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:

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) We now convert the exact coefficients to floating-point numbers:

feqs := map(eqs, map, float) The default pivoting strategy stabilizes floating-point 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 floating-point coefficients are involved. The following incorrect result is caused by internal round-off effects (“cancellation”):

numeric::linsolve(feqs, [x, y], Symbolic) 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) subs(eqs, Solution) The solution can be assigned to the unknowns via assign:

assign(Solution):
x, y, eqs delete eqs, Solution, Constraints, Pivots, x:

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 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:

Example 10

We demonstrate the difference between hardware and software arithmetic. The following problem is very ill-conditioned. The results, both with HardwareFloats as well as with SoftwareFloats, are marred by numerical round-off:

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:

Example 11

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 ab. 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 1-parameter solution x = 1 - y for a = b = c = d = 1.

delete eqs:

Parameters

 eqs A list, set, array, or matrix (Cat::Matrix) of linear equations or arithmetical expressions vars A list or set of unknowns to solve for. Unknowns may be identifiers or indexed identifiers or arithmetical expressions.

Options

Hard, HardwareFloats, Soft, SoftwareFloats

With Hard (or HardwareFloats), computations are done using fast hardware float arithmetic from within a MuPAD® session. Hard and HardwareFloats are equivalent. With this option, the input data are converted to hardware floats and processed by compiled C code. The result is reconverted to MuPAD floats and returned to the MuPAD session.

With Soft (or SoftwareFloats) computations are dome using software float arithmetic provided by the MuPAD kernel. Soft and SoftwareFloats are equivalent. SoftwareFloats is used by default if the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY.

Compared to the SoftwareFloats used by the MuPAD kernel, the computation with HardwareFloats may be many times faster. Note, however, that the precision of hardware arithmetic is limited to about 15 digits. Further, the size of floating-point numbers may not be larger than approximately 10308 and not smaller than approximately 10- 308.

If no HardwareFloats or SoftwareFloats are requested explicitly, the following strategy is used: If the current value of DIGITS is smaller than 16 or if the matrix A is a hardware float array of domain type DOM_HFARRAY, then hardware arithmetic is tried. If this is successful, the result is returned.

If the result cannot be computed with hardware floats, software arithmetic by the MuPAD kernel is tried.

If the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY, or if one of the options Soft, SoftwareFloats or Symbolic is specified, MuPAD computes the result with its software arithmetic without trying to use hardware floats first.

There may be several reasons for hardware arithmetic to fail:

• The current value of DIGITS is larger than 15.

• The data contains symbolic objects.

• The data contains numbers larger than 10308 or smaller than 10- 308 that cannot be represented by hardware floats.

If neither HardwareFloats nor SoftwareFloats is specified, the user is not informed whether hardware floats or software floats are used.

If HardwareFloats are specified but fail due to one of the reasons above, a warning is issued that the (much slower) software floating-point arithmetic of the MuPAD kernel is used.

Note that HardwareFloats can only be used if all input data can be converted to floating-point numbers.

The trailing digits in floating-point results computed with HardwareFloats and SoftwareFloats may differ.

Note

For ill-conditioned systems, the result is subject to round-off errors. The results returned with HardwareFloats and SoftwareFloats may differ significantly! SeeExample 10.

Symbolic

Prevents conversion of input data to floating-point numbers. This option overrides HardwareFloats and SoftwareFloats.

This option must be used if the coefficients of the equations contain symbolic parameters that cannot be converted to floating-point numbers.

Note

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

ShowAssumptions

Returns information on internal assumptions on symbolic parameters in eqs.

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

Note

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, the return value of Constraints is the empty list.

Pivots is a list of inequalities involving symbolic parameters in the coefficient matrix A of the linear system Ax = b represented by eqs. Internally, division by pivot elements occurs in the Gaussian elimination. The expressions collected in Pivots are the numerators of the pivot elements that contain symbolic parameters. If only numerical pivot elements were used, the return value of Pivots is the empty list.

Note

The option ShowAssumptions 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.

See  Example 2, Example 3, Example 4, and Example 11.

NoWarning

Suppresses warnings

If symbolic coefficients are found, numeric::linsolve automatically switches to the Symbolic mode with a warning. With this option, this warning is suppressed; numeric::linsolve still uses the symbolic mode for symbolic coefficients, i.e., exact arithmetic without floating-point conversions is used.

Return Values

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.

Mathematical Modeling with Symbolic Math Toolbox

Get examples and videos