Metoda Eliminacji Gaussa C++

[Boczek11]

Witam,

 

Mój problem wygląda następująco. Musze napisać program na rozwiązywanie układu liniowego metodą eliminacji Gaussa w języku C++. Znam algorytm(korzystałem z książki autorstwa D.Kincaid oraz W. Cheney "Analiza numeryczna"),lecz przekształcenie go w program jest dla mnie czarną magią .Bardzo proszę o jakąkolwiek pomoc.

[rezo_]

Moj program napisany w C. Macierz musi byc kwadratowa. W pierwszym wierszu pliku musisz podac rozmiar macierzy. Nazwe pliku w ktorym znajduje sie macierz podajesz w argumencie wywolania programu.

 

#include <stdio.h>#include <stdlib.h>void wypisz(float **mac, int n){										/*funkcja wypisujaca macierz*/  int i,j;	for (i=0; i<n; i++){		for (j=0; j<=n; j++){			printf("%f ",mac[i][j]);		}		printf("\n");	}}void wyzeruj(float *x){												 	int a=1;	if ((*x)<0) a=-1;	if ((a**x)<10e-6 && (a**x)>0) *x=0;}void schodkuj(float **mac, int i, int j, int size){					 /*funkcja rekurencyjnie sprowadza macierz do postaci schodkowej*/	int wie, m,n;	float tmp;	if (i<size){														/*warunek konca rekurencji*/		wie=i;		while (wie<size-1&&mac[wie][j]==0)wie++;						/*szukamy w danej kolumnie elementu roznego od zera i jesli taki istnieje to przechodzimy dalej przez IF'a*/		if (mac[wie][j]){			if (wie!=i){												/*jesli znaleziony wyzej wiersz z niezerowym elementem nie jest wierszem i-tym to zamieniamy te wiersze*/				for (m=j; m<=size; m++){				tmp=mac[wie][m];				mac[wie][m]=mac[i][m];				mac[i][m]=tmp;				}			}			for (m=i+1; m<size; m++){				tmp=mac[m][j]/mac[i][j];				wyzeruj(&tmp);				for (n=j; n<=size; n++){								/*Odejmujemy od m-tego wiersza macierzy tmp krotnosc wiersza i-tego.*/					mac[m][n]-=(mac[i][n])*tmp;					wyzeruj(&mac[m][n]);				}			}		}	   schodkuj(mac, i+1, j+1, size);								   /*wywolujemy rekurencyjnie ta funkcje dla macierzy pomniejszonej o pierwsza kolumne i pierwszy wiersz*/	}}float rozwiaz(float **mac, int i, int size){							/*Ta funkcja wyliczy rozwiazania x1 x2 itd*/	float a=0;	int m;	for (m=i+1; m<size; m++)		a+=(mac[i][m]*rozwiaz(mac,m,size));	return ((mac[i][size]-a)/mac[i][i]);}/*Funkcje odpowiadajace za sprawdzanie czy uklad posiada rozwiazania*/int czy_rowny(float **mac, int i, int size){							/*Jest to podfunkcja funkcji "sprawdz", wywolana w odpowiednim momencie pozwala wykryc sprzecznosc ukladu, kiedy i po co ja wywolujemy opisane zostalo nizej*/	int m,n,flaga;	for (m=0; m<size; m++){		if(m!=i){			flaga=1;			for (n=0; n<size; n++)				if (!(mac[m][n]==mac[i][n])) flaga=0;			if (flaga) return 1;		}	}	return 0;}/*DOTYCZY POSTACI SCHODKOWEJ!*/int sprawdz (float **mac, int size){									/*Funkcja zwraca: 0 gdy uklad jest sprzeczny, 1 gdy ma nieskonczenei wiele rozwiazan, 2 gdy istnieje jedno rozwiazanie*/	int m,n,accident=0,licznik=0;									   /*Ta funkcja jak i wszystkie jej wyrazanie i komentarze odnosza sie do POSTACI SCHODKOWEJ MACIERZY!!*/	float a;	a=(float)0/(float)2;	for (m=0; m<size; m++){		if ((mac[m][m])==0){			accident=1;												 /*jezeli wspolczynnik na przekatnej glownej macierzy jest rowny 0, to oznacza to sprzecznosc lub nieskonczona ilosc rozwiazan ukldu, co za chwile ustalimy*/			if (mac[m][size]!=0){									   /*jezeli wyraz wolny tego wiersza jest rozny od 0 to moze zajsc sprzecznosc*/				licznik=0;				for (n=m+1; n<size;n++)					if ((mac[m][n])!=0) licznik++;				if (!licznik) return 0;								 /*jezeli wszystkie wspolczynniki w tym wierszu sa rowne 0 a wyraz wolny rozny od 0 to uklad jest sprzeczny, dlatego tak jak wczesniej ustalilismy zwracamy 0*/				else if(czy_rowny(mac, m, size))return 0;}			  /*jezeli w tym wierszu byl co najmniej 1 niezerowy wspolczynnik to uklad moze byc sprzeczny jezeli istnieje inny wiersz o tych samych wspolczynnikach lecz innym wyrazie wolnym, co sprawdza funkcja czy_rowny*/		}	}	if (accident)return 1;											  /*jezeli ktorys element na przekatnej byl rowny 0 a nie wykrylismy wczesniej sprzecznosci to uklad bedzie posiadal nieskonczona ilosc rozwiazan i tak jak wczeniej ustalilismy zwracamy 1*/	return 2;														   /*jezeli wszystkie wspolczynniki na przekatnej byly rozne od 0 to istnieje jedno rozwiazanie, zwracamy 2*/}/*Koniec funkcji sprawdzajacych czy uklad posiada rozwiazania*/int main(int argc, char *argv[]){	FILE *plik;	char *nazwa=argv[1];	float **mac;	int a,n,i,j,rozwiazanie;	/*----------------PRACA Z PLIKIEM----------------------------*/	if (!(plik = fopen (nazwa,"r"))){								   /*Otwieramy plik jesli to mozliwe*/		printf("Nie mozna otworzyc pliku!");		exit (0);	}	if(!fscanf(plik,"%d",&n)||n<0){									 /*Sczytujemy wielkosc macierzy o ile jest poprawna*/		printf("Niepoprawny rozmiar macierzy!");		exit (0);	}	if (!(mac=(float**)(malloc(n*(sizeof(float*)))))){				  /*alokujemy dynamicznie tablice o ile starcza pamieci*/		printf("Za malo pamieci!");		exit (0);	}	for (i=0; i<n; i++){		if(!(mac[i]=(float*)(malloc((n+1)*(sizeof(float)))))){			printf("Za malo pamieci!");			exit (0);		}	}	for (i=0; i<n; i++){												/*wczytujemy dane do tablicy o ile dane w pliku sa poprawne*/		for (j=0; j<=n; j++){			if(((a=fscanf(plik,"%f",&mac[i][j]))==EOF)||!a){				printf("Niepoprawne dane w pliku!");				exit (0);			}		}	}	fclose(plik);	/*--------------KONIEC PRACY Z PLIKIEM---------------------*/	printf("Oto nasz uklad: \n\n");	wypisz(mac,n);	schodkuj(mac,0,0,n);	rozwiazanie=sprawdz(mac,n);	if (!rozwiazanie)printf("\nUklad sprzeczny!");	else if(rozwiazanie==1) printf("\nUklad posiada nieskonczona ilosc rozwiazan!");		 else {			 printf("\nUklad posiada dokladnie 1 rozwiazanie: \n\n");   /*Jesli sa rozwiazania to je wyznaczamy*/			 for(i=0; i<n; i++)				 printf("X[%d] = %f\n",i+1,rozwiaz(mac, i, n));		 }	for (i=0; i<n; i++) free(mac[i]);								   /*Zwalniamy pamiec*/	free(mac);	return 0;}

Edytowane 2 Marca 2009 przez krawetko

[Boczek11]

Dziękuje bardzo za pomoc.

[Boczek11]

Próbowałem przedstawić metodę eliminacji Gaussa w mathematice, lecz coś mi się "psuje" tzn. wynik dla x3 jest poprawny natomiast błędne wyniki wychodzą w przypadku x1 jaki i x2. Czy ktoś mógłby mi pomóc??

 

 

 

Gauss[a1_, b1_] := Module[{n, a, b, p, x}, a = a1;	b = b1;	n = Length[b1];	x = Table[0, {i, 1, n}];	Do[Do[p = a[[i, j]]/a[[j, j]]; b[[i]] = b[[i]] - p*b[[j]]; 	a[[i]] = a[[i]] - p*a[[j]], {i, j + 1, n}], {j, 1, n - 1}]; 	x[[n]] = b[[n]]/a[[n, n]];	Do[x[[n - i]] = (b[[n - i]] - 		Do[x[[n - k]]*a[[n - i, n - k]]{k, 0, i - 1, 1}])/a[[n - 		i, n - i]], {i, 1, n - 1}];	Return[x];	]Gauss[{{1/4, 1/5, 1/6}, {1/3, 1/4, 1/5}, {1/2, 1, 2}}, {9, 8, 8}]

Edytowane 3 Marca 2009 przez Boczek11