Lineare Gleichungssysteme lösen -> Problem

[Icarus]

Hallo,

ich versuche mich derzeit an einer eigenen Klasse zum lösen von linearen Gleichungssystemen der Form Ax=y, wobei A eine quadratische Matrix mit einer beliebigen Zeilen-/Spaltenanzahl n ist.

Mein Code führt auch zu den richtigen Ergebnissen, solange an den entscheidenden Stellen keine Nullen im Spiel sind.

(Anmerkung: Die Spalte ganz rechts repräsentiert den y Vektor)

(Konsolenausgabe) Matrix, die korrekt gelöst wird:

( 2.0 | 3.0 | 0.0 )
( 4.0 | 7.0 | 1.0 )

Ergebnis: [-1.5, 1.0]

(Konsolenausgabe) Matrix, bei der Probleme auftauchen:

( 0.0 | 2.0 | 1.0 )
( 2.0 | 1.0 | 1.0 )

Ergebnis: [NaN, 0.5]


Wie dieser Fehler zustande kommt ist mir klar: Es kommt zu einer Division von 0 / R (R: eine Reele Zahl), die bei einer "von Hand" ausgeführten Rechnung durch die richtige Sortierung der Zeilen verhindert würde.

Nur wie bringt man dem Computer bei, die Zeilen vorher korrekt zu sortieren?


Ich habe mich in dem Zusammenhang ein wenig durch das Web gegooglet und oft etwas von Pivotisierung und LR-Zerlegung gelesen, aber leider nicht wirklich verstanden, worum es dabei geht, bzw. wie es funktioniert. ???:L

Kennt jemmand vieleicht einen guten Link oder eine gute Erklärung?

Hoffe Ihr könnt mir weiterhelfen



Mein Code:

import java.util.*;

public class LGS 
{
	
	public static void main(String[] args) 
	{
//		LGS lTestLGS = new LGS( new double[][]{{2,3,0},{4,7,1}});
//		LGS lTestLGS = new LGS( new double[][]{{3,6,-2,-4},{3,2,1,0},{(double)3/2,5,-5,-9}});
		LGS lTestLGS = new LGS( new double[][]{{0,2,1},{2,1,1}});
		lTestLGS.printToConsole();
		System.out.print( Arrays.toString( lTestLGS.solve() ) );
	}
	
	private double lLGS[][];
	
	public LGS ( double pLGS[][] ) 
	{
		lLGS = pLGS;
		
	}
	
	public double[] solve( )
	{
		for( int i=0; i<lLGS.length-1; i++ )
		{
//			System.out.println("LGS Primary cursor is now in line " + i);
			for( int j=i+1; j<lLGS.length; j++ )
			{
//				System.out.println("LGS Secondary cursor is now in line " + j);
				double lFactor = lLGS[i][i] / lLGS[j][i];
				System.out.println("LGS Factor: " + lLGS[i][i] + " / " + lLGS[j][i] + " = " + lFactor);
				for( int k=i; k<lLGS[0].length; k++ )
				{
					System.out.print("LGS Calc: " + lLGS[j][k] + " * " + lFactor + " - " + lLGS[i][k] + " = " );
					lLGS[j][k] = lLGS[j][k] * lFactor - lLGS[i][k];
					System.out.print(lLGS[j][k] + "\r\n");
					
				}
				
//				this.printToConsole();
			}
		}
		
		double lSolutions[] = new double[lLGS.length];
		
		for( int i=lLGS[0].length-2; i>=0; i-- )
		{
			double lTemp = lLGS[i][lLGS[0].length-1];
			for( int j=lLGS[0].length-2; j>i; j--)
			{
//				System.out.print( lTemp );
				lTemp = lTemp - lLGS[i][j] * lSolutions[j];
//				System.out.print( " - " + lLGS[i][j] + " * " + lSolutions[j] + " = " + lTemp + "\r\n");
//				System.out.println("j " + j );
			}
//			System.out.println("i " + i );
			lSolutions[i] = lTemp / lLGS[i][i];
//			System.out.println( "LSG: " + lTemp + " / " + lLGS[i][i] + " = " + lSolutions[i] );
		}
		
		return lSolutions;
	
	}
	
	public void printToConsole()
	{
		String lOutput = new String();
		for( int i=0; i<lLGS.length; i++ )
		{
			lOutput = lOutput + "( ";
			for( int j=0; j<lLGS[i].length; j++ )
			{
				lOutput = lOutput + lLGS[i][j];
				if( j+1 <lLGS[i].length ) lOutput = lOutput + " | ";
			}
			lOutput = lOutput + " )\r\n";
		}
		System.out.println(lOutput);
	}
}

 

[0x7F800000]

Das gehört eigentlich ins Mathematik-unterforum...

Die zeilen darfst du vertauschen wie du willst, von der Reihenfolge der Gleichungen ändert sich die Lösungsmenge des Gleichungssystems ja nicht. Deswegen darfst du in jeder Spalte nach einem pivot element !=0 suchen, und dieses pivot-element auf die diagonale verschieben.

Es könnte aber auch sein, dass du kein pivotelement !=0 findest, wenn etwa die gesamte Spalte von Anfang an 0 gewesen ist.

Was du dann machst, hängt davon ab, was du machen willst.
Wenn du nur eine Partikuläre Lösung finden willst, dann solltest du besser PA=LR zerlegung nehmen, das ist schneller. Für die Anschauung find ich folgende Animation ziemlich gut.

Kommt es dir auf die partikuläre Lösung nicht so sehr an, und willst du stattdessen den Kern von A haben, dann weiß ich gar nicht, was normalerweise gemacht wird. Ich habe zB folgende methode implementiert, die zusätzlich den Kern berechnet:

Klasse zum abspeichern von simplen double-matrizen:

package fastnumerics;

import numbertheory.ZMath;

public class Matrix implements Cloneable{
	
	public double[][] elements;
	
	public Matrix(int rows, int columns){
		if(rows<=0 || columns<=0){
			throw new IllegalArgumentException("number of rows and columns must be positive");		
		}
		elements=new double[rows][columns];
	}
	
	public Matrix(double[][] _elements){
		elements=_elements;
		
		if(rows()<=0 || columns()<=0){
			throw new IllegalArgumentException("number of rows and columns must be positive");
		}
	}
	
	public Matrix(int rows, int columns, double ... _elements){
		elements=new double[rows][columns];
		for(int r=0; r<rows; r++){
			for(int c=0; c<columns; c++){
				elements[r][c]=_elements[r*columns+c];
			}
		}
	}
	
	public int rows(){
		return elements.length;
	}
	
	public int columns(){
		return elements[0].length;
	}
	
	public Matrix transpose(){
		Matrix result=new Matrix(columns(),rows());
		for(int i=0; i<rows(); i++){
		for(int k=0; k<columns(); k++){
			result.elements[k][i]=elements[i][k];
		}
		}
		return result;
	}
	
	public String toString(){
		StringBuffer buffer=new StringBuffer();
		if(columns()>1){
			// write as matrix
			for(int r=0; r<rows(); r++){
				for(int c=0; c<columns(); c++){
					buffer.append(elements[r][c]+"\t");
				}
				buffer.append("\n");
			}
		}else{
			//write as column vector
			buffer.append("(");
			for(int r=0; r<rows(); r++){
				buffer.append(elements[r][0]+(r==rows()-1?"":"\t"));
			}
			buffer.append(")^T");
		}
		return buffer.toString();
	}
	
	public Matrix clone(){
		Matrix clone=new Matrix(rows(),columns());
		for(int r=0; r<rows(); r++){
		for(int c=0; c<columns(); c++){
			clone.elements[r][c]=elements[r][c];
		}
		}
		return clone;
	}
	
	public Matrix multiply(Matrix other){
		Matrix result=new Matrix(this.rows(),other.columns());
		for(int r=0; r<result.rows(); r++){
		for(int c=0; c<result.columns(); c++){
			for(int k=0; k<this.columns(); k++){
				result.elements[r][c]+=this.elements[r][k]*other.elements[k][c];
			}
		}
		}
		return result;
	}
	
	public Matrix add(Matrix other){
		Matrix result=new Matrix(this.rows(),this.columns());
		for(int r=0; r<result.rows(); r++){
		for(int c=0; c<result.columns(); c++){
			result.elements[r][c]=this.elements[r][c]+other.elements[r][c];
		}
		}
		return result;
	}
	
	public static Matrix getIdentity(int n){
		Matrix result=new Matrix(n,n);
		for(int r=0; r<result.rows(); r++){
		for(int c=0; c<result.columns(); c++){
			result.elements[r][c]=ZMath.kronekerDelta(r, c);
		}
		}
		return result;
	}
}

...und hier das eigentliche Gauß'sche Eliminationsverfahren:

package fastnumerics;
import java.util.*;

public class GaussianElimination {
	
	/**
	 * Solves Ax=b, returns particular solution, writes Kernel-basis into the list
	 * warning: this method completely wrecks the matrices A and b !
	 * 
	 * @param a			Matrix A 	
	 * @param b			multiple column vectors b_1,...,b_r combined to one Matrix b
	 * @param basis 	(output) reference to an empty collection of matrices, that is filled with basis vectors of the kernel
	 * 
	 * @return			particular solutions for Ax_1=b_1, ... , Ax_r=b_r combined to one single matrix x so that Ax=b
	 */
	public static Matrix gaussianEliminationWithSideEffects(Matrix a,Matrix b,Collection<Matrix> basis){
		
		// #{basis vectors}+#pivotElements= dim(Ker)+dim(Img) = dim(Source)=number of columns
		/*  from the one side, the array is filled with positions of pivot elements
		 *  from the other side, the array is filled with indices of columns, that have no pivot element
		 */
		int[] pivotPositions_basisVectors=new int[a.columns()];
		int rang=0;
		
		for(int row=0,column=0; row<a.rows() && column<a.columns(); row++,column++){
			//search for a row with a_i0 != 0, get the maximal
			double bestPivot=0; int bestPivotIndex=row;
			for(int pivotSearch=row; pivotSearch<a.rows(); pivotSearch++){
				if(Math.abs(a.elements[pivotSearch][column])>bestPivot){
					bestPivotIndex=pivotSearch;
					bestPivot=Math.abs(a.elements[pivotSearch][column]);
				}
			}
			if(bestPivot>0){
				//save position of the pivot element
				pivotPositions_basisVectors[row]=column;
				rang++;
				//swap row and bestPivotIndex
				double[] temp=a.elements[row];
				a.elements[row]=a.elements[bestPivotIndex];
				a.elements[bestPivotIndex]=temp;
				//also in b
				temp=b.elements[row];
				b.elements[row]=b.elements[bestPivotIndex];
				b.elements[bestPivotIndex]=temp;
				
				double f;
				//eliminate all elements (except for the current row) in the current column
				for(int eliminating=0; eliminating<a.rows(); eliminating++){
					if(eliminating==row) continue;
					f=-a.elements[eliminating][column]/a.elements[row][column];
					a.elements[eliminating][column]=0;
					for(int h=column+1; h<a.columns(); h++){
						a.elements[eliminating][h]+=a.elements[row][h]*f;
					}
					//same for b
					for(int h=0; h<b.columns(); h++){
						b.elements[eliminating][h]+=b.elements[row][h]*f;
					}
				}
				//normalize row
				f=1/a.elements[row][column];
				a.elements[row][column]=1;
				for(int h=column+1; h<a.columns(); h++){
					a.elements[row][h]*=f;
				}
				//same for b
				for(int h=0; h<b.columns(); h++){
					b.elements[row][h]*=f;
				}
			}else{
				//pivot != 0 not found
				//save the column, which later will be transformed into vector of the basis
				pivotPositions_basisVectors[a.columns()-1-column+rang]=column;
				//look at the next column, stay in the same row
				row--;
				continue;
			}
		}
		
		//create basis vectors
		for(int i=rang; i<a.columns(); i++){
			
			Matrix basisVector=new Matrix(a.columns(),1);
			int vectorIndex=pivotPositions_basisVectors[i];
			
			basisVector.elements[vectorIndex][0]=-1;
			for(int k=0; k<rang; k++){
				basisVector.elements[pivotPositions_basisVectors[k]][0]=a.elements[k][vectorIndex];
			}
			
			basis.add(basisVector);
		}
		
		//create particular solution
		Matrix particularSolution=new Matrix(a.columns(),b.columns());
		for(int bRow=0; bRow<rang; bRow++){
			particularSolution.elements[pivotPositions_basisVectors[bRow]]=b.elements[bRow];
		}
		
		//check that rang(b)<=rang(a)
		for(int bRow=rang; bRow<b.rows(); bRow++){
		for(int h=0; h<b.columns(); h++){
			if(b.elements[bRow][h]!=0){
				return null; //no solution
			}
		}
		}
		
		return particularSolution;
	}
	
	/**
	 * Does the same as gaussianEliminationWithSideEffects, but works with clones of A and b.
	 * @param a			Matrix A
	 * @param b			multiple column vectors combined to one Matrix
	 * @param basis		(output) empty collection of matrices, that is filled with basis vectors of the kernel of A
	 * @return			particular solution x so that Ax=b
	 */
	
	public static Matrix gaussianElimination(Matrix a, Matrix b, Collection<Matrix> basis){
		return gaussianEliminationWithSideEffects(a.clone(),b.clone(),basis);
	}
	
	public static void main(String[] args){
		
		//* 
		Matrix a=new Matrix(3,3, 2,2,0, 2,4,1, 0,2,1);
		Matrix b=new Matrix(3,1, 4,7,3);
		//*/
		
		Collection<Matrix> basis=new LinkedList<Matrix>();
		
		System.out.println(a);
		System.out.println(b);
		
		System.out.println("particular="+gaussianEliminationWithSideEffects(a,b,basis));
		System.out.println("basis=");
		for(Matrix m:basis){
			System.out.println(m);
		}
	}
}

Hier werden die Spalten nicht hinundherpermutiert, es gibt auch keine Permutationsmatrix P (gut, die gibt es bei der anderen implementierung eigentlich auch nicht)
Stattdessen werden die positionen der pivot-elemente und der Spalten, wo keine pivot-elemente zu finden waren, in einem array abgespeichert (in einem wohlbemerkt: weil die größe des Arrays vorher eindeutig durch die dimension des Urbildes bestimmt ist. Für mathematische Anschauung ist das etwas besser, zudem spart es speicherplatz, aber es ist vielleicht nicht so gut zu lesen)
Am ende wird aus den spalten mit pivot-elementen die spezielle Lösung zusammengebastelt.
Aus den Spalten ohne pivot-elemente wird der Kern zusammengebaut.

Es funzt ganz gut. Aber alleine aus dem code wirst du das wohl auch nicht nachvollziehen können, rechne lieber selbst ein wenig herum, bis du dir das alles klargemacht hast.


edit: Achso, zum testen solltest du noch das ZMath package rauswerfen, brauchst du eh nicht. kroneker-delta wirst du wohl selbst hinbekommen oder schmeiß die methode, die das benötigt auch raus: für die gauß elimination braucht man das alles nicht.

 

[Icarus]

Hallo Andrey,

entschuldige bitte, dass ich das falsche Forum gewählt habe. Ich dachte, das Mathematikforum sei für mathematische Fragen, welche nichts mit Programmierproblemen zu tun haben.

Danke für den Link und Deinen Code.

Werde beides morgen mal durcharbeiten in der Hoffnung, dass ich dann mein Problem gelöst bekomme! :###

Ansonsten wird man nochmal von mir hören :wink:

 

[0x7F800000]

Ich dachte, das Mathematikforum sei für mathematische Fragen, welche nichts mit Programmierproblemen zu tun haben.

Nee, dafür ist das eher ungeeignet. kein Latex, kein MathML... :?
Ist im Prinzip nicht soo wichtig, solangs nicht bei EE o.ä. landet