http://qs1969.pair.com?node_id=347877

This is a little howto for using matrices to solve simultaneous linear equations. Don't let it be confusing, you can just plug your values into the examples to get answers. I have included a bit of theory, so that you are not "confused" by unexpected results.

When you solve simultaneous equations, it is easiest to use matrix equations.

Generally the solution goes like this: Say you have a matrix equation

Ax=b

where A is a 'n x n' matrix, x is a 'n x 1' matrix (vector), and b is a set of known values of 'n x 1' then the solution is to find the inverse of the matrix A and mutiply it by b

x = inv(A) * b

This all corresponds with solving simultaneous linear equations. You need n number of equations to solve for n variables.

There is only one glitch.... there is no inverse to a matrix, if there is not a unique solution. And there is no inverse, if the determinant=0. So checking the determinant is an easy way of telling whether you will get a unique solution. Just because you don't get a unique solution, dosn't mean there are not alot of good solutions. In Math::Matrix::Real, you can use the $dim variable to see what form your solutions is. If $dim=0, the solution is a unique point. If $dim=1, the solution is a line. If $dim=2, the solution is a plane. etc...up into mind-boggling spaces.

The matrix inverse method is "mathematically pure", but not always easy to implement in real solutions. Most matrix solutions use a method called LR-decomposition. It is well explained in perldoc Math::Matrix::Real.

So here are 2 examples, where I solve the same set of equations, with Math::Matrix, and Math::MatrixReal. The first is not unique, where the answer is a plane, and you can see where some confusion can be generated, when Math::Matrix gives different answers from Math::MatrixReal.So always check the determinant or the $dim variable, so see if the solution is unique. (I've mixed my matrix loading notation to show you the equivalence of the different available methods).

EXAMPLE 1

####################################################################

With Math::Matrix 

#!/usr/bin/perl
use warnings;
use Math::Matrix;

#say you have equations
# 1x   2y   3z   4w  = 17
# 5x   6y   7z   8w  = 18
# 9x   10y  11z  12w = 19
# 13x  14y  15z  16w = 20

$a = new Math::Matrix ([1,2,3,4],
                       [5,6,7,8],
                       [9,10,11,12],
                       [13,14,15,16]);

$x = new Math::Matrix ([17,18,19,20]);

print "det-> ",$a->determinant,"\n";

$a->print("A\n");
$E = $a->concat($x->transpose);
$E->print("Equation system\n");
$s = $E->solve;
$s->print("Solutions s\n");
$a->multiply($s)->print("A*s\n");

[download]
OUTPUT
det-> 0
A
1.00000    2.00000    3.00000    4.00000
5.00000    6.00000    7.00000    8.00000
9.00000   10.00000   11.00000   12.00000
13.00000  14.00000   15.00000   16.00000

Equation system
1.00000    2.00000    3.00000    4.00000   17.00000
5.00000    6.00000    7.00000    8.00000   18.00000
9.00000   10.00000   11.00000   12.00000   19.00000
13.00000  14.00000   15.00000   16.00000   20.00000

Solutions s
-16.50000
16.75000
0.00000
0.00000

A*s
17.00000
18.00000
19.00000
20.00000

[download]
###########################################################

with Math::MatrixReal 

#!/usr/bin/perl
use warnings;
use strict;
use Math::MatrixReal;

my ($a1, $b1, $c1, $d1, $e1,
$a2, $b2, $c2, $d2, $e2,
$a3, $b3, $c3, $d3, $e3,
$a4, $b4, $c4, $d4, $e4) =
(1,  2,  3,  4,  17,
 5,  6,  7,  8,  18,
 9,  10, 11, 12, 19,
 13, 14, 15, 16, 20);

my $A = Math::MatrixReal->new_from_string(<<EOF);
[ $a1 $b1 $c1 $d1 ]
[ $a2 $b2 $c2 $d2 ]
[ $a3 $b3 $c3 $d3 ]
[ $a4 $b4 $c4 $d4 ]
EOF

my $b = Math::MatrixReal->new_from_string(<<EOF);
[ $e1 ]
[ $e2 ]
[ $e3 ]
[ $e4 ]
EOF

print "det-> ",$A->det(),"\n";

my $LR = $A->decompose_LR();

my ($dim, $solution, $B) = $LR->solve_LR($b);

#print "dim->$dim\n\nsolution->\n$solution\n\nbase->\n$B\n\n";
print "dim->$dim\n\nsolution->\n$solution\n\nbase->\n$B\n\n";

my $out = $A * $solution;
print $out,"\n";

[download]
OUTPUT:
det-> 0
dim->2

solution->
[ -5.333333333333E+00 ]
[  0.000000000000E+00 ]
[  0.000000000000E+00 ]
[  5.583333333333E+00 ]


[  1.700000000000E+01 ]
[  1.800000000000E+01 ]
[  1.900000000000E+01 ]
[  2.000000000000E+01 ]

[download]
############################################################

So you can see some confusion can be generated here. 2 different answers are generated, but they both are good, since the solutions are not unique, as evidenced by $dim = 2, and determinant=0. This means that the answer is an entire set of points which make a plane.

##############################################################

So here is another example(similar) except the solution is unique.

EXAMPLE 2

########################################################

With Math::Matrix 

#!/usr/bin/perl
use warnings;
use Math::Matrix;

$a = new Math::Matrix ([1, -2,  3,  4],
                       [4, -5,  6,  7],
                       [7,  8, -9,  10],
                       [8,  9, -9, -11],
                       );

$x = new Math::Matrix ([17,18,19,20]);

print "det-> ",$a->determinant,"\n";

$a->print("A\n");
$E = $a->concat($x->transpose);
$E->print("Equation system\n");
$s = $E->solve;
$s->print("Solutions s\n");
$a->multiply($s)->print("A*s\n");

[download]
OUTPUT:
det-> -1104

A
1.00000   -2.00000    3.00000    4.00000
4.00000   -5.00000    6.00000    7.00000
7.00000    8.00000   -9.00000   10.00000
8.00000    9.00000   -9.00000  -11.00000

Equation system
1.00000   -2.00000    3.00000    4.00000   17.00000
4.00000   -5.00000    6.00000    7.00000   18.00000
7.00000    8.00000   -9.00000   10.00000   19.00000
8.00000    9.00000   -9.00000  -11.00000   20.00000

Solutions s
1.45652
18.03261
16.02899
0.88043

A*s
17.00000
18.00000
19.00000
20.00000

[download]
#######################################################

With Math::MatrixReal: 

#!/usr/bin/perl
use warnings;
use strict;
use Math::MatrixReal;

my ($a1, $b1, $c1, $d1, $e1,
$a2, $b2, $c2, $d2, $e2,
$a3, $b3, $c3, $d3, $e3,
$a4, $b4, $c4, $d4, $e4) =
(1,  -2,   3,   4,  17,
 4,  -5,   6,   7,  18,
 7,   8,  -9,  10,  19,
 8,   9,  -9, -11,  20);


my $A = Math::MatrixReal->new_from_string(<<EOF);
[ $a1 $b1 $c1 $d1 ]
[ $a2 $b2 $c2 $d2 ]
[ $a3 $b3 $c3 $d3 ]
[ $a4 $b4 $c4 $d4 ]
EOF

my $b = Math::MatrixReal->new_from_string(<<EOF);
[ $e1 ]
[ $e2 ]
[ $e3 ]
[ $e4 ]
EOF

print "det-> ",$A->det(),"\n";

my $LR = $A->decompose_LR();

my ($dim, $solution, $B) = $LR->solve_LR($b);

#print "dim->$dim\n\nsolution->\n$solution\n\nbase->\n$B\n\n";
print "dim->$dim\n\nsolution->\n\n$solution\n\n";

my $out = $A * $solution;
print $out,"\n";

[download]
OUTPUT: 
det-> -1104
dim->0

solution->

[  1.456521739130E+00 ]
[  1.803260869565E+01 ]
[  1.602898550725E+01 ]
[  8.804347826087E-01 ]


[  1.700000000000E+01 ]
[  1.800000000000E+01 ]
[  1.900000000000E+01 ]
[  2.000000000000E+01 ]

[download]
########################################################

Ok, now the solution is unique, and Math::Matrix gives the exact same solution as Math::MatrixReal

###########################################################

Finally I just want to demonstrate the validity of the original mathematical method of solving by multiplying both sides of the matrix equation by the matrix inverse. It is basic algebra, in that a matrix times it's inverse equals 1, or in matrix terms, a matrix of all zeros except the diagonal being 1's. Like this: 

[1  0  0  0],
[0  1  0  0],
[0  0  1  0],
[0  0  0  1]

[download]
Also we can use the matrix from the second example, since it's determinant is non-zero, and therefore it has an inverse. 
#!/usr/bin/perl
use warnings;
use Math::MatrixReal;

$a = Math::MatrixReal->new_from_rows ([
[1, -2,  3,   4],
[4, -5,  6,   7],
[7,  8, -9,  10],
[8,  9, -9, -11],
]);

print $a,"\n";

my $b = Math::MatrixReal->new_from_cols ([
[17,18,19,20]
]);
print $b,"\n";

print $a->det(),"\n";

$inva = $a->inverse();
print $inva,"\n";

$c = $inva * $b ;

print "solution->\n$c\n";

[download]
OUTPUT:
[  1.000000000000E+00 -2.000000000000E+00  3.000000000000E+00  4.00000
+0000000E+00 ]
[  4.000000000000E+00 -5.000000000000E+00  6.000000000000E+00  7.00000
+0000000E+00 ]
[  7.000000000000E+00  8.000000000000E+00 -9.000000000000E+00  1.00000
+0000000E+01 ]
[  8.000000000000E+00  9.000000000000E+00 -9.000000000000E+00 -1.10000
+0000000E+01 ]

[  1.700000000000E+01 ]
[  1.800000000000E+01 ]
[  1.900000000000E+01 ]
[  2.000000000000E+01 ]

-1104
[ -1.684782608696E-01  1.739130434783E-01  5.434782608696E-03  5.43478
+2608696E-02 ]
[  1.595108695652E+00 -6.304347826087E-01 -3.532608695652E-02  1.46739
+1304348E-01 ]
[  1.362318840580E+00 -4.492753623188E-01 -8.695652173913E-02  1.30434
+7826087E-01 ]
[  6.793478260870E-02 -2.173913043478E-02  4.619565217391E-02 -3.80434
+7826087E-02 ]

solution->
[  1.456521739130E+00 ]
[  1.803260869565E+01 ]
[  1.602898550725E+01 ]
[  8.804347826087E-01 ]

[download]
#################################################################################


So you can see, that the answers agree, but we can only use the inverse method when the determinant does not equal 0. Corrections, comments, and suggestions welcome. (I'm doing this from memory from my matrix theory 25 years ago :-) )
I'm not really a human, but I play one on earth. flash japh

janitored by ybiC: Balanced <readmore> tags around longish 'examples' section

Replies are listed 'Best First'.
Re: Solving Simultaneous Equations with Matrices
by thor (Priest) on Apr 24, 2004 at 20:51 UTC
    The thing of it is that in general, it is not a good idea to calculate the inverse of a matrix. IIRC, that isn an O(n**3) operation. Most times, one solves a system of linear equations by first reducing the n x n matrix to Row-Reduced Echelon Form through Gaussian Elimination. If the matrix has full rank (a fancy way of saying that the system is uniquely solvable), then it is trivial to solve the system.

    thor

[reply]
      The method that he suggested for inverting a matrix, using the determinant as described, takes n! time. Which is considerably worse than O(n**3), and furthermore with that many terms in the calculation, round-off errors can get bad, quickly.

      However inverting a matrix has a O(n**3) algorithm. What you do is write [M|I] where M is the n*n matrix that you want to invert and I is the n*n identity matrix. Then through row-reduction try to turn that into [I|M']. You will succeed iff M is invertable, and M' will be the inverse of M.

      To understand why this works, applying a row operation is the same as left-multiplying by a matrix. So the series of row operations that you did in the reduction is a series of matrix multiplications that turn M into I. In other words you multiplied by the inverse of M. Luckily you also recorded what happens if you start with the identity matrix and do the same sequence of multiplications. But the inverse of M times I is going to be the inverse of M.

      However one would only want to go through the work of calculating an inverse if you had to solve the same system of equations for many different values.

      The determinant formula for an inverse is called Cramer's rule. It has tremendous theoretical value, it is faster to calculate for 2x2 matrices, and it is about even with row reduction for 3x3 matrices. When you get to 4x4, though, row reduction wins significantly, and never looks back.

[reply]
[d/l]
[select]
Re: Solving Simultaneous Equations with Matrices
by gjb (Vicar) on Apr 25, 2004 at 07:41 UTC

    For a bit of background on this problem and the problems that can arise when solving sets of equations, check out Numerical recipes for C (or any of the other languages such as Fortran77, Fortran90 and C++ although the latter is not available online. The book is quite nice to read, actually.

    A nice book that's not available online is Real computing made real by Forman S. Acton.

    Just my 2 cents, -gjb-

[reply]
Re: Solving Simultaneous Equations with Matrices
by flyingmoose (Priest) on Apr 26, 2004 at 16:01 UTC
    obligatory reference to arbitrary alternative wheels: PDL.
[reply]
      You know, I looked at PDL when I wrote this, and it was nearly impossible to use compared to Math::Matrix or Math::Matrix::Real. PDL has a "simq" method to deal with simultaneous equations, but I couldn't figure out from the docs how to use it. The PDL module Math::Matrix:Ops has about 5 lines describing simq......the interface is very obscure.

      I also tried to use PDL to get the inverse of a matrix, which had no inverse, and was given results. Either I had the rows - columns setup wrong, or there is something wacky with PDL.

      I'm not really a human, but I play one on earth. flash japh
[reply]
        For future reference: a slight surprise with PDL is the first dimension is column, second is row, unlike with Fortran, normal maths, etc. So sequence(2,3) has 2 columns, not 3. Printing out the ndarray shows what's going on. PDL::LinearAlgebra has many useful matrix functions that wrap lightning-fast LAPACK routines. I'd recommend that over the routines in PDL::MatrixOps, but see https://metacpan.org/pod/PDL::MatrixOps#lu_backsub for a comparison of the various available ways to solve linear problems with PDL.
[reply]
[d/l]
Re: Solving Simultaneous Equations with Matrices
by tsee (Curate) on Apr 28, 2004 at 15:08 UTC

    Another alternative would be to use Math::Symbolic to solve certain classes of equations. The manual of Math::Symbolic::MiscAlgebra has the following on it: 

      linear_solve
     Calculates the solutions x (vector) of a linear
     equation system of the form "Ax = b" with "A"
     being a matrix, "b" a vector and the solution "x"
     a vector. Due to implementation limitations, "A"
     must be a quadratic matrix and "b" must have a
     dimension that is equivalent to that of "A".
     Furthermore, the determinant of "A" must be
     non-zero. The algorithm used is devised from
     Cramer's Rule and thus inefficient. The preferred
     algorithm for this task is Gaussian Elimination.
     If you have a matrix and a vector of real numbers,
     please consider using either Math::MatrixReal or
     Math::Pari instead.

    [download]

    Note that the complexity of Cramer's rule is O(n!) whereas Gaussian Elimination works in O(n^3) as noted by others in this thread.

    Steffen

[reply]
[d/l]
Re: Solving Simultaneous Equations with Matrices
by Anonymous Monk on Mar 09, 2014 at 10:39 UTC
    sir please tell how i can solve for the value of three variable using crammers rule for ex. i have equation 9a+3b+c=6 and so on.. or if i want to pass the value of with the help of a txt file instead of 9;3;1;6: please help
[reply]

Back to Meditations