Układ równań - met Gaussa

0

Witam wszystkich... mam następujący problem:
Chcę rozwiązać przykładowy układ równań met Gaussa:

X1 + X4 = 0
X1 + X5 = 0
X1 + X7 = 14
X2 + X6 = 0
X3 + X6 = 0
X3 + X7 = 0

Ponieważ w takiej sytuacji układ ma nieskończenie wiele rozwiązań za X1 podstawiam 0.

Korzystam z tego znalezionego programu:

// Program rozwiązuje układ równań liniowych o n niewiadomych
// za pomocą metody eliminacji Gaussa.
//-----------------------------------------------------------
// (C)2006 mgr J.Wałaszek                     I LO w Tarnowie

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>

using namespace std;

const double EPS   = 0.0000000001; // dokładność porównania z zerem
const int    MAXEQ = 100;          // maksymalna ilość równań w układzie

// Funkcja dokonująca eliminacji niewiadomych. Jeśli operacja
// się powiedzie, zwraca true. Inaczej zwraca false.

bool EliminujX(int n, double AB[][MAXEQ+1])
{
  int    i,j,k;
  double m;

  for(i = 0; i < n - 1; i++)
  {
    if(fabs(AB[i][i]) < EPS) return false;
    for(j = i + 1; j < n; j++)
    {
      m = -AB[j][i] / AB[i][i];
      for(k = i + 1; k <= n; k++) AB[j][k] += m * AB[i][k];
    }
  }
  return true;
}

// Funkcja oblicza kolejne niewiadome x z macierzy AB
// przetworzonej przez funkcję Eliminuj_X().
// Jeśli operacja się powiedzie, zwraca true. Inaczej
// zwraca false.

bool ObliczX(int n, double X[], double AB[][MAXEQ+1])
{
  int    i,j;
  double s;

  for(i = n - 1; i >= 0; i--)
  {
    if(fabs(AB[i][i]) < EPS) return false;
    s = AB[i][n];
    for(j = n - 1; j > i; j--) s -= AB[i][j] * X[j];
    X[i] = s / AB[i][i];
  }
  return true;
}

//-----------------------------------------------------
// Program główny
//-----------------------------------------------------

int main(int argc, char* argv[])
{
  ifstream fin;
  ofstream fout;
  int i,j,n;
  double AB[MAXEQ][MAXEQ+1], X[MAXEQ];

  cout.precision(5);     // 5 cyfr po przecinku
  cout.setf(ios::fixed); // format stałoprzecinkowy

  cout << "Uklad rownan liniowych  - metoda eliminacji Gaussa\n"
          "--------------------------------------------------\n"
          "(C)2006 mgr Jerzy Walaszek         I LO w Tarnowie\n\n";

// Dane dla programu odczytujemy z pliku tekstowego o nazwie in.txt,
// który musi się znajdować w tym samym katalogu co program
// Pierwszy wiersz pliku powinien zawierać liczbę n
// Następne n kolejnych wierszy powinno zawierać współczynniki ai dla
// tego wiersza, a na końcu współczynnik bi. Kolejne współczynniki
// oddzielone są od siebie przynajmniej jedną spacją.

  fin.open("in.txt");
  fin >> n;
  if(n <= MAXEQ)
  {
    for(i = 0; i < n; i++)
      for(j = 0; j <= n; j++) fin >> AB[i][j];
    fin.close();

// Dokonujemy eliminacji oraz obliczania niewiadomych x

    cout << "\n--------------------------------------------------\n"
            "Wyniki:\n\n";

    if(EliminujX(n,AB) && ObliczX(n,X,AB))
    {
      fout.open("out.txt");
      for(i = 0; i < n; i++)
      {
        cout << "x" << i + 1 << " = " << setw(12) << X[i] << endl;
        fout << X[i] << endl;
      }
      fout.close();
    }
    else cout << "Rozwiazanie ukladu rownan nie powiodlo sie\n";
  }
  else
  {
    fin.close();
    cout << "Zbyt wiele rownan!\n";
  }
  cout << "\n--------------------------------------------------\n";
  system("pause");
  return 0;
}

tak wygląda plik wejściowy:

6
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 0 1 14
1 0 0 0 1 0 0
0 1 0 0 1 0 0
0 1 0 0 0 1 0

czyli taki układ równań:

0X2 + 0X3 + 1X4 + 0X5 + 0X6 + 0X7 = 0
0X2 + 0X3 + 0X4 + 1X5 + 0X6 + 0X7 = 0
0X2 + 0X3 + 0X4 + 0X5 + 0X6 + 1X7 = 14
1X2 + 0X3 + 0X4 + 0X5 + 1X6 + 0X7 = 0
0X2 + 1X3 + 0X4 + 0X5 + 1X6 + 0X7 = 0
0X2 + 1X3 + 0X4 + 0X5 + 0X6 + 1X7 = 0

Próbowałem również innym: http://www.programmersheaven.com/download/29309/download.aspx
i też zima...

Nie rozumiem dlaczego nie rozwiązuje takiego prostego układu... na innych działa, przykład:

5
2 -2 2 -7 6 -4
7 -3 -2 7 2 11
2 2 -1 1 4 9
9 8 -2 2 -2 21
4 8 -3 3 -1 16

Macie pomysł dlaczego?

Pozdrawiam.

1

Krótki kurs myślenia by Bury Pajac.

Rozumiem, że pojawia ci się następujący komunikat:

cout << "Rozwiazanie ukladu rownan nie powiodlo sie\n";

Popatrzmy gdzie się znajduje ta instrukcja:

    if(EliminujX(n,AB) && ObliczX(n,X,AB))
    {
[..]
    }
    else cout << "Rozwiazanie ukladu rownan nie powiodlo sie\n";

czyli jeżeli EliminujX lub ObliczX zwróci false to wtedy dostajemy komunikat.
Looknijmy kiedy EliminujX i ObliczX zwracają false:

bool EliminujX(int n, double AB[][MAXEQ+1])
{
  [..]
  for(i = 0; i < n - 1; i++)
  {
    if(fabs(AB[i][i]) < EPS) return false; //<-------Tutaj, jedyny return false z tej funkcji
[..]
  }
  return true;
}
bool ObliczX(int n, double X[], double AB[][MAXEQ+1])
{
[..]
  for(i = n - 1; i >= 0; i--)
  {
    if(fabs(AB[i][i]) < EPS) return false; //<------- I tutaj, też jedyny return false z tej funkcji
[..]
  }
  return true;
}

Czyli w obu przypadkach warunek jest taki:

if(fabs(AB[i][i]) < EPS) return false;

,gdzie

 
const double EPS   = 0.0000000001; // dokładność porównania z zerem

To jak interpretować ten warunek?
Jeżeli na przekątnej jest wartość np. zero, to się nie udało rozwiązać układu.
A czym karmisz tego proga?

6
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 0 1 14
1 0 0 0 1 0 0
0 1 0 0 1 0 0
0 1 0 0 0 1 0

O, zera na przekątnej.

Popatrzmy kto jest autorem:

// (C)2006 mgr J.Wałaszek                     I LO w Tarnowie

Gość z wyższym wykształceniem i do tego uczy w szkole, czyli jest po jakiejś szkole pedagogicznej, fuck, to nie można go podać jako przykład, iż uczelnie techniczne nie uczą, bądź uczą źle. Chociaż jest możliwość, że gość skończył politechnikę, a potem podyplomowo zrobił pedagogikę.

0

Dzięki serdeczne za wyjaśnienia.

Szkoda, bo chciałem użyć tego kodu w programie rozwiązującym zagadnienie transportowe... w którym niemalże występują same tego typu równania.

Masz może (lub ktoś inny ma) jakieś wskazówki gdzie dalej szukać rozwiązania tego typu problemu? Może jakaś inna metoda?
Wydaje mi się że mógłbym sam napisać funkcje, która to rozwiązuje ale najpierw chciałbym poszukać i poznać gotowe rozwiązania jeśli istnieją (i jakoś lepiej spożytkować ten czas).

Jeszcze raz dziękuję za odp.

1
beny6.6.6 napisał(a)

Szkoda, bo chciałem użyć tego kodu w programie rozwiązującym zagadnienie transportowe... w którym niemalże występują same tego typu równania.

Wciąż możesz go uzyć, wystarczy naprawić to co jest spieprzone. Jak?
Problemem jest to, że ten kod potrzebuje mieć niezerowe wartości na przekątnej, a skoro potrzebuje to wystarczy mu to zapewnić. Znowu pytanie, jak?
Jeżeli na przekątnej jest zero, to niech znajdzie wiersz w którym ta wartość w tej kolumnie jest poprawna i niech te wiersze zamieni, czyli:

#include <cstring> 
//dla memcpy
[..]
bool fixit(int n, double AB[][MAXEQ+1],int start){ // start jest numerem wadliwego wiersza, trzeba wymyśleć inną nazwę
   for(int j=start+1;j<n;++j)       // od następnego wiersza do ostaniego
      if(fabs(AB[j][start])>EPS)     // jeżeli spełnia warunki zamień ten wiersz z tym wadliwym
      {    
	  double * temp = new double[n+1];
          size_t size_in_bytes =(n+1)*sizeof(double);

	  memcpy(temp,AB[start],size_in_bytes);
	  memcpy(AB[start],AB[j],size_in_bytes);
	  memcpy(AB[j],temp,size_in_bytes);

	  delete[] temp;
	  return true;                    // i zakończ funkcję
	}
    return false;                        // jeżeli w kolumnie są same zera, to równania nie można naprawić
} 
bool EliminujX(int n, double AB[][MAXEQ+1])
{
  int    i,j,k;
  double m;

  for(i = 0; i < n - 1; i++)
  {
    if (fabs(AB[i][i])<EPS)    //jeżeli na przekątnej jest zero
      if(!fixit(n,AB,i))           // i jeżeli nie można tego naprawić
	return false;             // to rozwiązywanie układu się nie powiodło
                                      // w przeciwnym przypadku jesteśmy dalej w grze
    for(j = i + 1; j < n; j++)
    {
      m = -AB[j][i] / AB[i][i];
      for(k = i + 1; k <= n; k++) AB[j][k] += m * AB[i][k];
    }
  }
  return true;
}

Ale to tylko głupia łata na jeden bug, a może ich być więcej.

beny6.6.6 napisał(a)

Wydaje mi się że mógłbym sam napisać funkcje, która to rozwiązuje ale najpierw chciałbym poszukać i poznać gotowe rozwiązania jeśli istnieją (i jakoś lepiej spożytkować ten czas).

Może gsl? http://www.gnu.org/software/gsl/manual/html_node/Linear-Algebra-Examples.html

1

mozna numerical recipes albo BLAS. Ten pierwszy jest banalny i jest z przykladami w c++/c/fortran/pascal w necie. A i sa porty na AS3 i c# i java

0

Dzięki za pomoc, chyba ten program + funkcja zamieniająca wiersze będzie dla mnie najlepszym rozwiązaniem.

Pozdrawiam.

1 użytkowników online, w tym zalogowanych: 0, gości: 1