2012-12-09 13:38:40Morris

[UVA][高斯消去][線性系統] 10109 - Solving Systems of Linear ...

 
Solving Systems of Linear Equations

The Problem

You may have solved linear equation early in the school. Problems involving solving sets of linear equation are very important in the field of Engineering and Mathematics.
Let us consider that we have a system of linear equations
                 a11x1 + a12x2 + a13x3 = c1
                 a21x1 + a22x2 + a23x3 = c2
                 a31x1 + a32x2 + a33x3 = c3
We can solve it by reducing technique:
            Step 1:
            x1 + a12/a11x2 + a13/a11x3 = c1/a11
            a21x1 + a22x2 + a23x3 = c2
            a31x1 + a32x2 + a33x3 = c3

            Step 2:
            x1 + a12/a11x2 + a13/a11x3 = c1/a11
                    (-a21* a12/a11)a22x2 +  (-a21* a13/a11)a23x3 =  (-a21* c1/a11)c2
                    (-a31* a12/a11)a22x2 +  (-a31* a13/a11)a23x3 =  (-a31* c1/a11)c2

Now do as step 1 for second row and so on.
This can be made more effective using matrix method. The set of equation for `n' unknowns can be written as
                 a11x1 + a12x2 + a13x3 + ... + a1nxn = c1
                 a21x1 + a22x2 + a23x3 + ... +  a2nxn = c2
                 a31x1 + a32x2 + a33x3 + ... +  a3nx = c3
                         ...   ...   ...   ...   ...   ...   ...   ...   ...   ...   ...   ...   ...   ...
                 am1x1 + am2x2 + am3x3 + ... +  amnx = cn

In matrix form
Equations in matrix form: [A] * {X} = {C}
Compactly        [A] * {X} = {C}
From this we can solve values of X's. The matrix [AC] is called an augmented (see example below) matrix. If after elimination process the rank of matrix [A] and rank of matrix [AC] not equals, the system is called inconsistent and it does not have a solution. If the matrix is consistent and number of unknowns is greater then rank of matrix then the matrix system has arbitarily many solutions containing (NumberOfUnknowns-rank) arbitary constants. Rank of a matrix is defined as the number of non zero rows of a matrix system. Otherwise if the rank and number of unknows equals then the system has been solved.

For example let a system of equations be
                  9x1 + 4x2 + x3 = -17
                  x1 - 2x2 - 6x3 = 14
                  x1 + 6x2  = 4
This sets of equation can be written as
9    4      1   -17
1  -2   -6      14
1    6      0        4
So the steps involving the solution is

Step : 1
              1             -2             -6             14
              1              6              0              4
              9              4              1            -17

Step : 2
              1             -2             -6             14
              0              8              6            -10
              0             22             55           -143

Step : 3
              1             -2             -6             14
              0              1              3/4             -5/4
              0              0             77/2           -231/2

Step : 4
              1             -2             -6             14
              0              1              3/4             -5/4
              0              0              1             -3

Step : 5
              1             -2              0             -4
              0              1              0              1
              0              0              1             -3

Step : 6
              1              0              0             -2
              0              1              0              1
              0              0              1             -3

x1 = -2
x2 = 1
x3 = -3
Again consider this system
        2x1 + 2x2 + 2x3 = 2
        4x1 + 4x2 + 4x3 = 4
        16x1 + 16x2 + 16x3 = 16
Steps are:
Step : 1
              2              2              2              2
              4              4              4              4
             16             16             16             16

Step : 2
              1              1              1              1
              0              0              0              0
              0              0              0              0
This system has number of unknowns 3 and rank is 1. So this system has arbitarily many solutions containing (3-1) = 2 arbitary constants.
Another system
x + y = 10
x + y = 20
x + y = 50
Steps are:
Step : 1
              1              1             10
              1              1             20
              2              2             50

Step : 2
              1              1             10
              0              0             10
              0              0             30

Step : 3
              1              1             10
              0              0             10
              0              0             30

Step : 4
              1              1             10
              0              0             10
              0              0              1

Step : 5
              1              1              0
              0              0              0
              0              0              1

Step : 6
              1              1              0
              0              0              0
              0              0              1

As rank of [A] (in this case: rank(A) = 1) is not equal to the rank of augmented matrix [AC] (in this case: rank(AC) = 2) , the system has no solution.

However though there are other methods to compute this solution for the matrix system, the main problem occurs are
        1. Round off errors or computational error due to the use of floating point number
        2. Error due to wrong order of the given equation.

To prevent round off error due to floating point number an approach can be used, similar to the process of doing fractional number. So we may use 1/3 as a expression of two integer, the numerator and the denominator, instead of .333333 (with loss of precision). Thus we can prevent this kind of error.

Consider this set of equations
                  5x3 = 10
                  3x2 - 3x3 = 3
                  2x1 - x2 + 2x3 = 7
This set of equations can be written as
0     0     5    10
0     3   -3     3
2   -1     2     7
Now how will you evaluate this matrix without ordering?

The Input

The first line of input is the number of the problem. The next line contains two integers - NumberOfUnknowns and NumberOfEquations (none of these is less then or equal to 0 and greater then 50). The next lines contains the matrix for the system of linear equations. There are number of rows equal to the NumberOfEquations and number of column equal to the NumberOfUnknowns+1. The numbers may be fractional, that is there may be numbers like 1/3 or 6/8. An problem number zero indicates the end of input.

The Output

First print (without the quotation mark)
"Solution for Matrix System # N"
Here `N' is the problem number as taken from input. Then on the next line, for each system of equations output the solution (if exists) expressed in the fractional form in each line. You may assume each of the numerator and denominator part will not exceed the limit of data type long long (64 bit). If there are many solutions as described above print (without the quotation mark)
"Infinitely many solutions containing n arbitrary constants."
(here `n' is the number as described above) , and if there is no solutions print (without the quotation mark)
"No Solution."
.Print a blank line between two systems of linear equations.

Sample Input

1
3 3
9  4  1 -17
1 -2 -6 14
1  6  0  4

2
3 3
2 2 2 2
4 4 4 4
16/1 16/1 16/1 16/1

3
2 3
1 1 10
1 1 20
2 2 50

4
1 1
3 10

0

Sample Output

Solution for Matrix System # 1
x[1] = -2
x[2] = 1
x[3] = -3

Solution for Matrix System # 2
Infinitely many solutions containing 2 arbitrary constants.

Solution for Matrix System # 3
No Solution.

Solution for Matrix System # 4
x[1] = 10/3

Suman Mahbub
01-04-2001

#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
class Frac {
    public:
        long long a, b;
        Frac() {
            a = 0, b = 1;
        }
        Frac(int x, int y) {
            a = x, b = y;
            reduce();
        }
        Frac operator+(const Frac &y) {
            long long ta, tb;
            tb = this->b/gcd(this->b, y.b)*y.b;
            ta = this->a*(tb/this->b) + y.a*(tb/y.b);
            Frac z(ta, tb);
            return z;
        }
        Frac operator-(const Frac &y) {
            long long ta, tb;
            tb = this->b/gcd(this->b, y.b)*y.b;
            ta = this->a*(tb/this->b) - y.a*(tb/y.b);
            Frac z(ta, tb);
            return z;
        }
        Frac operator*(const Frac &y) {
            long long tx, ty, tz, tw, g;
            tx = this->a, ty = y.b;
            g = gcd(tx, ty), tx /= g, ty /= g;
            tz = this->b, tw = y.a;
            g = gcd(tz, tw), tz /= g, tw /= g;
            Frac z(tx*tw, ty*tz);
            return z;
        }
        Frac operator/(const Frac &y) {
            long long tx, ty, tz, tw, g;
            tx = this->a, ty = y.a;
            g = gcd(tx, ty), tx /= g, ty /= g;
            tz = this->b, tw = y.b;
            g = gcd(tz, tw), tz /= g, tw /= g;
            Frac z(tx*tw, ty*tz);
            return z;
        }
    private:
        static long long gcd(long long x, long long y) {
            if(!y)  return x;
            if(x < 0)   x *= -1;
            if(y < 0)   y *= -1;
            long long t;
            while(x%y)
                t = x, x = y, y = t%y;
            return y;
        }
        void reduce() {
            long long g = gcd(a, b);
            a /= g, b /= g;
            if(b < 0)   a *= -1, b *= -1;
        }
};
ostream& operator<<(ostream& out, const Frac&x) {
    out << x.a;
    if(x.b != 1)
        out << '/' << x.b;
    return out;
}
int main() {
    int n, m, i, j, k, N;
    char NUM[100], first = 0;
    long long X, Y;
    Frac matrix[100][100];
    while(scanf("%d", &N) == 1 && N) {
        scanf("%d %d", &n, &m);
        for(i = 0; i < m; i++) {
            for(j = 0; j <= n; j++) {
                scanf("%s", NUM);
                if(sscanf(NUM, "%lld/%lld", &X, &Y) == 2) {
                    matrix[i][j].a = X;
                    matrix[i][j].b = Y;
                } else
                    sscanf(NUM, "%lld", &matrix[i][j].a), matrix[i][j].b = 1;
            }
        }
        Frac tmp, one(1,1);
        int idx = 0, rank = 0;
        for(i = 0; i < m; i++) {
            while(idx < n) {
                int ch = -1;
                for(j = i; j < m; j++)
                    if(matrix[j][idx].a) {
                        ch = j;
                        break;
                    }
                if(ch == -1) {
                    idx++;
                    continue;
                }
                if(i != ch)
                    for(j = idx; j <= n; j++)
                        swap(matrix[ch][j], matrix[i][j]);
                break;
            }
            if(idx >= n) break;
            tmp = one/matrix[i][idx];
            rank++;
            for(j = idx; j <= n; j++)
                matrix[i][j] = matrix[i][j]*tmp;
            for(j = 0; j < m; j++) {
                if(i == j)  continue;
                tmp = matrix[j][i]/matrix[i][idx];
                for(k = idx; k <= n; k++) {
                    matrix[j][k] = matrix[j][k] - tmp*matrix[i][k];
                }
            }
            idx++;
        }
        if(first)   puts("");
        first = 1;
        printf("Solution for Matrix System # %d\n", N);
        int sol = 0;
        for(i = 0; i < m; i++) {
            for(j = 0; j < n; j++) {
                if(matrix[i][j].a)
                    break;
            }
            if(j == n) {
                if(matrix[i][n].a == 0 && sol != 1)
                    sol = 0; // INFINITELY
                else
                    sol = 1; // No Solution.
            }
        }
        if(rank == n && sol == 0) {
            for(i = 0; i < n; i++) {
                printf("x[%d] = ", i+1);
                cout << matrix[i][n] << endl;
            }
            continue;
        }
        if(sol == 1)
            puts("No Solution.");
        else
            printf("Infinitely many solutions containing %d arbitrary constants.\n", n-rank);
    }
    return 0;
}