C Gauß-Elimination - Unterscheidung keine/unendlich viele Lösungen

[simpsonsfan]

Servus zusammen,

ich hab mir für die Lösung eines LGS Ax=b das Gauß-Verfahren in C geschrieben:

Spoiler: Code

#include <stdio.h>
#include <stdlib.h>
#include <math.h> //für fabs()

void swap_lines(double **matrix, int line1, int line2){
	double *temp;
	temp=matrix[line1];
	matrix[line1]=matrix[line2];
	matrix[line2]=temp;
}

void add_lines(double **matrix, int changeline, int addline, double factor_k, int dimension){
	// Adds factor_k*addline to changeline
	int j;
	for(j=0;j<dimension+1;j++){
		matrix[changeline][j] += factor_k*matrix[addline][j]; 
	}
}

void pivotisierung(double **matrix, int dimension){
	// Hinweis: Da nur eindeutig lösbare Systeme berechnet werden sollen, fallen Systeme mit Nullspalten automatisch aus der Betrachtung, womit sich die Pivotisierung hier auf Zeilenvertauschungen beschränken kann (keine Spaltenvertauschungen)
	
	int i, tauschzeile;
	double hoechster_wert=0;
	for(i=1;i<dimension;i++){
		if ((fabs(matrix[i][0])) > hoechster_wert){
			hoechster_wert=fabs(matrix[i][0]);
			tauschzeile=i;
		}
	}
	swap_lines(matrix, 0, tauschzeile);
}

#define EINDEUTIG 0
#define UNENDLICH 1
#define KEINE 2
int gauss_elim(double ***LGS, int dimension){
	int i,j, status;
	double *temp, k;
	
	if (dimension==1) return EINDEUTIG;
	temp=malloc(dimension*sizeof(double));
	//Prüfung auf Nullspalten
	for(j=0;j<dimension;j++){
		temp[j]=0;
		for(i=0;i<dimension;i++){
			temp[j]+=fabs((*LGS)[i][j]);
		}
		if (temp[j]==0) return UNENDLICH; //FALSCH! Nullspalten führen zu keinem eindeutigen Ergebnis, es kann allerdings immer noch vorkommen, dass keine Lösung exisitiert.
	}
	//Prüfung auf Nullzeilen
	for(i=0;i<dimension;i++){
		temp[i]=0;
		for(j=0;j<dimension;j++){
			temp[i]+=fabs((*LGS)[i][j]);
		}
		if (temp[i]==0){
			if ((*LGS)[i][dimension]==0){
				return UNENDLICH;
			}else{
				return KEINE;
			}
		}
	}
	
	//
	if ((*LGS)[0][0]==0) pivotisierung((*LGS), dimension);
	
	//vorwärtselimination
	for(i=1;i<dimension;i++){
		k= -((*LGS)[i][0])/((*LGS)[0][0]);
		add_lines((*LGS), i, 0, k, dimension);
	}
	
	//Pointer in Diagonale weiter setzen
	for(i=0;i<dimension;i++){
		((*LGS)[i])++;
	}
	(*LGS)++;
	status=gauss_elim(LGS, dimension-1);
	
	//Pointer rücksetzen auf Anfang
	(*LGS)--;
	for(i=0;i<dimension;i++){
		((*LGS)[i])--;
	}
	return status;
}

int gauss_einsetzen(double **LGS, double *loesung, int dimension){
	int i, j;
	for(i=dimension-1;i>=0;i--){
		loesung[i]=LGS[i][dimension];
		for(j=dimension-1;j>i;j--){
			loesung[i]-=(LGS[i][j])*(loesung[j]); //LGS[n][dimension]=b[n]
		}
		loesung[i]/=(LGS[i][i]);
	}
}

int main(void){
	char dateiname[80], temp;
	FILE *datei, *ausgabedatei;
	int anzahl_zeilen, i,j, c;
	double **input, *loesung;
	
	//datei oeffnen
	printf("Bitte Dateinamen eingeben: ");
	scanf("%s", dateiname);
	datei=fopen(dateiname, "r");
	if (datei){
		//zeilenzahl ermitteln
		anzahl_zeilen=1;
		while((temp=getc(datei))!=EOF){
			if (temp=='\n') anzahl_zeilen++;
		}
	}else{
		printf("Die Datei wurde nicht gefunden.");
		return EXIT_FAILURE;
	}
	
	input=malloc(anzahl_zeilen*sizeof(double*));
	for(i=0;i<anzahl_zeilen;i++){
		input[i]=malloc((anzahl_zeilen+1)*sizeof(double));
	}

	//Daten einlesen
	fseek(datei, 0L, SEEK_SET); //Dateizeiger auf Anfang setzen
	for(i=0;i<anzahl_zeilen;i++){
		for(j=0;j<anzahl_zeilen+1;j++){
			fscanf(datei, "%le", &(input[i][j]));
		}
	}
	fclose(datei);
	
	switch (gauss_elim(&input, anzahl_zeilen)){
		case EINDEUTIG :
			loesung=malloc(anzahl_zeilen*sizeof(double));
			gauss_einsetzen(input, loesung, anzahl_zeilen);
			for(i=0;i<anzahl_zeilen;i++){
				printf("%le\n", loesung[i]);
			}
		break;
		case KEINE : printf("null");
		break;
		case UNENDLICH : printf("unendlich");
		break;
		default : printf("Dieser Fall kann nicht auftreten");//Dieser Fall darf nie erreicht werden;
	}
	
	return EXIT_SUCCESS;
}

Da eine explizite Lösung nur im Falle einer 0-dimensionalen Lösung verlangt ist und sonst nur die Ausgabe "unendlich viele Lösungen" bzw. "keine Lösungen",
dachte ich, ich kann es mir einfach machen und diese Fälle direkt aussieben. Nun stellt sich allerdings das Problem dar, dass bei Vorhandensein einer Nullspalte zwar eine eindeutige Lösung auszuschließen ist, allerdings ist es doch noch möglich, dass keine Lösung existiert. Die Zuweisung in Zeile 49 ist also nicht korrekt.

Nun meine Frage, welchen Weg seht ihr, um hier eine Unterscheidung zu treffen?
Muss ich doch die Vorwärtselimination in diesem Fall bis zum Schluss durchführen und dann prüfen, ob keine oder unendlich viele Lösungen exisitieren? Was ja wiederum bedeuten würde, dass ich bei der Pivotisierung auch Spaltenvertauschungen durchführen müsste.
Da ich allerdings auf die explizite Lösung im 1-oder-merhdimensionalen Fall verzichten kann, müsste ich mir wenigstens die Spaltenvertauschungen nicht merken.
Trotzdem würde ich, wenn möglich, einen einfacheren Weg bevorzugen, momentan komme ich aber nicht drauf, wie man das alternativ vorher schon sehen könnte.

Vielleicht habt ihr ja eine Idee?

 

[convexus]

Die Entscheidung im singulären Fall kannst du anhand der entsprechende Komponente des Bildvektors machen. Steht im [transformierten] Bildvektor eine Null gibt es unendliche viele Lösungen ansonsten keine.

 

[simpsonsfan]

Dazu muss ich dann wohl aber erst die Vorwärtselimination bis zum Ende durchführen. Womit ich bei einer Nullspalte am Anfang aber auch ein Pivot aus einer anderen Spalte brauche.
Na gut, dann wird das halt noch mit eingebaut... Es sei denn, jemand wüsste doch noch eine Vereinfachung.