<- previous index next ->
Simultaneous equations are multiple equations involving the same
variables. In general, to get a unique solution we need the
same number of equations as the number of unknown variables,
and the equations mutually linearly independent.
A sample set of three equations in three unknowns is:
eq1: 2*x + 3*y + 2*z = 13
eq2: 3*x + 2*y + 3*z = 17
eq3: 4*x - 2*y + 2*z = 12
One systematic solution method is called the Gauss-Jordan reduction.
We will reduce the three equations such that eq1: has only x,
eq2: has only y and eq3: has only z making the constant yield
the solution. Operations are always based on the latest version
of each equation. The numeric solution will perform the same
operations on a matrix.
Reduce the coefficient of x in the first equation to 1, dividing by 2
eq1: 1*x + 1.5*y + 1*z = 6.5
Eliminate the variable x from eq2: by subtracting eq2 coefficient of
x times eq1:
eq2: becomes eq2: - 3* eq1:
eq2: (3-3)*x + (2-4.5)*y + (3-3)*z = (17-19.5) then simplifying
eq2: 0*x - 2.5*y + 0*z = -2.5
Eliminate the variable x from eq3: by subtracting eq3 coefficient of
x times eq1:
eq3: becomes eq3: - 4* eq1:
eq3: (4-4)*x (-2 -6)*y + (2-4)*z = (12-26) then becomes
eq3: 0*x - 8*y - 2*z = -14
The three equations are now
eq1: 1*x + 1.5*y + 1*z = 6.5
eq2: 0*x - 2.5*y + 0*z = - 2.5
eq3: 0*x - 8.0*y - 2*z = -14.0
Reduce the coefficient of y in eq2: to 1, dividing by -2.5
eq2: 0*x + 1*y + 0*z = 1
Eliminate the variable y from eq1: by subtracting eq1 coefficient of y times eq2:
eq1: becomes eq1: -1.5* eq2:
eq1: 1*x + 0*y + 1*z = (6.5-1.5) then becomes
eq1: 1*x + 0*y + 1*z = 5
Eliminate the variable y from eq3: by subtracting eq3 coefficient of y times eq2:
eq3: becomes eq3: -8* eq2:
eq3: 0*x + 0*y - 2*z = (-14 +8) then becomes
eq3: 0*x + 0*y - 2*z = -6
The three equations are now
eq1: 1*x + 0*y + 1*z = 5
eq2: 0*x + 1*y + 0*z = 1
eq3: 0*x + 0*y - 2*z = -6
Reduce the coefficient of z in eq3: to 1, dividing by -2
eq3: 0*x + 0*y + 1*z = 3
Eliminate the variable z from eq1: by subtracting eq1 coefficient of z times eq2:
eq1: becomes eq1: -1* eq2:
eq1: 1*x + 0*y + 0*z = (5-3) then becomes
eq1: 1*x + 0*y + 0*z = 2
Eliminate the variable z from eq2: by subtracting eq2 coefficient of z times eq2:
eq2: becomes eq3: 0* eq2:
eq2: 0*x + 1*y + 0*z = 1
The three equations are now
eq1: 1*x + 0*y + 0*z = 2 or x = 2
eq2: 0*x + 1*y + 0*z = 1 or y = 1
eq3: 0*x + 0*y + 1*z = 3 or z = 3 the desired solution.
The numerical solution simply places the values in a matrix and uses the same
reductions shown above.
Given the equations:
eq1: 2*x + 3*y + 2*z = 13
eq2: 3*x + 2*y + 3*z = 17
eq3: 4*x - 2*y + 2*z = 12
Create the matrix | 2 3 2 13 | having 3 rows and 4 columns
A = | 3 2 3 17 |
| 4 -2 2 12 |
The following code, using n=3 for these three equations, computes the
same desired solution as the manual method above.
for k=1:n
for j=k+1:n
A(k,j) = A(k,j)/A(k,k)
end
for i=1:n
if(i not k)
for j=1:n
A(i,j) = A(i,j) - A(i,k)*A(k,J)
Now we must consider the possible problem of A(k,k) being zero
for some value of k. It is rather obvious that the order of the
equations does not matter. The equations can be given in any
order and we get the same solution. Thus, simply interchange
any equation where we are about to divide by a zero A(k,k)
with some equation below that would not result in a zero A(k,k).
It turns out that we get better accuracy if we always pick
the equation that has the largest absolute value for A(k,k).
If the largest value turns out to be zero then there is no
unique solution for the set of equations.
We generally want numerical code to run efficiently and thus we
will not physically interchange the equations but rather keep
a row index array that tells us where the k th row is now.
The code for the final algorithm is given in the links below.
The instructor understands that some students have a strong
prejudice in favor of, or against, some programming languages.
After about 50 years of programming in about 50 programming
languages, the instructor finds that the difference between
programming languages is mostly syntactic sugar. Yet, since
students may be able to read some programming languages
easier than others, these examples are presented in "C",
Fortran 95, Java and Ada 95. The intent was to do a quick
translation, keeping most of the source code the same,
for the different languages. Style was not a consideration.
Some rearranging of the order was used when convenient.
The numerical results are almost exactly the same.
The same code has been programmed in "C", Fortran 95, Java and Ada 95
as shown below with file types .c, .f90, .java and .adb:
simeq.c "C" language source code
simeq.h "C" header file
simeq.f90 Fortran 95 source code
simeq.java Java source code
simeq.adb Ada 95 source code
real_arrays.ads Ada 95 source code
real_arrays.adb Ada 95 source code
simeq.m MATLAB language source code
simeq_m.out MATLAB output
Throughout this course you will see variations of this source code,
tailored for specific applications. The packaging will change with
"C" files having code inside with 'static void', Fortran 95 code using
modules and, Java and Ada code using packages.
It should be noted that the algorithm is exactly the same for sets
of equations with complex values. The code change is simply
changing the type in Fortran 95, Java, and Ada 95. The Java class
'Complex' is on my WEB page. The "C" code requires a lot of changes.
I wrote the first version of this program for the IBM 650 in assembly
language as an electrical engineering student. The program was for
complex values and solved for node voltages in alternating current
circuits. A quick and dirty version is ac_circuit.java
that needs a number of java packages:
Matrix.java
Complex.java
ComplexMatrix.java
ac_analysis.java an improved version
Then an even more complete version that plots up to eight node voltages.
ac_plot.java simple Java plot added
ac_plot.dat capacitor, inductor and tuned circuits
Output of java myjava.ac_plot.java < ac_plot.dat
There are systems of equations with no solutions:
eq1: 1*x + 0*y + 0*z = 2
eq2: 2*x + 0*y + 0*z = 2
eq3: 4*x - 2*y + 3*z = 5
Some may ask: What about solving |A| * |X| = |Y| for X, given A and Y
using |X| = |Y| * |A|^-1 (inverse of matrix A) ?
The reason this is not a good numerical solution is that slightly
more total error will be in the inverse |A|^-1 and then a little
more error will come from the vector times matrix multiplication.
The code for matrix inverse is very similar to the code for solving
simultaneous equations. Added effort is needed to find the
maximum pivot element and there must be both row and column
interchanges. An example that shows the increasing error with the
increasing size of the matrix, on a difficult matrix, is shown below.
Note that results of a 16 by 16 matrix using 64-bit IEEE Floating
point arithmetic that is ill conditioned may become useless.
inverse.f90
test_inverse.f90
test_inverse_f90.out
Extracted form test_inverse_f90.out is
initializing big matrix, n= 2 , n*n= 4
sum of error= 1.84748050191530E-16 , avg error= 4.61870125478825E-17
initializing big matrix, n= 4 , n*n= 16
sum of error= 2.19971263426544E-12 , avg error= 1.37482039641590E-13
initializing big matrix, n= 8 , n*n= 64
sum of error= 0.00000604139304982709 , avg error= 9.43967664035483E-8
initializing big matrix, n= 16 , n*n= 256
sum of error= 83.9630735209012 , avg error= 0.327980755941020
initializing big matrix, n= 32 , n*n= 1024
sum of error= 4079.56590417946 , avg error= 3.98395107830025
initializing big matrix, n= 64 , n*n= 4096
sum of error= 53735.8765782488 , avg error= 13.1191104927365
initializing big matrix, n= 128 , n*n= 16384
sum of error= 85784.2643647822 , avg error= 5.23585597929579
initializing big matrix, n= 256 , n*n= 65536
sum of error= 1097119.16168229 , avg error= 16.7407098645368
initializing big matrix, n= 512 , n*n= 262144
sum of error= 1.36281435213093E+7 , avg error= 51.9872418262837
initializing big matrix, n= 1024 , n*n= 1048576
sum of error= 1.24247404738082E+9 , avg error= 1184.91558778841
Very similar results from the C version:
inverse.c
test_inverse.c
test_inverse.out
Extracted form test_inverse.out is
initializing big matrix, n=1024, n*n=1048576
sum of error=1.24247e+09, avg error=1184.92
A case study using 32-bit IEEE floating point and 50, 100, and 200
digit multiple precision are shown in Lecture 3a
<- previous index next ->