Beispiel: Torsionsproblem, FEM-Lösung der Poissonschen Differenzialgleichung
Zur Berechnung des Torsionsträgheitsmoments für einen Rechteckquerschnitt (Aufgabe aus dem Kapitel "Methode der finiten Elemente") wird ein 6-Knoten-Dreieckselement aus Matlab-Femset verwendet (nebenstehend links ein FEM-Netz mit 64 Elementen). Das Matlab-Script, das dort verwendet wird , muss leicht modifiziert werden, damit das entstehende Gleichungssystem untersucht werden kann (dieses Beispiel wird auch auf der Seite "Testrechnungen mit iterativen Verfahren (2)" untersucht).
Nachfolgend ist das Matlab-Script zu sehen, mit dem die Rechnungen ausgeführt werden. Das Gleichungssystem wird mit einer symmetrischen und positiv definiten Koeffizientenmatrix als "Sparse matrix" bereitgestellt (rechts ist das Belegungsmuster der Matrix für das links oben zu sehende FEM-Modell zu sehen).
Der Algorithmus zum Aufbau des Gleichungssystems interessiert für die hier angestellten Betrachtungen nicht, allerdings beeinflusst die Feinheit der Finite-Elemente-Diskretisierung die Größe des Gleichungssystems. Dies wird im Script durch die beiden Werte in Zeile 5 festgelegt.
- In Zeile 19 wird das Gleichungssystem mit dem Standard-Solver von Matlab (Backslash-Operator) gelöst.
- Die alternative Möglichkeit, mit der auch die positive Definitheit der Koeffizientenmatrix überprüft wird, muss in zwei Schritten realisiert werden: Zunächst wird die Cholesky-Zerlegung mit der Function chol ausgeführt (Zeile 23), anschließend wird das Vorwärts-Rückwärteinsetzen mit der zweifachen Anwendung des Backslash-Operators realisiert (Zeile 24).
Verglichen werden die erforderlichen Rechenzeiten, wobei die Relationen doch überraschend sind: Nach den Angaben im Matlab-Helpsystem prüft der Backslash-Operator, ob die Koeffizientenmatrix symmetrisch und positiv definit ist, um bei positiver Antwort dann das besonders effektive Cholesky-Verfahren zu verwenden. Auf diesem Weg könnte der Backslash-Operator allerdings nicht wesentlich schneller zur Lösung kommen als der direkte Aufruf des Cholesky-Verfahrens. Genau das aber weist die folgende Tabelle aus.
Der Backslash-Operator arbeitet offensichtlich nach einem in jeder Hinsicht optimierten Algorithmus. Dazu gehört auch ein "Fill-in-Optimierer", ohne den die Function chol nicht konkurrenzfähig ist. Deshalb wurde in den Zeilen 28 bis 32 die Cholesky-Zerlegung noch einmal mit vorgeschaltetem "Fill-in-Optimierer" und nachgeschaltetem Rücksortierer programmiert (die Beschreibung dieser Strategie findet man auf der Seite "Reduzierung von Bandweite und Fill-in"). In der Tabelle sieht man, dass damit die Function chol auch in den Rechenzeiten konkurrenzfähig wird (die geringfügig größeren Zeiten sind dem Overhead für den Aufruf mehrerer Functions zuzuschreiben, wobei der wesentliche Anteil dem zweifachen Aufruf des Backslash-Operators für das Vorwärts- und Rückwärtseinsetzen anzulasten ist, denn dieser muss jeweils erst herausfinden dass eine Dreicksmatrix ankommt).
|
Anmerkung: Wenn z. B. bei einem FEM-Gleichungssystem die Koeffizientenmatrix des Gleichungssystems den Test auf positive Definitheit besteht, heißt das natürlich nicht, dass sie richtig ist. Aber die typischen Fehler im Berechnungsmodell führen sehr oft auf nicht positive Matrizen, so dass mit dem positiven Ausgang des Test zwar nur ein wertvolles Indiz für die Richtigkeit, bei negativem Ausgang allerdings ein sicherer Beweis für die Fehlerhaftigkeit geliefert wird.
Download
Das oben abgebildete Matlab-Script CholTest.m, die von ihm aufgerufenen Functions TorsnetD6.m, SymBand2Sparse.m und die DLL syswbc_m.dll sind zum Download verfügbar.
