Numerical Methods: Gauss Seidel

• Tweet
• 

        This solves systems of equations using the Gauss Seidel Method. 

/* 
 * File:   main.cpp
 * Author: Bangonkali
 *
 * Created on March 23, 2012, 8:14 PM
 */
#include <iostream>
#include <cmath>
#include <cstdlib>


// example taken from wikipedia 
// http://j.mp/GV1mPz

using namespace std;

void printMatrix(int N, double A[20][21]) {
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j <= N; j++)
        {
            cout << A[i][j] << '\t';
        }
        
        cout << endl;
    }
    
    cout << endl;
}

int gauss_seidel(int N, double A [20] [21], double r[20], double error[20]) {
    double r_old[N-1];
    
    for (int i = 0; i < N; i++) {
        r[i] = 0; // make initial guess to 0
        r_old[i] = 0; // make initial guess to 0
        error[i] = 0; // make all errors to 0
    }
    
    double max_absolute_error = 0;
    int iter = 0;
    for (iter = 0; iter < 1E12; iter++) {
        // get each iteration result
        max_absolute_error = 0;
        for (int i = 0; i < N; i++) {
            r_old[i] = r[i];
            r[i] = A[i][N] / A[i][i];
            for (int j = 0; j < N; j++) {
                if (j == i ) continue;
                r[i] = r[i] - (A[i][j] * r[j]) / A[i][i];
            }
            
            error[i] = abs((r[i] - r_old[i]) / r[i]) * 100;
            if (error[i] > max_absolute_error) {
                max_absolute_error = error[i];
                
            }
        }
//        for (int i=0; i<3; i++) cout << r[i] << "  ";
//        cout << "max_absolute_error = " << max_absolute_error << endl;
        if (max_absolute_error < 1E-12) {
            cout << "error break" << endl;
            return iter;
        }
    }
    
    return iter;
}
/*
 * 
 */
int main(int argc, char** argv) {
    double A[20][21];
    double X[20];
    double E[20];
    
    bool err;

    A[0][0] = 7; 
    A[0][1] = -2; 
    A[0][2] = 1; 
    A[0][3] = 2;
    A[0][4] = 3;
    
    A[1][0] = 2; 
    A[1][1] = 8; 
    A[1][2] = 3; 
    A[1][3] = 1; 
    A[1][4] = -2; 
    
    A[2][0] = -1;
    A[2][1] = 0; 
    A[2][2] = 5; 
    A[2][3] = 2; 
    A[2][4] = 5; 
    
    A[3][0] = 0;
    A[3][1] = 2; 
    A[3][2] = -1; 
    A[3][3] = 4; 
    A[3][4] = 4; 

    printMatrix(4,A);
    int iterations = gauss_seidel(4, A, X, E);
    for (int i=0; i<3; i++) cout << X[i] << "  ";
    cout << endl;
    for (int i=0; i<3; i++) cout << E[i] << "  ";
    cout << endl;
    cout << "iterations = " << iterations << endl;
    return 0;
}

Tags:
• gauss
• Methods
• Numerical
• seidel
• systems