Hi,
CORRECTLY COMBINING EQUATIONS
When you add two equations you have to add *ALL* the columns not just the ones on the right or the ones on the left or some other subset.
For example, if you add the two equations:
y=-6x+1
y=-2x+2
you do not get:
y=-8x+3
you actually get:
2y=-8x+3
and dividing both sides by 2 we get:
y=-4x+1.5
and plugging in x=-0.25 we get:
y=2.5
So adding or subtracting equations means to add the ENTIRE equation to the other, not just some terms.
GENERAL IDEA OF GAUSSIAN ELIMINATION
The whole idea with GE is to reduce the equation coefficients to a triangular matrix so that the 'last' equation has a solution like:
x=2
To do this, we multiply some equations by a constant so that when we add or subtract them to another row that row develops a zero coefficient. Once N-1 variables have zero coefficients in any one row we have a solution for that one variable and then we can use back substitution to get the solutions for the other variables.
The coefficient matrix for the set:
A*x+B*y=R1
C*x+D*y=R2
is simply:
|A,B|
|C,D|
and the constant vector is:
|R1,R2|
We want to get the matrix in the form:
|A,B|
|0,E|
where E came from some operations on D using A, B and C. The zero in the matrix means that the last equation was reduced to one unknown
Ey=F
so y is simply:
y=F/E
and then we can use back substitution to get x.
NUMERICAL STABILITY
The whole idea with numerical stability stems from the fact that the usual numerical processor has limited precision. This is true of any computer really because there is an approximation used to represent the numbers due to the fact that we can not store the entire number. For example it would take an infinite number of digits to store the true value of pi. And this is true of other numbers too like 1/3 which is:
0.333333333333333333333333333 etc.
Now storing a number like 1.001 would not be a problem really (it may come down to an approximation too but lets say it doesnt for now). But when we go to MULTIPLY it with another number and then use it with still yet another calculation it can quickly lead to errors.
So lets say we have that 1.001 and we want to multiply it times another number 1.001 (happens to be the same to make this easy).
We get:
1.002001
But notice that we started with 4 digits, and we ended up with 7 digits (that's almost twice as many digits). Well what would happen if we only had four digits (this simulates a limited storage type like float or double but with even more limited precision just for the illustration). If we only had four digits to work with all along, we would have to store that number as:
1.002
so we would have lost some precision. The error is small 0.000001, but it is nonetheless still ever present. And lets say we multiply that by 1.001 again, we get:
1.003002
but have to store it as:
1.003
Now if we had a larger storage register we would (with all the calculations so far) have gotten:
1.003003001
so now the error is:
0.000003001
when before it was only:
0.000001
so we can see the error grew.
But the real problem comes in when we go to add or subtract, especially subtract two numbers with built in errors from previous calculations.
If we subtract two numbers that are close to each other, we end with with a significant portion of the result being the error itself. This leads to gross errors sometimes.
THE PIVOT
Lets take a quick look at a simple set of equations and use GE to solve them. The best way to list them is with the variables on one side and the constants on the other side like so:
2x+4y=8
4x+2y=10
Now we want to use GE to solve this, so we start by noting that we want to end up with a form like this:
2x+4y=8
0x+Ay=B
where we have a zero before the x in the last equation.
Noting that if we take the first equation and multiply it by 2 we get the x coefficient the same as the x coeff in the last equation, so doing that:
2x+4y=8 times (2) gives us:
4x+8y=16
Now we replace the top equation with that and get:
4x+8y=16
4x+2y=10
and then subtract the bottom one from the top and get:
0x+6y=6
Now we've eliminated the x coefficient so the solution for y is very clear:
6y=6
so
y=1
Now we go back to the first equation:
2x+4y=8
and substitute y=1 into it and get:
2x+4=8
and so
2x=8-4
so:
2x=4
so:
x=2
and that's that.
But what if the original set was written as (i know this is trivial but just for illustration):
0x+6y=6
2x+4y=8
Now we can not get the equations into the form:
Ax+By=R1
0x+Dy=R2
because whatever we multiply the coefficient in the first equation by we always come out with 0 for the x coeff so it will never change the second equation at all and this means we can not get a zero coefficient in the second equation this way.
When this happens we have to 'rotate' the whole 'matrix' of coefficients so that we push the 0 to one side and replace it with a non zero number:
By+Ax=R1
Dy+0x=R2
Now we have B which is non zero so we can start the process. We note that D=B*a where a is some number so when we subtract we force D to be zero. In this very simple example we would not have to actually do this, but this is the way it works when we run into this situation and actually have to perform these steps.
Now there is the question of numerical stability while doing this task. We are forced to use number crunchers that have limited precision so this introduces some new problems.
[LATER]
A NUMERICAL INVESTIGATION OF THE EFFECT OF MAGNITUDE OF THE PIVOT ON THE RESULT
I was looking at the equations (Set 1):
1.001x+2.001y=3.002
1.002x+2.002y=3.004
and comparing these to the equations with rotated x and y columns (Set 2):
2.001y+1.001x=3.002
2.002y+1.002x=3.004
The idea here is that because the pivot in Set 1 (1.001) is lower than the pivot in Set 2 (2.001) the result of the GE calculated variable in Set 2 will be numerically more stable than the GE calc'd var in Set 1 which means it will be a better approximation.
To check this out using the modern computer however we have to impose a rounding function like:
v=floor(N*v)/N
where v is the result of a calculation. We have to do this to see the effect rounding has on an actual calculation using more digits where it happens automatically within the numeric processor.
With N=1000 we see some interesting results. This limits the precision of the calculations so we end up seeing the effect on the solution set x and y. For Set 1 i was seeing y=0.333 which is not correct at all, off by quite a bit, and for Set 2 i saw x=1 which is correct.
I want to go over this to make sure i did everything right, but you get the general idea. Supposedly we get better results with a larger pivot.
Also written in some literature which i havent investigated fully yet is that the stability also gets better with the introduction of a scaling factor where s=sqrt(Akk) where Akk are the coefficients along the diagonal. What i had found so far was that any introductions of any scaling factor meant that there was yet another multiplication involved (that's what we have to do to all the coefficients first) so the accuracy always seemed to get worse. My only guess at this point was that the scaling factor was introduced when computers could only do fixed point precision. With floating point it actually looks like it gets worse.
I havent looked at tons of matrixes with this in mind, but have looked a quite a few so far and all the results either do not change at all or get worse. So i cant recommend a scaling factor unless you're using fixed point (such as maybe with a microcontroller).