Gleichungssystem lösen...

Goldenflash

Hallo,
dies ist mein erster Eintrag in diesem Forum und ich hoffe mal, dass ich hier richtig bin.
Also mein Problem sieht folgendermaßen aus: Ich habe eine sehr große Matrix mit sehr vielen Nulleinträgen und zwei Vektoren.

http://free.pages.at/osdevdevu/bsp.jpg

Um ein Gleichungsystem dieser Art zu lösen würde ich gerne lapack oder imsl verwenden, allerdings habe ich leider gar keine Ahnung wie ich sowas in ein C++ Programm einbindet . Gibt es irgendwo ein Tutorial oder ein Beispiel wo das ganze für kleine Matrizen gemacht wurde? Ich hab ein paar Seiten auf englisch gefunden die ich aber leider nicht verstanden hab.

Danke schonmal

Mario Sandler

Doch nicht...

Jan

Dein Gleichungssystem scheint mir ein bisschen komisch zu sein. Die 3 Zeilen mit a, b und c sind doch völlig überflüssig?

Goldenflash

Da hast du natürlich recht!

Hab das ganze mal verbessert!
Also a,b,c und x&y sollen natürlich unbekannt sein.

Ben04

Jetzt mal unabhängig von der Frage wie man die von dir angesprochenen Bibliotheken einbindet. Das System hat doch gar keine Lösung. Die erste, dritte und letzte Gleichungen sind unabhängig von der Wahl der Parameter nicht lösbar:
1=0
2=0
3=0
Oder stehe ich gerade auf dem Schlauch?

Michael E.

Nein, das stimmt schon so, wie du das siehst.

=> leere Lösungsmenge

Jan

Goldenflash schrieb:

Hab das ganze mal verbessert!

Verschlimmbessert... Jetzt ist es nicht mehr nur sinnlos, sondern auch noch unsinnig.

Goldenflash

Ja ok, von mir aus hätte ich mir beim Ausdenken eines Beispiels mehr mühe geben sollen.
Is mal wieder verbessert!
Was aber nichts an dem Problem an sich ändert!!!

Ben04

Wieso musst du überhaupt unbedingt auf diese Bibs zurückgreifen? Unter Umständen ist es einfach einen Löser selbst zu schreiben.

Dein aktuelles Beispiel lässt sich zu
23x + 26y=0
8x + 3y=0
vereinfachen.

Goldenflash

Ben04 schrieb:

Wieso musst du überhaupt unbedingt auf diese Bibs zurückgreifen? Unter Umständen ist es einfach einen Löser selbst zu schreiben.

Dein aktuelles Beispiel lässt sich zu
23x + 26y=0
8x + 3y=0
vereinfachen.

Ok, dass mir hier irgendwer sagt wie man etwas so triviales löst hab ich befürchtet! Hatte mir ja auch schon mit meinem Beispiel nicht so viel Mühe gegeben, da es nur ein BEISPIEL ist!

Aber ich werd dir erklären wozu ich es brauche. Ich wollte mir ein kleines FEM Programm schreiben, dabei ist die Matrix meine assemblierte Elementsteifigkeitsmatrix, welche auf die nötigen Zeilen reduziert wurde. Ich multipliziere sie mit dem Kraftvektor und setze nur die bekannten Komponenten des Verschiebungsvektor ein. Standard(in Maple funktioniert es ohne Probleme)! Wenn ich meinen Kraftvektor raus hab kann ich ihn mit der inversen der vollbesetzten Elementsteifigkeitsmatrix multiplizieren und ich hab meinen Verschiebungsvektor mit allen Komponenten.
Wenn du dir jetzt vorstellst, dass mein Problem sagen wir 10000 Freiheitsgrade hat (was nichtmal viel ist) dann heißt das, dass es 20000 Knoten und damit eine 40000*40000 Einträge Matrix ist. Die Vektoren sind entsprechend!

Walli

Und warum machst Du dann so ein dämliches Beispiel und sagst nicht sofort worum es geht? Dann hätte man Dir eher helfen können. Du suchst also eine Möglichkeit dünnbesetzte Matrizen effizient zu lösen. Da deine Matrix symmetrisch und positiv semi-definit ist solltest Du bei großen Systemen einen MINRES-Löser benutzen. Mit LAPACK++, einem Wrapper für die F77-Bib LAPACK, sollte das kein Problem sein den zu implementieren.

Goldenflash

Wir kommen der Sache langsam näher...
An Lapack hatte ich auch gedacht also ich hier geschrieben haben allerdings gibt es ja auch noch andere Löser die ich nicht ausschließen wollte.
Was ich suche ist ein Tutorial wie ich lapack++ in C++ benutze! Am liebsten einfach ein paar Beispiele dann sollte ich das verstehen.

Walli

Was meinst Du mit 'benutzen'? Wie Du die Lib einbindest oder wie Du damit arbeitest? Ich habe selber nie was mit LAPACK gemacht. Ich nutze ganz gerne boost::ublas für solche Sachen. Als C++-Programmierer sollte man sich ohnehin mit boost auseinandersetzen . Beispiele gibt es in der Doku ein paar glaub ich. Ich hab auch mal spaßhalber das CG-Verfahren implementiert. Das solltest Du aber nicht benutzen, da es von strikt positiver Definitheit ausgeht, was Deine Matrix aber nicht ist. MINRES ist aber auch nicht viel schwieriger zu implementieren. Mit nem gescheiten Vorkonditionierer könnte die Implementierung nur noch ein bisschen kniffliger werden.

// conjugate gradient solver
	unsigned cgSolver(SparseMatrix const& A, SparseVector const& f, Vector& x, double& res)
	{
		using namespace boost::numeric::ublas;
		if(A.size1() != A.size2())
			return 0;

		Vector r = f - prod(A, x), p = r;

		for(unsigned nit = 0; ; ++nit)
		{
			Vector Ap = prod(A, p);
			double rTr = inner_prod(r, r);
			double a = rTr / inner_prod(Ap, p);

			x += a * p;
			r -= a * Ap;

			// test for convergence
			double res0 = res;
			res = norm_2(r);
			if(res < eps || res0 >= res)
				return nit;

			double b = inner_prod(r, r) / rTr;
			p = r + b * p;
		}
	}

SparseMatrix usw. sind einfach nur typedefs auf die entsprechenden Typen in uBLAS.