Matrix- und Tensor-Objekte höherer Stufe für numerische Simulationen

Mittels numerischer Simulationen werden in Naturwissenschaft und Technik aus spezifischen Computermodellen wichtige Erkenntnisse über die Gestaltung von natürlichen und technischen Prozessen gewonnen. Dafür sind oft Algorithmen eingesetzt, die auf der linearen Algebra traditioneller Matrixoperationen aufbauen und in entsprechenden Bibliotheken der prozeduralen Programmierung vorzufinden sind. Nach den Prinzipen der objektorientierten Modellierung und Programmierung hat sich dazu inzwischen seit den 1990er Jahren eine Alternative von Matrixklassen herausgebildet, welche auch mehrdimensionale Matrixobjekte erfaßt und die Standardoperatoren von arithmetischen und funktionalen Operationen für entsprechend komplex aufgebaute Objekte bereitstellt. Dazu dient im folgenden eine ausgewählte Hierarchie von Matrix- und Tensorobjekten, deren Struktur und Eigenschaften nach objektorientiertem Prinzip eingekapselt und intern an Unterklassen vererbbar sind. Die Kommunikation mit den Objekten erfolgt dabei ausschließlich über öffentliche Mitgliedsfunktionen, die zusammen mit den überladenen Operatoren der indexbasierten Matrix- und Tensorarithmetik eine intuitiv verständliche Benutzungsoberfläche bilden.

Einleitung

Die Entwicklung der objektorientierten Analyse und Modellierung hat bedeutenden Einfluss auf die Anwendungen der kooperativen Planung und Konstruktion des Ingenieurwesens, bis hin zur gegenwärtigen Produkt- und Bauwerkinformationsmodellierung (BIM). Dabei spielen speziell numerische Modelle für die Simulation von physikalischen Vorgängen in der Festkörper- und Strömungsmechanik eine entscheidende Rolle für das Verständnis komplexer Wechselwirkungen in Natur und Technik. Im Allgemeinen werden die physikalischen Grundlagen mit den Methoden der Tensoranalysis, des Matrixkalküls und der numerischen Programmstrukturen in zuverlässige Algorithmen unterschiedlicher Approximationsmethoden, wie z. B. die Differenzen- oder Finite-Elemente-Methoden, überführt.

Dabei sind Algorithmen in objektorientierten Strukturen von hoher Relevanz, weil durch die zugrunde liegenden Paradigmen der Kapselung und der Vererbung die Transparenz der Formulierung und die Zuverlässigkeit der Methoden abgesichert und die Kommunikation mit den Objekten eindeutig definiert wird. So eröffnet vor allem die Objektorientierung neue Möglichkeiten, die theoretischen Grundlagen der Tensoranalysis mit den numerischen Vorteilen des Matrixkalküls und den Algorithmen der Rechenprogramme zu einer methodischen Einheit zu verknüpfen.

Denn die klassische Matrizenrechnung basiert auf einem Konzept mit einem Schema der zweidimensionalen Anordnung von Rechengrößen, das aus einer Zeit stammt, als es noch keine digitalen Rechner gab. Inzwischen ist es in der Programmierung jedoch Standard, mehrdimensionale Felder (arrays) zu verwenden, um Daten auch dynamisch zu speichern und zu verarbeiten. Um diese Ressourcen optimal zu nutzen ist es notwendig, die Matrizenrechnung auf mehrdimensionale Matrizen auszuweiten und eine geeignete Syntax und Semantik zu formulieren, welche die Matrixarithmetik diesem Konzept anpasst. Dazu bietet sich in idealer Form die bewährte Indexnotation der Tensorrechnung an, wie das folgende Beispiel zeigt

$C_{ikl}=A_{ikl}+B_{ikl}$ bzw. C = A + B

oder

$f^{ilk}=k^{isl}*{v^{k}}_{s}$ bzw. f ≠ k * v$^{T}$

und bei dem die symbolische Matrixnotation leider versagt. Durch das objektorientierte Klassenkonzept hingegen werden geeignete Strukturen für neue Tensor- und Matrixklassen geschaffen, die es ermöglichen, indexbasierte Matrixoperationen für mehrdimensionale Matrizen zu definieren, um die Matrizenrechnung auf Objekte höhere Stufe zu verallgemeinern und die Nomenklatur von Tensoren, Matrizen und Feldern zu vereinheitlichen. Dabei stärkt auch das Überladen der arithmetischen und funktionalen Standardoperatoren für die Matrixobjekte die Benutzerfreundlichkeit der Schnittstellen in ganz erheblichem Maße.

Der vorliegende Beitrag verdeutlicht dieses bewährte Konzept in den nachfolgenden Kapiteln mittels überschaubarer Beispiele der linearen Algebra anhand von übersichtlich gegliederten Prototypen der Programmstrukturen. Damit werden die Abläufe in transparenter Weise nachvollziehbar dargestellt und numerisch verifiziert. Die Methoden der objektorientierten Programmierung sind in ausführlich beschrieben. Eine umfassende Einführung in die Programmiersprache C++ findet man u. a. bei.

Matrixobjekte und überladene Operatoren

Es wird im Folgenden davon ausgegangen, dass Matrizen 2. Stufe spaltenweise gespeichert werden, wie das in naturwissenschaftlich-technischen Anwendungen und algorithmischen Programmiersprachen – wie FORTRAN – üblich ist.

$\mathbf {M} =\left({\begin{array}{ccc}a_{0,0}&a_{1,1}&...&a_{1,N_{k-1}}\\...\\...\\...\\...\\a_{N_{i-1},0}&a_{N_{i-1},1}&...&a_{N_{i-1},N_{k-1}}\end{array}}\right)=M_{ik}$ mit den Dimensionen $(N_{i}\times N_{k})$.

Dabei läuft der erste Index $i$ über die Zeilen und der zweite Index $k$ über die Spalten der Matrix $M_{ik}$.

Die interne Speicherung in Matrixobjekten erfolgt durch Zuweisung der Matrixelemente in sequentieller Reihenfolge zu einem eindimensionalen Vektorfeld:

$\mathbf {V} :=\left({\begin{array}{ccc}a_{0,0}&a_{1,0}&...&a_{N_{i-1},0}&a_{0,1}&a_{1,1}&...&a_{N_{i-1},1}&...&...&a_{0,N_{i-1}}&a_{1,N_{i-1}}&...&a_{N_{i-1},N_{k-1}}\end{array}}\right)$

$=\left({\begin{array}{ccc}v_{0}&v_{1}&...&v_{N_{i-1}}&...&...&v_{N_{i}*N_{k}-1}\end{array}}\right)=V_{z}$ mit der Dimension $N_{z}=N_{i}*N_{k}$,

wie es auch intern auf von Neumann-Rechnerarchitekturen geschieht.

Die entsprechenden bijektiven Abbildungen erfolgen durch die Vorschriften für

die Abbildung

$V_{z}:=M\underbrace {_{ik}} _{z}$ durch die Zuordnung $z=k*N_{i}+i$,

und die Rückabbildung

$M\underbrace {_{ik}} _{z}:=V_{z}$ durch die Zuordnungen $k={\text{Entier}}(z/N_{i})$ und $i=z-k*N_{i}$.

Der Deskriptor des Objektes für die Darstellung der Matrix ist offenbar unabhängig von der internen sequenziellen Speicherung der Elemente und steuert lediglich nach außen die spezielle Darstellung der Matrixform.

Reduziert man die Matrix $\mathbf {M}$ auf eine einzelne Zeile, so entsteht ein Zeilenvektor $M_{0k}$, entsprechendes gilt für einen Spaltenvektor $M_{i0}$. Daraus erkennt man, dass es bei der Indexschreibweise möglich ist, in die Nomenklatur 0-Indizes einzufügen, ohne dadurch die Speicherung von Vektoren oder Matrizen zu ändern.

Dies gilt auch für den Vektor $V_{z}:=V_{0z}:=V_{z0}$ wie auch für beliebige Matrizen: $M_{ik}:=M_{0ik}:=M_{i0k}:=M_{ik0}$.

Über den Konstruktor $M(N_{i},N_{k},N_{l},N_{m})$ wird in objektorientierter Weise die mathematische Grundform auf Matrizen höherer Stufe erweitert und durch die Wahl der Dimensionen festgelegt. Zulässig sind zunächst bis zu vier Dimensionen, die für ingenieurmäßige Anwendungen in der Regel ausreichen, zumal auch mehrere benachbarte Indizes mittels der obigen Abbildungen zusammengefasst werden können, wovon im Programmcode Gebrauch gemacht wird.

Anhand des folgenden C++ Programmbeispiels wird zuerst erläutert, wie die neuen Matrizenklassen semantisch aufgebaut sind, um entsprechende Matrixobjekte (Instanzen) zu erzeugen und diese mit einer indexbasierten Arithmetik mittels überladener Standardoperatoren syntaktisch zu verarbeiten. Mit den Konstruktoren der Klassen Vektor, MATRIX, delta_MATRIX, e_MATRIX können statische oder dynamische Objekte mit Elementen vom Typ double erzeugt werden. Über ein Feld (array) in Notation der Programmiersprache C ist es auch möglich, die Vektor- und Matrixelemente mit Zahlenwerten vorzubelegen. Die Ein-/Ausgabe erfolgt über die überladenen Standardoperatoren >> bzw. <<.

Prototypen neuer Klassen und Methoden

cout << endl
     << "1.1 - Matrizenobjekte und überladene Operatoren" << endl
     << "===============================================" << endl <<endl;
//
// C-Feld mit Werten vorbelegen
// ----------------------------
double F[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
// cout << "C-Vektor F(10) ="; for (int j = 0; j < 10; j++) cout << " " << F[j]; // drucken
cout << endl << endl;
//
// Erzeugen einer statischen Instanz der Klasse VEKTOR mittels Konstruktor
// -----------------------------------------------------------------------
VEKTOR V(F, 6);                                                                  // (1x6)-Vektor V wird mit F vorbelegt
cout << "(1x6)-Vektor V := F = " << V << endl;                                   // überladener Operabisr '<<'
//
// Erzeugen einer statischen Instanz der Klasse MATRIX mittels Konstruktor
// -----------------------------------------------------------------------
MATRIX M(3, 3);                                                                  // (3x3)-matrix, assignment of V bis M
//
// Ausnahme: Typumwandlung (casting) durch M := V zur Demonstration
M = V;                                                                           // überladener Operator '='
cout << "(3x3)-Matrix M := V = " << M << endl;                                   // drucken
//
// Erzeugen einer statischen Instanz der Klasse delta_MATRIX mittels Konstruktor
// -----------------------------------------------------------------------------
delta_MATRIX d(2, 2);                                                            // (2x2)-delta-matrix
cout << "(2x2) delta-Matrix d = " << d << endl;                                  // drucken
//
// Erzeugen einer statischen Instanz der Klasse MATRIX mittels Konstruktor
// -----------------------------------------------------------------------
MATRIX A(2, 2);                                                                  // (2x2)-matrix
cout << "(2x2)-Matrix A = " << A << endl;                                        // drucken

Test von Effektivität und Genauigkeit

C-Vektor F(10) = 0 1 2 3 4 5 6 7 8 9

(1x6) Vektor V := F = [ (1x6) * (1x1) ]-MATRIX:

Spalten 0 bis 5

0.000e+00  1.000e+00  2.000e+00  3.000e+00  4.000e+00  5.000e+00

Pause -> Zeichen eingeben:�

(3x3) Matrix M := V = [ (3x3) * (1x1) ]-MATRIX:

Spalten 0 bis 2

0.000e+00  3.000e+00  0.000e+00
1.000e+00  4.000e+00  0.000e+00
2.000e+00  5.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

(2x2) delta-Matrix d = [ (2x2) * (1x1) ]-MATRIX:

Spalten 0 bis 1

1.000e+00  0.000e+00
0.000e+00  1.000e+00

Pause -> Zeichen eingeben:�

(2x2) Matrix A = [ (2x2) * (1x1) ]-MATRIX:

Spalten 0 bis 1

0.000e+00  0.000e+00
0.000e+00  0.000e+00

Die zu Anfang festgelegten Dimensionen eines Matrixobjektes $\mathbf {M}$ können im vorliegenden Konzept jedoch auch zur sinnvollen Ausführung der überladenen Operatoren temporär abgeändert werden. Dies geschieht durch Aufruf der an das Objekt $M$ gebundenen Mitgliedsfunktion M.dimension(N_i, N_k, N_l, N_m) und gilt lediglich bis zum Ende der jeweiligen Operation!

Die Funktion dimension ist speziell wichtig für die Realisierung der allgemeinen Matrizenmultiplikation – wie noch gezeigt wird.

Das folgende Beispiel zeigt zunächst elementare Anwendungen der arithmetischen Mitgliedsfunktionen der Addition und der Subtraktion mit den überladenen Operatoren =, + und - sowie die Mitgliedsfunktion determinante mit der deklarierten Bindung an die definierten Klasseninstanzen (Objekte).

Prototypen neuer Klassen und Methoden

cout << endl
     << "1.2 - Addition und Subtraktion" << endl
     << "==============================" << endl << endl;
//
// Matrix M wird durch die Mitgliedsfunktion dimension temporaer verkleinert!
//                         --------------------------
// Matrizen-Addition M(i,k) + A(i,k)
// ---------------------------------
A = M.dimension(2,2) + d;                                                        // überladener Operator '+'
cout << "Matrix A = A + d = " << A << endl;                                      // drucken
cout << "Original (3x3)-Matrix preserved: M = " << M << endl;                    // drucken
//
// Matrizen-Subtraktion A(i,k) - d(i,k)
// ------------------------------------
A = A - d;                                                                       // überladener Operator '-'
cout << "Matrix A = A - d = " << A << endl;                                      // drucken
//
// Determinante von Matrix A
// -------------------------
double S;
S = A.determinante(2);
cout << "Determinante det(A) = " << S << endl << endl;

Test von Effektivität und Genauigkeit

Pause -> Zeichen eingeben:�

Matrix A = A + d = [ (2x2) * (1x1) ]-MATRIX:

Spalten 0 bis 1

1.000e+00 2.000e+00
1.000e+00 4.000e+00

Pause -> Zeichen eingeben:�

Original (3x3)-Matrix erhalten: M = [ (3*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

0.000e+00  3.000e+00  0.000e+00
1.000e+00  4.000e+00  0.000e+00
2.000e+00  5.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

Matrix A = A – d = [ (2x2) * (1x1) ]-MATRIX:

Spalten 0 bis 1

0.000e+00 2.000e+00
1.000e+00 3.000e+00

Pause -> Zeichen eingeben:�

Determinante det(A) = -2

Kartesische Basis und Matrix-Vektor-Multiplikationen

Für die Multiplikation von zwei Rechteckmatrizen A und B bestehen i.a. mehrere Varianten, wie A*B, A^T*B, A*B^T, A^T*B^T, die üblicherweise unterschiedliche Prozeduren erfordern. Diese Varianten werden hier in der allgemeinen Matrizenmultiplikation mit der vierdimensionalen Ergebnismatrix

$\mathbf {C} (i,k,l,m)=\mathbf {A} (i,{\textbf {s}},k)*\mathbf {B} (l,\mathbf {s} ,m)$

zusammengefasst, wobei über den Summationsindex s zu summieren ist.

Die entsprechende Anweisung im C++ Programm hat mit den im Konstruktor passend vereinbarten Matrixdimensionen und dem überladenen arithmetischen Operator * die einfache Form

$\mathbf {C} =\mathbf {A} *\mathbf {B} ;$

Für allgemeine Fälle, bei denen die Dimensionen angepasst werden müssen, ergibt sich die Modifikation dieser Anweisung i.a. durch Verwendung der Mitgliedsfunktion dimension wie folgt:

$\mathbf {C.} dimension$$(N_{i},N_{k},N_{l},N_{m})$ = $\mathbf {A.} dimension$$(N_{i},N_{s},N_{k})$ * $\mathbf {B.} dimension$$(N_{l},N_{s},N_{m})$,

wobei deren Dimensionen mit $1$ vorbelegt sind und die Anwendung optional erfolgt.

Für die Verwendung dieses Algorithmus auf Vektoren und Matrizen mit unterschiedlichen Dimensionen ist es dabei wichtig, die folgenden Regeln zu beachten:

- wie bereits erläutert, können benachbarte Indizes mit der Funktion dimension temporär – nur für die Ausführung dieser Operationen – zusammengefasst werden, ohne die Originaldimensionen oder die interne Speicherung der Matrixelemente zu verändern. Damit ist es möglich, beliebige Matrizen diesem Schema anzupassen. Das gilt auch für den Summationsindex;

- für einzelne Indizes kann auch der 0-Index verwendet werden. Was bedeutet, dass an dieser Stelle die Dimension 1 angesetzt und der Matrix zugeordnet wird. Auch diese Maßnahme ändert nicht das interne Speicherschema;

- mit beiden Maßnahmen wird erreicht, dass der Summationsindex stets in die mittlere Position gebracht werden kann;

- wegen der temporären Änderung der Dimensionen mit der Funktion dimension und auch wegen der Abfolge von mehreren Multiplikationen, die der Compiler bestimmt, sollte diese allgemeine Multiplikation nicht verkettet werden und muss in der Regel in einer einzelnen Anweisung ausgeführt werden!

- Vektoren v(N_i) der Klasse VEKTOR sind aus der Klasse MATRIX abgeleitet und werden intern stets als Zeilenvektoren v(1,N_i) definiert. Diese Maßnahme dient der Vereinfachung der Skalarprodukte, weil dadurch der Summationsindex stets auf die mittlere Position voreingestellt ist – was die Handhabung optimiert. Auch hier ist die Anwendung der Funktion dimension unbenommen.

Die Basisvektoren eines kartesischen Bezugssystems (i.a. als x-y-z-System bezeichnet)

${\vec {e}}{^{i}}={\vec {e}}{_{i}}=({\vec {e}}{_{0}}\;{\vec {e}}{_{1}}\;{\vec {e}}{_{2}})$ mit den Komponenten :${e_{0}}^{i}=[1,0,0]),\;{e_{1}}^{i}=[0,1,0]$, und ${e_{2}}_{i}=[(0,0,1]\;$

sind im Folgenden mit lateinischen Indizes gekennzeichnet.

In diesem räumlichen Bezugssystem werden im Beispiel zwei physikalische Vektoren wie folgt angeordnet:

${v^{0}}_{i}=(6,0,0)$ und ${v^{1}}_{k}=(0,2,0)\;.$

Das Skalarprodukt der Vektoren ergibt sich damit zu

$s={v^{0}}_{i}*{v^{1}}_{i}\;,$

während das dyadische Produkt aus der Multiplikation

$M_{ik}={v^{0}}_{i}*{v^{1}}_{k}$

folgt.

Das Levi-Civita-Symbol – im kartesischen System auch als e-Matrix bezeichnet – ist als Permutationssymbol wie folgt durch die

schiefsymmetrische $(2\times 2)$-Matrix bzw. $(3\times 3\times 3)$-Matrix mit den Zahlenwerten ±1:

wie folgt definiert:

$e_{ik}=\left({\begin{array}{lc}0&1\\-1&0\\\end{array}}\right)$ bzw. $e_{ik0}=\left({\begin{array}{lcr}0&0&0\\0&0&1\\0&-1&0\\\end{array}}\right)$ $e_{ik1}=\left({\begin{array}{lcr}0&0&-1\\0&0&0\\1&0&0\\\end{array}}\right)$ $e_{ik2}=\left({\begin{array}{lcr}0&1&0\\-1&0&0\\0&0&0\\\end{array}}\right).$

Mit Hilfe der e-Matrix lässt sich im kartesischen System das Kreuzprodukt der Vektoren ${e_{0}}^{k}$ und ${e_{1}}^{l}$ wie folgt ermitteln:

$e{_{2}}_{i}=e_{ikl}*e{_{0}}^{k}*e{_{1}}^{l}\;.$

\newline

Als Beispiele von typischen Matrizenmultiplikationen werden folgende Operationen mit $(2\times 2)$-Rechteckmatrizen ausgewählt:

$C_{ik}=A_{is}*B_{sk}$ and $C_{ik}=A_{si}*B_{sk}.$

Ein weiteres Beispiel betrifft die Graßmann identity, cf. (s. auch Tensorkalkül)

${\delta _{ik}}^{lm}=e_{iks}*e^{lms}\;,$

welche auf die delta-Matrix 4.Stufe führt.

Der folgende Programmcode erläutert diese Beispiele anhand von Multiplikationen mit entsprechenden VEKTOR- und MATRIX-Objekten.

Prototypen neuer Klassen und Methoden

cout << endl
     << "2.1 - " Cartesische Basis und Matrix-Vektor-Multiplikations"<< endl
     << "===========================================================" << endl << end;
//
// => Vektoren v(i) werden intern als Matrizen v(0,i) mit N(0) = 1 behandelt,
//    um Skalarprodukte einfach bilden zu können!
// -----------------------------------------------------------------
//
// Erzeugen von drei kartesischen Bezugsvektoren
VEKTOR e0(3), e1(3), e2(3);                                                      // (1x3)-Vektoren
e0(0) = 1.0;  e1(1) = 1.0;
//
std::cout << "(1x3) Vektor e0 = " << e0 << endl;                                 // drucken
std::cout << "(1x3) Vektor e1 = " << e1 << endl;                                 // drucken
//
// Erzeugen von zwei physikalischen Vektoren
VEKTOR v0(3), v1(3);                                                             // (1x3)-Vektors
v0(0) = 6.0; v0(1) = 2.0;
std::cout << "(1x3) Vektor v0 = " << v0 << endl;                                 // drucken
//
v1(0) = 1.0; v1(1) = 2.0;
std::cout << "(1x3) Vektor v1 = " << v1 << endl;                                 // drucken
//
// Berechnung des Skalarproduktes s = v0(0,s) * v1(0,s)
// ----------------------------------------------------
double s;
s = v0 * v1;
//
std::cout << "Skalar s = v0 * v1 = " << s << endl <<                             // drucken
//
// Berechnung des dyadischen Produktes  M(i,k) = v0(i) * v1(k)
// -----------------------------------------------------------
M = v0.dimension(3) * v1.dimension(3);
//
std::cout << "Dyadisches Produkt M = v0 * v1 = " << M << endl;                   // drucken
//
// Berechnung des Vektorproduktes e2(0,i) = e(ikl) * e0(0,k) * e1(0,l)
// -------------------------------------------------------------------
// Erzeugen einer statischen Instanz der Klasse e_MATRIX mittels Konstruktor
e_MATRIX e(3,3,3);                  // (3x3x3)-e-Matrix
//
std::cout << "(3x3x3) e-Matrix e = " << e << endl;                               // drucken
//
// verkettete Multiplikation ohne die Mitgliedsfunktion dimension
e2 = (e * e0) * e1;
//
std::cout << "(1x3) Vektor e2 = e0 x e1 = " << e2 << endl;                       // drucken
//
// allgemeines Matrizen-Produkt C(i,k,l,m) = A(i,s,k) * B(l,s,m)
// -------------------------------------------------------------
//    Mittels der Mitgliedsfunktion M.dimension(N{i}, N{k}, N{l}, N{m})
// => kann Matrix M(s,k) um einen 0-Index erweitert werden
//    M(s,k) =: M(0,s,k) mit N{0} = 1,
// => können Matrizenindizes M(i,k,l) zu M(IK,s,l) zusammengefasst werden
//    M(i,k,s,l) =: M(IK,s,l)  mit N{IK} = N{i} * N{k},
// => um den Summationsindex s in die mittlere Position M(i,s,l) zu bringen !
// --------------------------------------------------------------------
std::cout << "(2x2) Matrix A = " << A << endl;                                  // drucken
//
// Erzeugen der Matrizen C and B
MATRIX C(2,2), B(2,2);                                                          // (2x2) Matrizen
B = d;   B(0,1) = 0.5;
//
std::cout << "(2x2) Matrix B = " << B << endl;                                  // drucken
//
// Matrixmultiplikation C(i,k) = A(i,s) * B(0,s,k)
// -----------------------------------------------
C = A * B.dimension(1, 2, 2);
//
std::cout << "(2x2) Matrix C = A * B  = " << C << endl;                         // drucken
//
// Matrixmultiplikation C(i,k) = A(0,s,i) * B(0,s,k)
// -------------------------------------------------
C = A.dimension(1,2,2) * B.dimension(1, 2, 2);
//
std::cout << "(2x2) Matrix C = A(T) * B  = " << C << endl;                      // drucken
//
// Grassmann Identity
// ------------------
// delta-Matrix(i,k,l,m) = e-Matrix(i,k,s) * e-Matrix(l,m,s)
// Erzeugen einer statischen Instanz der Klasse MATRIX mittels Konstruktor
MATRIX D4(3,3,3,3);                 // (3x3x3x3) Matrix
//
D4 = e.dimension(3 * 3, 3) * e.dimension(3 * 3, 3);
//
std::cout << "(3x3x3x3) delta-Matrix 4. Stufe D4(i,k,l,m) = " << D4 << endl;     // drucken
//
// Probe
// ----
// Erzeugen einer dynamischen Instanz der Klasse delta_MATRIX mittels des Konstruktors
//*** delta_MATRIX& d4 = *new delta_MATRIX(3, 3, 3, 3);                          // (3x3x3x3) delta-Matrix
//
//*** std::cout<<"Probe: (3x3x3x3) delta-Matrix 4. Stufe d4(i,k,l,m) = " << d4 << endl; // drucken
//
// Freigabe der dynamischen Matrix d4
//*** delete& d4;                                                                // löschen mittels Destruktor '~'
// ------------------------------------------------------------------------

Test von Effektivität und Genauigkeit

Pause -> Zeichen eingeben:�

(1x3) Vektor e0 = [ (1*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

1.000e+00  0.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

(1x3) Vektor e1 = [ (1*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

0.000e+00  1.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

(1x3) Vektor v0 = [ (1*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

6.000e+00  2.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

(1x3) Vektor v1 = [ (1*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

1.000e+00  2.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

Skalar s = v0 * v1 = 10

Dyadic product M = v0 * v1 = [ (3*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

6.000e+00  1.200e+01  0.000e+00
2.000e+00  4.000e+00  0.000e+00
0.000e+00  0.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

(3x3x3) e-Matrix e = [ (3*3) * (3*1) ]-MATRIX:

Spalten 0 bis 2

 0.000e+00  0.000e+00  0.000e+00
 0.000e+00  0.000e+00 -1.000e+00
 0.000e+00  1.000e+00  0.000e+00
 0.000e+00  0.000e+00  1.000e+00
 0.000e+00  0.000e+00  0.000e+00
-1.000e+00  0.000e+00  0.000e+00
 0.000e+00 -1.000e+00  0.000e+00
 1.000e+00  0.000e+00  0.000e+00
 0.000e+00  0.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

(1x3) Vektor e2 = e0 x e1 = [ (1*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

0.000e+00  0.000e+00  1.000e+01

Pause -> Zeichen eingeben:�

(2x2) Matrix A = [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

0.000e+00  2.000e+00
1.000e+00  3.000e+00

Pause -> Zeichen eingeben:�

(2x2) Matrix B = [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

1.000e+00  5.000e-01
0.000e+00  1.000e+00

Pause -> Zeichen eingeben:�

(2x2) Matrix C = A * B = [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

0.000e+00  2.000e+00
1.000e+00  3.500e+00

Pause -> Zeichen eingeben:�

(2x2) Matrix C = A(T) * B = [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

0.000e+00  1.000e+00
2.000e+00  4.000e+00

Pause -> Zeichen eingeben:�

(3x3x3x3) delta-Matrix of 4. Order D4(i,k,l,m) = [ (3*3) * (3*3) ]-MATRIX:

Spalten 0 bis 5

0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00
0.000e+00  1.000e+00  0.000e+00 -1.000e+00  0.000e+00  0.000e+00
0.000e+00  0.000e+00  1.000e+00  0.000e+00  0.000e+00  0.000e+00
0.000e+00 -1.000e+00  0.000e+00  1.000e+00  0.000e+00  0.000e+00
0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00
0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  1.000e+00
0.000e+00  0.000e+00 -1.000e+00  0.000e+00  0.000e+00  0.000e+00
0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00 -1.000e+00
0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00

Spalten 6 bis 8

 0.000e+00  0.000e+00  0.000e+00
 0.000e+00  0.000e+00  0.000e+00
-1.000e+00  0.000e+00  0.000e+00
 0.000e+00  0.000e+00  0.000e+00
 0.000e+00  0.000e+00  0.000e+00
 0.000e+00 -1.000e+00  0.000e+00
 1.000e+00  0.000e+00  0.000e+00
 0.000e+00  1.000e+00  0.000e+00
 0.000e+00  0.000e+00  0.000e+00

Basisvektoren und Metriken im euklidischen Raum

Für das Tensorkalkül werden zwei duale, schiefwinklige nichtnormierte Bezugssysteme mit kovarianten Basisvektoren

${\vec {g}}{_{\alpha }}=\left({\begin{array}{ccc}{\vec {g}}{_{0}}&{\vec {g}}_{1}&{\vec {g}}_{2}\end{array}}\right)={g{^{i}}}_{\alpha }*{\vec {e}}{_{i}}$

einerseits und mit kontravarianten Basisvektoren

${\vec {g}}{^{\beta }}=\left({\begin{array}{ccc}{\vec {g}}{^{0}}&{\vec {g}}{^{1}}&{\vec {g}}{^{2}}\end{array}}\right)=g{_{i}}^{\beta }*{\vec {e}}{^{i}}$

andererseits in das kartesische Bezugssystem :${\vec {e}}{_{i}}={\vec {e}}{^{i}}$ eingebettet, dessen Komponenten mit lateinischen Indizes gekennzeichnet sind.

Die Bezeichnungen der schiefwinkligen Objekte werden hingegen zwecks systematischer Unterscheidung mit griechischen Buchstaben markiert.

Die beiden schiefwinkligen Bezugssysteme sind über die Einheitsmatrix :$\delta {_{\alpha }}^{\beta }$ zueinander orthogonal ausgerichtet

${\vec {g}}_{\alpha }\cdot {\vec {g}}^{\beta }=g{^{i}}{_{\alpha }}*g{_{i}}{^{\beta }}=\delta {_{\alpha }}^{\beta }.$

Die Zahlenwerte werden hier im Programmcode von den physikalischen Vektoren ${\vec {g_{0}}}={\vec {v0}}$ and ${\vec {g_{1}}}={\vec {v1}}$ übernommen.

Im nächsten Schritt wird damit die kovariante Metrik (Maßtensor)

$g_{\alpha \beta }=g{^{i}}{_{\alpha }}*g{^{i}}{_{\beta }}$

ermittelt. Die Volumenform g gewinnt man aus der Determinante

$g=det(g_{\alpha \beta })\;.$

Die kontravariante Metrik erfolgt durch Inversion der kovarianten Metrik

$g^{\alpha \gamma }*g_{\gamma \beta }=\delta {^{\alpha }}_{\beta }\;.$

Daraus werden abschließend die kontravarianten Basisvektoren

$g{_{i}}^{\alpha }=g^{\alpha \gamma }*g{^{i}}_{\gamma }.$

ermittelt.

Der folgende Programmcode erläutert die Tensorklasse BASIS und zeigt typische Anwendungsbeispiele der erzeugten Objekte.

Prototypen neuer Klassen und Methoden

cout << endl
     << "3.1 - Basisvektoren und Metriken im Euklidischen Raum" << endl
     << "=====================================================" << endl << endl;
//
// Kovariante Basisvektoren
VEKTOR g0(3), g1(3);
g0 = v0;                                                                        //.dimension(3);
g1 = v1;                                                                        //.dimension(3);
//
// Erzeugen einer dynamischen Instanz der Klasse BASIS mittels Konstruktor
// -----------------------------------------------------------------------
BASIS& basis = *new BASIS(g0,g1);
//
// Kovariante Basisvektoren
// ------------------------
MATRIX cov_basis(3,2);                                                          // (3x2)-Matrix
kov_basis = basis.kovariante_basis();
//
std::cout << "Kovariante Basisvektoren: (3x2) Matrix kov_basis = "
          << kov_basis << endl;                                                 // drucken
//
// Kovariante Metrik
// -----------------
MATRIX kov_metrik(2,2);                                                         // (2x2) Matrix
double g = basis.kovariante_metrik(kov_metrik);
//
std::cout << "Kovariante Metrik: (2x2) Matrix kov_metrik = "
          << kov_metrik << endl;                                                // drucken
//
// Kontravariante Metrik
// ---------------------
MATRIX kontra_metrik(2,2);                                                      // (2x2) Matrix
basis.kontravariante_metrik(kontra_metrik);
//
std::cout << "Kontravariante Metric: (2x2) Matrix Kontra_metrik "
          << kontra_metrik << endl;                                             // drucken
//
// Probe: -> Einheitsmatrix
// ------------------------
MATRIX& Einh = *new MATRIX(2, 2);
Einh = kov_metrik.dimension(1,2,2) * kontra_metrik.dimension(1,2,2);
//
std::cout << "Probe: Einheitsmatrix = kov_metrik * kontra_metrik "
          << Einh << endl;                                                      // drucken
//
// Kontravariante Basisvektoren
// ----------------------------
MATRIX kontra_basis(3, 2);                                                      // (3x2) Matrix
kontra_basis = basis.kontravariante_basis();
//
std::cout << "Kontravariante Basisvektoren: (3x2) Matrix kontra_basis = "
          << kontra_basis << endl;                                              // drucken
//
// Probe: -> Einheitsmatrix
// ------------------------
dEinh = *new MATRIX(2, 2);           // dynamic (2x2) Matrix
dEinh = kov_basis.dimension(1,3,2) * kontra_basis.dimension(1,3,2);
//
std::cout << "Probe: Einheitsmatrix = kov_basis * kontra_basis "
     << dEinh << endl;                                                          // drucken
//
// Freigaben von dynamischen Objekten
// ---------------------------------
delete& basis;                                                                  // loeschen mittels Destruktor '~'
delete& dEinh;                                                                  // loeschen mittels Destruktor '~'
// ------------------------------------------------------------------------

Test von Effektivität und Genauigkeit

Pause -> Zeichen eingeben:�

Kovariante Basisvektoren: (3x2) Matrix kov_Basis = [ (3*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

6.000e+00  1.000e+00
2.000e+00  2.000e+00
0.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

Kovariante Metrik: (2x2) Matrix kov_metrik = [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

4.000e+01  1.000e+01
1.000e+01  5.000e+00

Pause -> Zeichen eingeben:�

Kovariante Metrik: Determinante g = 100

Kontravariante Metrik: (2x2) Matrix kontra_metrik [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

 5.000e-02 -1.000e-01
-1.000e-01  4.000e-01

Pause -> Zeichen eingeben:�

Probe: Einheitsmatrix = kov_metrix * kontra_metric [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

1.000e+00  0.000e+00
2.220e-16  1.000e+00

Pause -> Zeichen eingeben:�

Kontravariante Basisvektoren: (3x2) Matrix kontra_basis = [ (3*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

 2.000e-01 -2.000e-01
-1.000e-01  6.000e-01
 0.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

Probe: Einheitsmatrix = kov_basis * kontra_basis [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

1.000e+00 -4.441e-16
5.551e-17  1.000e+00

Tensortransformationen

Die Transformation der Tensorkomponenten zwischen dem kartesischen Bezugssystem und den schiefwinkligen Bezugssystemen erfolgt über die Komponenten

$g{^{i}}_{\alpha }$ bzw. $g{_{k}}^{\beta }$

der kovarianten bzw. kontravarianten Basisvektoren

${\vec {g}}{_{\alpha }}={\vec {e}}{_{i}}*g{^{i}}{_{\alpha }}$ resp. ${\vec {g}}{^{\beta }}={\vec {e}}{^{k}}*g{_{k}}^{\beta }.$

Nach dieser Maßgabe wird jeder Index eines Tensors nach Bedarf einzeln transformiert. So transformieren sich Tensoren n-ter Stufe

$v{_{\alpha }}=v_{i}*g{^{i}}{_{\alpha }}$ bzw. $v{^{\beta }}=v^{k}*g{_{k}}^{\beta }$

nach dieser Regel

$E_{\alpha \beta \delta \gamma }=E_{iklm}*{g{^{i}}_{\alpha }}*{g{^{k}}_{\beta }}*{g{^{l}}_{\delta }}*{g{^{m}}_{\gamma }}$ bzw. $E^{\alpha \beta \delta \gamma }=E^{iklm}*{g{_{i}}^{\alpha }}*{g{_{k}}^{\beta }}*{g{_{l}}^{\delta }}*{g{_{m}}^{\gamma }}\;$

n mal.

Die kovarianten bzw. kontravarianten Levi-Civita-Tensoren (ε-Tensoren) im schiefwinkligen Bezugssystem ergeben sich daraus, jedoch im Ergebnis auch vereinfacht aus den e-Matrixen durch Multiplikation mit der Volumenform (s. Tensorkalkül) zu

$\varepsilon _{\alpha \beta }=e_{\alpha \beta }*{\sqrt {g}}$ bzw. $\varepsilon ^{\alpha \beta }=e^{\alpha \beta }*(1/{\sqrt {g}})$

und

$\varepsilon _{\alpha \beta \delta }=e_{\alpha \beta \delta }*{\sqrt {g}}$ bzw. $\varepsilon ^{\alpha \beta \delta }=e^{\alpha \beta \delta }*(1/{\sqrt {g}})\;.$

Der folgende Programmcode erläutert einige typische Beispiele für solche Transformationen von Tensoren durch entsprechende Multiplikationen mit VEKTOR- und MATRIX-Objekten.

Prototypen neuer Klassen und Methoden

cout << endl
     << "4.1 - Tensortransformationen" << endl
     << "============================" << endl << endl;
//
Vektor v(3);                                                                     // (1x3)-Vektor
// Vektor im kartesischen Bezugssystem e0-e1-e2
v(0) = 0.0;    v(1) = 5.0;
//
cout << "kartesisches System: Vektor v(i) = " << v << endl;                         // drucken
//
// Vektor v(i) ins schiefwinklige Bezugssystem g0-g1 transformieren
// ----------------------------------------------------------------
Vektor v_kov(2);                                                                 // (1x2) Vektor
// v_kov(0,alpha) = kov_basis(0,i,alpha) * v(0,i)
//
v_kov = kov_basis.dimension(1,3,2) * v;
//
cout << "schiefwinkliges System: kovarianter Vektor v(alpha) = " << v_kov << endl;  // drucken
//
// Vektor v_schief(alpha) vom schiefwinkligen ins kartesische System zurücktransformieren
// --------------------------------------------------------------------
// v(0,i) = contra_basis(i,alpha) * v_cov(0,alpha)
v = kontra_basis * v_kov;
//
cout << "Probe: Kartesisches System: Vektor v(i) = " << v << endl;               // drucken
//
// Tensor 4.Stufe transformieren
// -----------------------------
// Elastiztaestensor E im kartesischen Bezugssystem e0-e1 aufbauen
// C-Feld mit fiktiven Werten zeilenweise (!) vorbelegen
// ---------------------------------------------------------------
double F16[16] = { 1,0,0,0, 0,0.5,0.5,0, 0, 0.5, 0.5, 0, 0,0,0,1 };
//
// Erzeugen einer statischen Instanz der Klasse MATRIX mittels Konstruktor
// -----------------------------------------------------------------------
MATRIX E(F16,2, 2, 2, 2);       // (2x2x2x2) Matrix wird mit F16 vorbelegt
//
cout << "(2x2x2x2) Elastizitaetstensor E := F16 = " << E <<endl;                 // drucken
//
MATRIX E_contra(2, 2, 2, 2);                                                     // (2x2x2x2)-Matrix
//
// Tensor E(i,k,l,m) ins schiefwinklige Bezugssystem g0-g1 transformieren
// ----------------------------------------------------------------------
// kontravariante Basis auf (2x2)-Matrix kuerzen - dritte Zeile streichen
//
kontra_basis = kontra_basis.untermatrix_heraus(0,2,0,2);
cout << "Kontravariante Basisvektoren: (2x2) Matrix kontra_basis = "
     << kontra_basis.dimension(2,2)     << endl;                                 // drucken
//
// E_kontra(i,k, gamma,delta) =
//         (E(IK,l,m) * kontra_basis(1,l,gamma)) * kontra_basis(0,m,delta)
//                                                     mit N{IK} = N{i}*N{k}
// E_kontra =        E.dimension(2 * 2, 2, 2) * kontra_basis.dimension(1, 2, 2);
//
E_kontra = E_kontra.dimension(2 * 2, 2, 2) * kontra_basis.dimension(1, 2, 2);
//
// E_kontra(alpha,beta, gamma,delta) =
//   kontra_basis(0,i,alpha) * kontra_basis(0,k,beta)
//                           * E_kontra(i,k,GD)
//                                   with N{GD} = N{gamma}*N{delta}
// E_kontra = kontra_basis.dimension(1, 2, 2) * E_kontra.dimension(2, 2, 2 * 2);
//
E_kontra = kontra_basis.dimension(1, 2, 2) * E_kontra.dimension(2, 2, 2 * 2);
//
cout <<"schiefwinkliges System: Tensor "`
     << "E_kontra(alpha,beta,gamma,delta) = " << E_kontra << endl;               // drucken
//
// kontravarianten (2x2)-Levi-Civita-Tensor mittels Klasse e_MATRIX erzeugen
// --------------------------------------------------------------------
e_MATRIX epsilon_Tensor(2,2,sqrt(g));
//
cout<< "schiefwinkliges System: kontravariant epsilon-Tensor = "
                             << epsilon_Tensor << endl;                          // drucken
// -------------------------------------------------------------------

Test von Effektivität und Genauigkeit

Pause -> Zeichen eingeben:�

kartesisches System: Vektor v(i) = [ (1*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

0.000e+00  5.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

schiefwinkliges System: kovarianter Vektor v(alpha) = [ (1*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

1.000e+01  1.000e+01

Pause -> Zeichen eingeben:�

Probe: kartesisches System: Vektor v(i) = [ (1*3) * (1*1) ]-MATRIX:

Spalten 0 bis 2

-4.441e-16 5.000e+00 0.000e+00

Pause -> Zeichen eingeben:�

(2x2x2x2) Elastizitaetstensor E := F16 = [ (2*2) * (2*2) ]-MATRIX:

Spalten 0 bis 3

1.000e+00  0.000e+00  0.000e+00  0.000e+00
0.000e+00  5.000e-01  5.000e-01  0.000e+00
0.000e+00  5.000e-01  5.000e-01  0.000e+00
0.000e+00  0.000e+00  0.000e+00  1.000e+00

Pause -> Zeichen eingeben:�

kontravariante Basisvektoren (2x2) Matrix kontra_Basis = [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

 2.000e-01 -2.000e-01
-1.000e-01  6.000e-01

Pause -> Zeichen eingeben:�

schiefwinkliges System: Tensor E_kontra(alpha,beta,gamma,delta) = [ (2*2) * (2*2) ]-MATRIX:

Spalten 0 bis 3

 2.500e-03 -5.000e-03 -5.000e-03  1.000e-02
-5.000e-03  1.500e-02  1.500e-02 -4.000e-02
-5.000e-03  1.500e-02  1.500e-02 -4.000e-02
 1.000e-02 -4.000e-02 -4.000e-02  1.600e-01

Pause -> Zeichen eingeben:�

schiefwinkliges System: kontravarianter epsilon-Tensor = [ (2*2) * (1*1) ]-MATRIX:

Spalten 0 bis 1

 0.000e+00  1.000e-01
-1.000e-01  0.000e+00

Zusammenfassung

Das bewährte Konzept bietet einen geeigneten Rahmen für die indexbasierte Tensor- und Matrix-Algebra, die nahtlos in eine benutzerfreundliche objektorientierte Formulierung der Matrixarithmetik in C++ Algorithmen überführt wird. Speziell eröffnet dabei die deklarierte allgemeine Matrizenmultiplikation flexible Operationen für mehrdimensionale Matrizen, wobei als Vorteil die Kommutativität der Matrixoperanden nutzbar bleibt.

Es ermöglicht darüber hinaus auch durch die indexbasierten Benutzerschnittstellen der verschiedenen Klassen besonders flexible Matrix- und Tensor-Operationen, die in jeder Hinsicht erweiterbar sind. Ergänzende Software-Quellen können im WWW für den akademischen Gebrauch und die private Nutzung heruntergeladen werden Damit wird auf synoptische Weise ein weiter Bereich von Anwendungen des Matrixkalküls erschlossen, der es auf zuverlässige Weise ermöglicht, konsistente Operationen und Transformationen von Tensoren in euklidischen Räumen durchzuführen.

Erst durch die Integration von Tensoranalysis, Matrizenalgebra und numerischen objektorientierten Methoden für Objekte höherer Stufe wird schließlich diese kohäsive Benutzerschnittstelle geschaffen. Der Artikel bietet damit, über die spezifischen Kapitel verteilt, eine umfassende Darstellung von linearer Algebra und geeigneter Datenverwaltung für indexbasierte mathematische Objekte sowie prototypische Lösungsansätze durch objektorientierte Programmstrukturen. Zum Nachweis von Effektivität und Präzision der numerischen Algorithmen werden typische Beispiele der linearen Algebra gezeigt.

Literatur

Diese artikel "Matrix- und Tensor-Objekte höherer Stufe für numerische Simulationen" ist von Wikipedia The list of its authors can be seen in its historical and/or the page Edithistory:Matrix- und Tensor-Objekte höherer Stufe für numerische Simulationen.

📰 Article(s) of the same category(ies)

• ORDA
• Arbeitskreis philosophierender Ingenieure und Naturwissenschaftler
• Tree Engineering
• Think ING.
• DFMA

---

 Facebook Page

 Follow us on Twitter !

Read or create/edit this page in another language

• Matrix- und Tensor-Objekte höherer Stufe für numerische Simulationen in English