Metoda Jacobiego przy obliczaniu układu równań

Prosze o pomoc bo nie moge sobie dać rady...

 
#include <iostream>
#include<windows.h>
using namespace std;


void wypelnij_tablice(double** tab, int size);
void wyswietl_tablice(double** tab, int size);
void wypelnij_wektor(double* tab, int size);
void wyswietl_wektor(double* tab, int size);
void wypelnij_diagonalna(double** tab, double** wyj, int size);
void zeruj_tab(double** tab, int size);
void wypelnij_LU(double** tab, double** wyj, int size);
void wypelnij_T(int size, double** LU, double** D, double** T);
void wypelnij_C(int size, double** D, double* B, double* C);
void wypelnij_X1(int size, double** T, double* X, double* C, double* X1);
void color(int n);

int main()
{
/*          tablica A dynamicznie           */ 
    int size;
    cout << "Podaj rozmiar tablicy: ";
    cin >> size;
    double **A;
    A = new double *[size];
    for(int i=0; i<size; i++)
    {
        A[i]= new double[size];
    }
    cout<<"Podaj rownania(tablica A): \n\n";
    wypelnij_tablice(A,size);
    cout<<"Twoja macierz :\n\n";
    color(11);
    wyswietl_tablice(A,size);
color(15);
cout<<"\n";
    double *B = NULL;
    B = new double[size];
    cout<<"Podaj rozwiazania (tablica B):\n\n";
    wypelnij_wektor(B,size);
    cout<<"Twoje odpowiedzi :\n\n";
    color(11);
    wyswietl_wektor(B,size);
    color(15);
    double **D;
    D = new double *[size];
    for(int i=0; i<size; i++)
    {
        D[i]= new double[size];
    }
    
    zeruj_tab(D,size);
    wypelnij_diagonalna(A, D, size);  // diagonalna do potegi -1
    cout<<"Tablica D stworzona z przekatnej podanej talblicy podniesiona do potegi -1 \n\n";
    color(11);
    wyswietl_tablice(D, size);
    color(15);
    double **LU;
    LU = new double *[size];
    for(int i=0; i<size; i++)
    {
        LU[i]= new double[size];
    }
    
    zeruj_tab(LU,size);
    wypelnij_LU(A, LU, size);
    cout<<"\nTablica LU stworzona z podanej tablicy z wylaczeniem jej przekatnej : \n\n";
    color(11);
    wyswietl_tablice(LU, size);
    color(15);
    double **T;
    T = new double *[size];
    for(int i=0; i<size; i++)
    {
        T[i]= new double[size];
    }

    zeruj_tab(T,size);
    wypelnij_T(size, LU, D, T); 
    cout << "\nTablica T(-D^(-1)*LU:" << endl << endl;
    color(11);
    wyswietl_tablice(T, size);
    color(15);
    double *C = NULL;
    C = new double[size];
    
    wypelnij_C(size, D, B, C); 
    cout << " \nTablica C(D*B) = \n" << endl << endl;
    color(11);
    wyswietl_wektor(C, size);
    color(15);
    double *X = NULL;
    X = new double[size];
    
    for(int i=0; i<size; i++)
    {
            X[i] = 1;
    }
    
    
    double *X1 = NULL;
    X1 = new double[size];   
    
    int licznik = 5;
    

    wypelnij_X1(size, T, X, C, X1); cout << endl << "Tablica X1:" << endl;
    cout<<"Wynik: "<<endl;
    color(11);
    wyswietl_wektor(X1, size);
    color(15);
    
    
    


    system("PAUSE");
    return 0;
}
void wypelnij_tablice(double** tab, int size)
{
    for(int i=0; i<size; i++)
    {
        for(int j=0; j<size; j++)
        {
            cout << "Podaj macierz [" << i << "][" << j << "]: ";
            cin >> tab[i][j];
            cout << endl;
        }
    }
}
void wyswietl_tablice(double** tab, int size)
{
    for(int i=0; i<size; i++)
    {
        for(int j=0; j<size; j++)
        {
            cout << tab[i][j] << "  ";
        }
        cout << endl;
    }
}
void wypelnij_wektor(double* tab, int size)
{
        for(int i=0; i<size; i++)
        {
            cout << "Podaj macierz [" << i << "]: ";
            cin >> tab[i];
            cout << endl;
        }
}
void wyswietl_wektor(double* tab, int size)
{
        for(int i=0; i<size; i++)
        {
            cout << tab[i] << endl;
        }
}
void wypelnij_diagonalna(double** tab, double** wyj, int size)
{
    for(int i=0; i<size; i++)
    {
        for(int j=0; j<size; j++)
        {
            if(i==j)
            {
                wyj[i][j] = tab[i][j];
                wyj[i][j] = 1 / wyj[i][j]; 
            }
        }
    }
}
void zeruj_tab(double** tab, int size)
{
  for(int i=0; i<size; i++)
    {
        for(int j=0; j<size; j++)
        {
            tab[i][j] = 0;
        }
    } 
}
void wypelnij_LU(double** tab, double** wyj, int size)
{
    for(int i=0; i<size; i++)
    {
        for(int j=0; j<size; j++)
        {
            if(i!=j)
            {
                wyj[i][j] = tab[i][j];
            }
        }
    }
    for(int i=0; i<size; i++)
     {
             for(int j=0; j<size; j++)
             {
                     wyj[i][j] = wyj[i][j]*(-1);
             }
     }
}

void wypelnij_T(int size, double** LU, double** D, double** T)
{
     for(int i = 0; i < size; ++i)
     {
             for(int j = 0; j < size; ++j)
             {
                     T[i][j] = 0;
                     for(int k = 0; k < size; ++k)
                     T[i][j] += LU[k][j] * D[i][k];
             }
     }

} 
void wypelnij_C(int size, double** D, double* B, double* C)
{
    for(int i=0; i<size; i++)
    {
            for(int j=0; j<size; j++)
            {
                   C[i] = C[i] + D[i][j]*B[i];   //nie wiem czy if dla każdej macierzy zadziała
            }
    }
     
}
void wypelnij_X1(int size, double** T, double* X, double* C, double* X1)
{
        for(int i=0; i<size; i++)
    {
            for(int j=0; j<size; j++)
            {
                    X1[i] = X1[i] + (T[i][j]*X[i]);   
            }
    } 

        for(int i=0; i<size; i++)
    {

                    X1[i] +=  C[i];
    } 
}
void color(int n)
{
    SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE), n);
}

Mam problem z kodem... Dla macierzy 3x3
(najpierw wprowadza sie macierz A potem b)
a=
[3 2 1]
[2 5 7]
[1 2 1]
b=
[11]
[33]
[8]
liczy dobrze... ale dla macierzy
z przykładu z Wikipedii <ort>http://en.wikipedia.org/wiki/Jacobi_method </ort>liczy źle.. i nie wiem gdzie jest błąd
i mam problem jeszcze z tym aby ten algorytm chodził w kółko.
na Wikipedii jest wyjaśnione lepiej...

Zapomniałem dodać że nie moge nigdzie znaleźć informacji na temat kiedy tan układ będzie sprzeczny albo nieoznaczony..
Wiem tylko ze na przekątnej macierzy A nie mogą być 0

Za dużo tego "wypełnij_*" zredukuj kod do tego co dotyczy problemu, żebyśmy nie widzieli ściany kodu (ja widzę u ciebie głównie inne metody).

Ja bym zrobił tak:

void nextJacobiIteration(const double **a, const double *b, const double *x1, double *x2, int n) {
    for(int i=0; i<n; ++i) {
         x2[i] = 0;
         for(int j=0; j<i; ++j)
              x2[i] += a[i][j]*x1[j];
         for(int j=i+1; j<n; ++j)
              x2[i] += a[i][j]*x1[j];
         x2[i] = (b[i] - x2[i])/a[i][i];
    }
}

Dziękuje bardzo zajebiście działa. [PIWO]