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.
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
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");
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
###########################################################
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";
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 ]
############################################################
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");
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
#######################################################
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";
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 ]
########################################################
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]
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";
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 ]
#################################################################################
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 :-) )