Previous Page Next Page Contents

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

Option: vars

Option: Symbolic

Option: ShowAssumptions

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




Do you have questions or comments?


Copyright © SciFace Software GmbH & Co. KG 2000