Układ równań 200x200

Witam

Piszę układ równań 200x200 i chce go rozwiązać dwiema metodami: LU i Gaussa Seidla ale coś mi tu nie gra:

Aij=1/3
Aii=2.2+2*i

Bi=suma ( po j ) Aij*j

algorytm LU:

 #include <iostream>
#include<fstream>
using namespace std;

int main()
{
    const int W=200,K=200; double A[W][K]; double L[W][K]; double U[W][K]; double Y[W];double X[W]; double B[W];
    int i,j,k; double w;
    for(i=0;i<W;i++) //Wprowadzanie danych do macierzy g³ównej
    {
        for(j=0;j<K;j++)
        {
            if(i==j)A[i][j]=1;
            else A[i][j]=(2.2+2*2*i)*3;
        }
    }
    for(i=0;i<W;i++) //Macierz L
    {
        for(j=0;j<K;j++)
        {
            if(i==j)L[i][j]=1;
        }
    }
    for(i=0;i<W;i++) //Macierz U
    {
        for(j=0;j<K;j++)
        {
          U[i][j]=0;
        }
    }
    for(i=0;i<W;i++)
    {
        for(j=0;j<K;j++)
        {
            double s=0;
            for(k=0;k<i;k++) //suma
            {
                s+=L[i][k]*U[k][j];
            }
            U[i][j]=A[i][j]-s;
        }
        for(j=0;j<K;j++)
        {
            double s=0;
            for(k=0;k<i;k++) //suma
            {
                s+=L[j][k]*U[k][i];
            }
            L[j][i]=(A[j][i]-s)/U[i][i];
        }

    }
    for(i=0;i<W;i++)
    {
        double s=0;
        for(j=0;j<K;j++)
        {
         s+=A[i][j]*j;
        }
        B[i]=3*s;
    }
    for(i=0;i<W;i++)
    {

        double s=0;
        for(j=0;j<i;j++)
        {
          s+=L[i][j]*Y[j];

        }
     Y[i]=B[i]-s;
    }
    for(i=W-1;i>=0;i--)
    {
        double s=0;
        for(j=i+1;j<K;j++)
        {
          s+=U[i][j]*X[j];
        }
        X[i]=(Y[i]-s)/U[i][i];
    }
    //for(i=0;i<W;i++)
    //{
       //cout<<""<<X[i]<<endl;
    //}
         std::ofstream plik("wynik_niewiadome.txt");
         for(i=0; i<W; i++)
    {

           plik<<X[i]/3<<"\n\n";

         }
plik.close();
return 0;
}

Czy to jest poprawnie realizowany algorytm??

    for(i=0;i<W;i++) //Macierz L
    {
        for(j=0;j<K;j++)
        {
            if(i==j)L[i][j]=1;
        }
    }

Na przekątnej wstawiasz jedynki, lecz pozostałe pola nie zostały niewypełnione - posiadają losowe wartości.

Wypełnić je zerami?

Dziękuję za szybką odpowiedź. W metodzie Gaussa Seidla jest wszystko ok? Na wasze oko? Iteracji przyjąłem 10.

 #include <iostream>
#include <math.h>
#include <fstream>

using namespace std;

const double eps=0.01;

int main()
{
    int const W=200,K=200; double A[W][K],B[W],X[W],G[W], XX[W]; int i,j;

    for(i=0;i<W;i++)
    {
        for(j=0;j<K;j++)
        {
            if(i==j)A[i][j]=2.2+2*i;
            else A[i][j]=1/3;
        }
    }


    for(i=0;i<W;i++)
    {
        double s=0;
        for(j=0;j<W;j++)
        {
            s+=A[i][j]*j;
        }
        B[i]=s;
    }



    for(i=0;i<W;i++) X[i]=0;

    //Wykonujemy iteracje:

//do{
for(int k=1;k<10;k++)
{


    for(i=0;i<W;i++)
    {
        double s=0;
        for(j=0;j<W;j++)
        {
         if(j>i) s+=A[i][j]*X[j];
         if(j<i) s+=A[i][j]*X[j];
        }
        XX[i]=X[i]; //do warunku stopu

        X[i]=(B[i]-s)/A[i][i];
     }
}

         std::ofstream plik("wynik.txt");
         for(i=0; i<W; i++)
         {

           plik<<X[i]<<"\n\n";;
         }
plik.close();

return 0;
}

Poczatkujacy21 napisał(a)

W metodzie Gaussa Seidla jest wszystko ok?
Zrób test to się dowiesz.

Poczatkujacy21 napisał(a)

Iteracji przyjąłem 10.
Tak się nie robi. Zobacz http://pl.wikipedia.org/wiki/Metoda_Gaussa-Seidla#Warunek_zako.C5.84czenia_iteracji