Benutzer:Fritzudo/Object-oriented models of higher-order matrices and tensors

Object-oriented models of higher-order matrices and tensors for numerical simulations

Through numerical simulations, important insights into the design of natural and technical processes are gained in the fields of natural sciences and technology using specific computer models. This often involves algorithms based on the linear algebra of traditional matrix operations, found in corresponding libraries of procedural programming. Following the principles of modern software design, object-oriented modeling and programming have developed an alternative to matrix classes since the 1990s, which also capture multidimensional matrix objects and provide standard operators for arithmetic and functional operations on complex objects. This is achieved through a selected hierarchy of matrix and tensor objects, whose structure and properties are encapsulated according to object-oriented principles and can be inherited internally by subclasses. Communication with these objects is done exclusively through public member functions, which, along with overloaded operators of index-based matrix and tensor arithmetic, form an intuitively understandable user interface.

Introduction

The development of object-oriented analysis and modeling has had a significant impact on the applications of cooperative planning and construction in engineering, extending to current product and building information modeling (BIM). Specifically, numerical models for simulating physical processes in solid-state and fluid mechanics play a crucial role in understanding complex interactions in nature and technology. In general, the physical principles are translated into reliable algorithms of various approximation methods, such as finite differences or finite element methods, using tensor analysis, matrix calculus, and numerical program structures.

Algorithms in object-oriented structures are highly relevant because the underlying paradigms of encapsulation and inheritance ensure transparency of formulation, reliability of methods, and clearly defined communication with objects. Object orientation opens up new possibilities to integrate the theoretical foundations of tensor analysis with the numerical advantages of matrix calculus and algorithmic computational programs into a methodical unity.

The classical matrix computation is based on a concept with a two-dimensional arrangement of computational quantities that dates back to a time when digital computers did not exist. However, it is now standard in programming to use multidimensional arrays to dynamically store and process data. To optimally utilize these resources, it is necessary to extend matrix computation to multidimensional matrices and formulate suitable syntax and semantics that adapt matrix arithmetic to this concept. The proven index notation of tensor calculus lends itself ideally to this task, as shown in the following example

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

or

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

where symbolic matrix notation unfortunately fails. Through the object-oriented class concept, appropriate structures for new tensor and matrix classes are created, enabling the definition of index-based matrix operations for multidimensional matrices to generalize matrix computation to higher-level objects and unify the nomenclature of tensors, matrices, and arrays. Overloading arithmetic and functional standard operators for matrix objects significantly enhances user-friendliness of interfaces.

This established concept is illustrated in subsequent chapters through manageable examples of linear algebra using well-structured prototypes of program structures. This allows processes to be transparently demonstrated and numerically verified. The methods of object-oriented programming are extensively described explicitly. Comprehensive introductions to the C++ programming language can be found in various sources. .

Matrix objects and overloaded operators

The following assumes that rectangular matrices are stored \textbf{\textit{column-wise}}, as is customary in scientific and technical applications and algorithmic programming languages such as FORTRAN.

${\textbf {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}$ with dimensions ($N_{i}\times N_{k})$.

Here, the first index $i$ runs over the rows, and the second index $k$ runs over the columns of matrix $M_{ik}$.

Internally, matrices are stored by assigning matrix elements in sequential order to a one-dimensional vector field:

${\textbf {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}$ with dimension N_{z} = N_{i} * N_{k},

similar to how it is done internally in von Neumann computer architectures.

The corresponding bijective mappings are defined as follows:

The mapping

$V_{z}:=M\underbrace {_{ik}} _{z}$ is given by $z=k*N_{i}+i$,

and the inverse mapping

$M\underbrace {_{ik}} _{z}:=V_{z}$ is given by $k={\text{Entier}}(z/N_{i})$ and $i=z-k*N_{i}$.

The object descriptor for representing the matrix is independent of the internal sequential storage of elements and only controls the external display of the matrix form.

If matrix M is reduced to a single row, it becomes a row vector $M_{0k}$; the same applies for a column vector $M_{i0}$. This demonstrates that it is possible to incorporate \textbf{0-indices} into the notation without changing the storage of vectors or matrices.

This holds true for vector V as well: $V_{z}:=V_{0z}:=V_{z0}$ and for any matrices M:

$M_{ik}:=M_{0ik}:=M_{i0k}:=M_{ik0}.$

The MATRIX M constructor $(N_{i},N_{k},N_{l},N_{m})$ defines the mathematical form of a rectangular matrix in the program code by choosing the dimensions. Currently, up to four dimensions are permitted, which are generally sufficient for engineering applications. Additionally, adjacent indices can be combined using the mappings, as utilized in the program code.

The following C++ program example explains how the matrix classes are handled to create corresponding matrix objects (instances). The constructors of the VEKTOR, MATRIX, delta_MATRIX, e_MATRIX classes can be used to create static or dynamic objects with elements of type 'double'. Using a C-style array notation, vector and matrix elements can be preloaded with numerical values. Input/output is facilitated through overloaded standard operators '>>' and '<<'.

Prototyping of new object-oriented classes and methods: Creation of Matrix Objects and Handling of overloaded Operators

cout << endl

    << "1.1 - Matrix Objects and overloaded operators" << endl
    << "=============================================" << endl <<endl; 

// Preallocate of C-array with values

// ----------------------------------

double F[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

cout << "C-Vector F(10) ="; for (int j = 0; j < 10; j++) cout << " " << F[j]; // printing

cout << endl << endl;

// Creating a static instance of the class VEKTOR using constructor

// ----------------------------------------------------------------

VEKTOR V(F, 6); // (1x6)-vector V is preloaded with F

cout << "(1x6) Vector V := F = " << V << endl; // overloaded operator '<<'

// Creating a static instance of the class MATRIX using constructor

MATRIX M(3, 3); // (3x3)-matrix, assignment of V to M

// Exception: type casting by M := V for demonstration

M = V; // overloaded operator '='

cout << "(3x3) Matrix M := V = " << M << endl; // printing

// Creating a static instance of class delta_MATRIX using constructor

delta_MATRIX d(2, 2); // (2x2)-delta-matrix

cout << "(2x2) delta-Matrix d = " << d << endl; // printing

// Creating a static instance of the MATRIX class using constructor

MATRIX A(2, 2); // (2x2)-matrix

cout << "(2x2) Matrix A = " << A << endl; // printing

Testing effectiveness and accuracy: Matrix Objects and elementary Operators

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

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

Columns 0 to 5

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

Pause -> Enter a character:�

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

Columns 0 to 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 -> Enter a character:� (2x2) delta-Matrix d = [ (2x2) * (1x1) ]-MATRIX:

Columns 0 to 1

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

Pause -> Enter a character:�

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

Columns 0 to 1

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

The initially defined dimensions of a matrix object M can, however, be temporarily modified in the present concept for the meaningful execution of overloaded operators. This is done by calling the member function M.dimension $(N_{i},N_{k},N_{l},N_{m})$ bound to the object, and applies only until the end of the respective operation!

This function is especially important for implementing the general matrix multiplication - as will be shown later.

The following example first demonstrates elementary applications of the overloaded arithmetic member functions '+', '-', '=', as well as the determinant member function with the declared binding to the defined class instances (objects).

Prototyping of new object-oriented classes and methods: Addition and Subtraction

cout << endl

    << "1.2 - Addition und Subtraction" << endl
    << "==============================" << endl << endl;

// Matrix M is temporarily resized by the member function dimension !

// Matrix Addition M(i,k) + A(i,k)

A = M.dimension(2,2) + d; // overloaded operator '+'

cout << "Matrix A = A + d = " << A << endl; // print

cout << "Original (3x3) Matrix preserved: M = " << M << endl; // print

//

// Matrix Subtraction A(i,k) - d(i,k)

// ----------------------------------

A = A - d; // overloaded operator '-'

cout << "Matrix A = A - d = " << A << endl; // print

// Determinant of Matrix A

double S;

S = A.determinant(2);

cout << "Determinant det(A) = " << S << endl << endl;

Testing effectiveness and accuracy: Addition and Subtraction

Pause -> Enter a character:�

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

Columns 0 to 1

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

Pause -> Zeichen eingeben:�

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

Columns 0 to 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 -> Enter a character:�

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

Columns 0 to 1

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

Pause -> Enter a character:�

Determinant det(A) = -2

Cartesian basis and matrix vector multiplications

For the multiplication of two rectangular matrices A and B, there are generally several variants, such as A * B, A{^T} * B, A * B^T, A^T * B^T, which usually require different procedures.

These variants are summarized here in the general matrix multiplication with the four-dimensional result matrix C

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

where the summation index \textbf{s} is to be summed over.

The corresponding instruction in the C++ program has the simple form, with the matrix dimensions properly agreed upon in the constructor and the overloaded arithmetic operator '*':

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

For general cases where the dimensions need to be adjusted, the modification of this instruction typically involves using the \textbf{member function dimension} as follows:

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

where its dimensions are preloaded with $1$ and the application is optional.

When using this algorithm on vectors and matrices with different dimensions, it is important to observe the following rules:

as already explained, adjacent indices can be temporarily combined using the dimension function - only for the execution of these operations - without changing the original dimensions or the internal storage of the matrix elements. This allows any matrices to be adapted to this schema. The same applies to the summation index,

for individual indices, the 0-index can also be used. This means that at this point, Dimension 1 is set and assigned to the matrix. This measure also does not change the internal storage scheme,

both measures ensure that the summation index s is always placed in the middle position,

- due to the temporary change of dimensions with the dimension function and also due to the sequence of multiple multiplications determined by the compiler, this general multiplication should not be chained and typically needs to be executed with a single statement,

- Vectors of the VEKTOR class v(N_i) are derived from the MATRIX class and are always internally defined as row vectors v(1,N_i). This measure simplifies scalar products, as it presets the summation index to the middle position - optimizing handling. The use of the dimension function is also permissible here.

The basis vectors of a Cartesian reference system (usually referred to as an x-y-z system) are denoted by Latin indices:

${\vec {e}}{^{i}}={\vec {e}}{_{i}}=({\vec {e}}{_{0}}\;{\vec {e}}{_{1}}\;{\vec {e}}{_{2}})$ with components :${e_{0}}^{i}=[1,0,0]),\;{e_{1}}^{i}=[0,1,0]$, and ${e_{2}}_{i}=[(0,0,1]\;$

are arranged in the following manner in this spatial reference system:

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

The \textbf{scalar product} of the vectors is then

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

while the \textbf{dyadic product} is obtained from the multiplication:

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

The Levi-Civita symbol - also known as the e-matrix in the Cartesian system - is defined as follows by the antisymmetric $(2\times 2)$-matrix or $(3\times 3\times 3)$-matrix with numerical values ±1:

$e_{ik}=\left({\begin{array}{lc}0&1\\-1&0\\\end{array}}\right)$ or $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).$

Using the e-matrix, the cross product of vectors ${e_{0}}^{k}$ and ${e_{1}}^{l}$ in the Cartesian system can be determined as follows:

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

\newline

As examples of typical matrix multiplications, the following operations with $(2\times 2)$ rectangular matrices are selected:

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

Another example concerns the Grassmann identity, cf. Levi-Civita_Symbol#Tensor_Calculus.

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

which leads to the delta matrix of the 4th order.

The following program code illustrates these examples using multidimensional multiplications with corresponding VEKTOR and MATRIX objects.

Prototyping of new object-oriented classes and methods: Cartesian basis and matrix vector multiplications

cout << endl

    << "2.1 - " Cartesian basis and matrix vector multiplications"<< endl
    << "=========================================================" << endl << end; 

// => Vectors v(i) are internally treated as matrices v(0,i) with N(0) = 1,

// in order to build scalar products rather easy !

// -----------------------------------------------------------------

// Creating three Cartesian reference vectors

VEKTOR e0(3), e1(3), e2(3); // (1x3)-Vektoren

e0(0) = 1.0; e1(1) = 1.0;

std::cout << "(1x3) Vektor e0 = " << e0 << endl; // print

std::cout << "(1x3) Vektor e1 = " << e1 << endl; // print

// Creating two physical vectors

VEKTOR v0(3), v1(3); // (1x3)-Vectors

v0(0) = 6.0; v0(1) = 2.0;

std::cout << "(1x3) Vektor v0 = " << v0 << endl; // print

v1(0) = 1.0; v1(1) = 2.0;

std::cout << "(1x3) Vektor v1 = " << v1 << endl; // print

// Calculating the scalar product s = v0(0,s) * v1(0,s)

// ----------------------------------------------------

double s;

s = v0 * v1;

std::cout << "Scalar s = v0 * v1 = " << s << endl << // print

// Calculating the dyadic product M(i,k) = v0(i) * v1(k)

// ------------------------------------------------------

M = v0.dimension(3) * v1.dimension(3);

std::cout << "Dyadic product M = v0 * v1 = " << M << endl; // print

// Calculating the vector product // e2(0,i) = e(ikl) * e0(0,k) * e1(0,l)

// ----------------------------------------------------------------

// Creating a static instance of the class e_MATRIX // using the constructor

e_MATRIX e(3,3,3); // (3x3x3)-e-Matrix

std::cout << "(3x3x3) e-Matrix e = " << e << endl; // print

// Chained multiplication without member function dimension

e2 = (e * e0) * e1;

std::cout << "(1x3) vector e2 = e0 x e1 = " << e2 << endl; // print

// General matrix product C(i,k,l,m) = A(i,s,k) * B(l,s,m)

// -------------------------------------------------------

// Using the member function M.dimension(N{i}, N{k}, N{l}, N{m}),

// => Matrix M(s,k) can be temporarily extended by a 0-Index

// M(s,k) =: M(0,s,k) mit N{0} = 1,

// => Matrix indices M(i,k,l) can be bundled up to M(IK,s,l)

// M(i,k,s,l) =: M(IK,s,l) mit N{IK} = N{i} * N{k},

// => in order to bring the index s in the centered position M(i,s,l) !

// --------------------------------------------------------------------

std::cout << "(2x2) Matrix A = " << A << endl; // print

// Creating the matrices C and B

MATRIX C(2,2), B(2,2); // (2x2) Matrices

              B = d;   B(0,1) = 0.5;

std::cout << "(2x2) Matrix B = " << B << endl; // print

// Matrix multiplication 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; // print

// Matrix multiplication 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; // print

// Grassmann Identity

// ------------------

// delta-Matrix(i,k,l,m) = e-Matrix(i,k,s) * e-Matrix(l,m,s)

// Creating a static instance of the class MATRIX using the constructor

MATRIX D4(3,3,3,3); // (3x3x3x3) Matrix

D4 = e.dimension(3 * 3, 3) * e.dimension(3 * 3, 3);

std::cout << "(3x3x3x3) delta-Matrix of 4. order D4(i,k,l,m) = " << D4 << endl; // print

// Test

// ----

// Creating a dynamic instance of the class delta_MATRIX using the constructor

//*** delta_MATRIX& d4 = *new delta_MATRIX(3, 3, 3, 3); // (3x3x3x3) delta-Matrix

//*** std::cout<<"Test: (3x3x3x3) delta-Matrix of 4. order d4(i,k,l,m) = " << d4 << endl; // print

// Release the dynamic matrix d4

//*** delete& d4; // delete using destructor '~'

// -------------------------------------------------------------------

Testing effectiveness and accuracy: Cartesian basis and matrix vector multiplications

Pause -> Enter a character:�

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

Columns 0 to 2

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

Pause -> Enter a character:�

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

Columns 0 to 2

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

Pause -> Zeichen eingeben:�

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

Columns 0 to 2

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

Pause -> Zeichen eingeben:�

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

Columns 0 to 2

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

Pause -> Enter a character:�

Skalar S = e0 * e1 = 10

Dyadic product M = e0 * e1 = [ (3*3) * (1*1) ]-MATRIX:

Columns 0 to 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 -> Enter a character:�

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

Columns 0 to 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 -> Enter a character:�

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

Columns 0 to 2

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

Pause -> Enter a character:�

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

Columns 0 to 1

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

Pause -> Enter a character:�

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

Columns 0 to 1

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

Pause -> Enter a character:�

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

Columns 0 to 1

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

Pause -> Enter a character:�

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

Columns 0 to 1

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

Pause -> Enter a character:�

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

Columns 0 to 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 

Columns 6 to 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

Basis vectors and metrics in Euclidian space

For tensor calculus, two dual, non-normalized orthogonal reference systems are used, each with covariant basis vectors

${\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}}$

on one hand, and contravariant basis vectors

${\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}}$

on the other hand, embedded in the Cartesian reference system :${\vec {e}}{_{i}}={\vec {e}}{^{i}}$.

All notations of the orthogonal objects are now marked with Greek letters for distinction from the Cartesian system. \newline

The two orthogonal reference systems are oriented orthogonally to each other via the unit matrix :$\delta {_{\alpha }}^{\beta }$

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

The numerical values are taken here in the program code from the physical vectors ${\vec {g_{0}}}={\vec {v0}}$ and ${\vec {g_{1}}}={\vec {v1}}$.

Next, the covariant metric (measure tensor)

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

is determined. The volume form g is obtained from the determinant

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

The contravariant metric is obtained by inverting the covariant metric

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

From this, the contravariant basis vectors are finally determined as

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

The following program code explains typical examples of BASIS objects.

Prototyping of new object-oriented classes and methods: Calculation of Basis Vectors and Metrics in Euclidian Space

cout << endl

    << "3.1 - Basis Vectors and Metrics in Euklidian Space" << endl
    << "==================================================" << endl << endl;

// Covariant basis vectors

VEKTOR g0(3), g1(3);

g0 = v0; //.dimension(3);

g1 = v1; //.dimension(3);

// Create a dynamic instance of the class BASIS using constructor

// --------------------------------------------------------------

BASIS& basis = *new BASIS(g0,g1);

// Covariant basis vectors // ------------------------

MATRIX cov_basis(3,2); // (3x2)-Matrix

// cov_basis = covariante_basis(); // upcoming

   cov_basis = basis.kovariante_basis();

std::cout << "Covariant basis vectors: (3x2) Matrix cov_basis = "

    << cov_basis << endl;                                     // print

// Covariant metric

// -----------------

MATRIX cov_metrik(2,2); // (2x2) Matrix

// double g = basis.covariante_metrik(cov_metrik); // upcoming

   double g = basis.kovariante_metrik(cov_metrik);

std::cout << "Covariant metric: (2x2) Matrix cov_metrik = "

    << cov_metrik << endl;  // print

std::cout << "Covariant metric: Determinant g = " << g << endl << endl;

                                                             // print

// Contravariant metric

// ---------------------

MATRIX contra_metrik(2,2); // (2x2) Matrix

// basis.contravariante_metrik(contra_metrik); // upcoming

   basis.kontravariante_metrik(contra_metrik);

std::cout << "Contravariant metric: (2x2) Matrix contra_metrik "

    << contra_metrik << endl;                                // print

// Test -> Identity matrix

// -----------------------

MATRIX& Einh = *new MATRIX(2, 2);

Einh = cov_metrik.dimension(1,2,2) * contra_metrik.dimension(1,2,2);

std::cout << "Test: Identity matrix = cov_metrik * contra_metrik "

    << Einh << endl;                                         // print

// Contravariant basis vectors

// ----------------------------

MATRIX contra_basis(3, 2); // (3x2) Matrix

// contra_basis = basis.contravariante_basis(); // upcoming

   contra_basis = basis.kontravariante_basis();

std::cout << "Contravariant basis vectors: (3x2) Matrix contra_basis = "

    << contra_basis << endl;                                 // print

// Test -> Identity matrix

// ------------------------

Einh = *new MATRIX(2, 2); // dynamic (2x2) Matrix

Einh = cov_basis.dimension(1,3,2) * contra_basis.dimension(1,3,2);

std::cout << "Test: Identity matrix = cov_basis * contra_basis "

    << Einh << endl;                                          // print

// Release dynamic objects

// -----------------------

delete& basis; // delete using destructor '~'

delete& Einh; // delete using destructor '~'

// -------------------------------------------------------------------

Testing effectiveness and accuracy: Basis Vectors and Metrics in Euclidian Space

Pause -> Enter a character:�

Covariant basis vectors: (3x2) Matrix cov_Basis = [ (3*2) * (1*1) ]-MATRIX:

Columns 0 to 1

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

Pause -> Enter a character:�

Covariant metric: (2x2) Matrix cov_metrik = [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

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

Pause -> Enter a character:�

Covariant metric: Determinant g = 100

Contravariant metric: (2x2) Matrix contra_metrik [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

5.000e-02 -1.000e-01 

-1.000e-01 4.000e-01

Pause -> Enter a character:�

Test: Identity matrix = cov_metrix * contra_metric [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

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

Pause -> Enter a character:�

Contravariant basis vectors: (3x2) Matrix contra_basis = [ (3*2) * (1*1) ]-MATRIX:

Columns 0 to 1

2.000e-01 -2.000e-01 

-1.000e-01 6.000e-01

0.000e+00  0.000e+00 

Pause -> Enter a character:�

Test: Identity matrix = cov_basis * contra_basis [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

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

Tensor transformations

The transformation of tensor components between Cartesian reference systems and skew-angled reference systems is performed using the components

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

of the covariant and contravariant basis vectors

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

Following this rule, each index of a tensor is transformed individually as needed. Thus, tensors of nth order transform according to this rule as

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

and

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

n times.

The covariant and contravariant Levi-Civita tensors ($\varepsilon$-tensor) in the skew-angled reference system are derived accordingly, but simplified from the e-matrices by multiplication with the volume form

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

The following program code explains typical examples of tensor transformations.

Prototyping of new object-oriented classes and methods: Tensor Transformation

cout << endl

    << "4.1 - Tensor Transformations" << endl
    << "============================" << endl << endl;

VEKTOR v(3); // (1x3)-Vektor

// Vector in cartesian reference system e0-e1-e2 v(0) = 0.0; v(1) = 5.0;

cout << "Cartesian system: Vector v(i) = " << v << endl; // print

// Transform vector v(i) into skew-angled reference system g0-g1

// -------------------------------------------------------------

VEKTOR v_cov(2); // (1x2) vector

// v_cov(0,alpha) = cov_basis(0,i,alpha) * v(0,i) //

v_cov = cov_basis.dimension(1,3,2) * v;

cout << "Skew-angled system: Covariant vector v(alpha) = " << v_cov << endl; // print

// Transform vector v_schief(alpha) from skew-angled to Cartesian system

// --------------------------------------------------------------------

// v(0,i) = contra_basis(i,alpha) * v_cov(0,alpha)

v = contra_basis * v_cov;

cout << "Test: Cartesian system: Vector v(i) = " << v << endl; // print

// Transform 4th-order tensor

// --------------------------

// Elasticity tensor E in Cartesian coordinate system e0-e1

// Allocated C-field with fictious values row by row (!)

// --------------------------------------------------------

double F16[16] = { 1,0,0,0, 0,0.5,0.5,0, 0, 0.5, 0.5, 0, 0,0,0,1 };

// Create a static instance of the MATRIX class using constructor

// ----------------------------------------------------------------

MATRIX E(F16,2, 2, 2, 2); // (2x2x2x2) matrix pre-filled with F16

cout << "(2x2x2x2) Elasticity tensor E := F16 = " << E <<endl; // print

MATRIX E_contra(2, 2, 2, 2); // (2x2x2x2)-Matrix

// Transform tensor E(i,k,l,m) into skew-angled reference system g0-g1

// -------------------------------------------------------------------

// Reduce contravariant basis to (2x2) matrix - eliminate third row

// contra_basis = contra_basis.submatrix_out(0,2,0,2); // upcoming

  contra_basis = contra_basis.untermatrix_heraus(0,2,0,2);    

cout << "Contravariant basis vectors: (2x2) Matrix contra_basis = " << contra_basis.dimension(2,2) << endl; // print

// E_contra(i,k, gamma,delta) = (E(IK,l,m) * contra_basis(1,l,gamma)) * contra_basis(0,m,delta) // with N{IK} = N{i}*N{k}

E_contra = E.dimension(2 * 2, 2, 2) * contra_basis.dimension(1, 2, 2);

E_contra = E_contra.dimension(2 * 2, 2, 2) * contra_basis.dimension(1, 2, 2);

// E_contra(alpha,beta, gamma,delta) = contra_basis(0,i,alpha) * contra_basis(0,k,beta) * E_contra(i,k,GD) // with N{GD} = N{gamma}*N{delta}

E_contra = contra_basis.dimension(1, 2, 2) * E_contra.dimension(2, 2, 2 * 2);

E_contra = contra_basis.dimension(1, 2, 2) * E_contra.dimension(2, 2, 2 * 2);

cout <<"Skew-angled system: Tensor E_contra(alpha,beta,gamma,delta) = " << E_contra << endl; // print

// Generate contravariant (2x2) Levi-Civita tensor using class e_MATRIX

// --------------------------------------------------------------------

e_MATRIX epsilon_Tensor(2,2,sqrt(g));

cout<< "Skew-angled system: Contravariant epsilon-Ttensor = " << epsilon_Tensor << endl; // print

// -------------------------------------------------------------------

Testing effectiveness and accuracy: Tensor Transformations

Pause -> Enter a character:�

Cartesian system: Vector v(i) = [ (1*3) * (1*1) ]-MATRIX:

Columnn 0 to 2

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

Pause -> Enter a character:�

Skew-angled system: Covariante vector v(alpha) = [ (1*2) * (1*1) ]-MATRIX:

Columnn 0 to 1

1.000e+01  1.000e+01 

Pause -> Enter a character:�

Test: Cartesian system: Vector v(i) = [ (1*3) * (1*1) ]-MATRIX:

Columnn 0 to 2

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

Pause -> Enter a character:�

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

Columnn 0 to 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 -> Enter a character:�

Contravariant basis vectors: (2x2) Matrix contra_Basis = [ (2*2) * (1*1) ]-MATRIX:

Columnn 0 to 1

2.000e-01 -2.000e-01 

-1.000e-01 6.000e-01

Pause -> Enter a character:�

Skew-angled system: Tensor E_contra(alpha,beta,gamma,delta) = [ (2*2) * (2*2) ]-MATRIX:

Columnn 0 to 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 -> Enter a character:�

Skew-angled system: Contravariant epsilon-Tensor = [ (2*2) * (1*1) ]-MATRIX:

Columnn 0 to 1

0.000e+00  1.000e-01 

-1.000e-01 0.000e+00

Conclusion

The approach presented herein offers a framework for linear tensor and matrix algebra using index notation, which is then transfered into an object-oriented and user-friendly formulation of matrix arithmetic accessible in C++ program routines. Specifically, the generalized matrix multiplication enables versatile operations on multidimensional matrices while maintaining the commutativity of matrix operands.

Furthermore, it also allows for particularly flexible matrix and tensor operations through the index-based user interfaces of the various classes, which are expandable in every respect. Additional software sources may be downloaded from the WWW for academic and private use. This opens up a wide range of applications of matrix calculus in a synoptic manner, enabling consistent operations and transformations of tensors in Euclidean spaces reliably.

Through this integration of tensor calculus, matrix algebra, and numeric utilizing object-oriented methods for mathematical objects of higher order, a cohesive framework is established. The article provides a comprehensive exposition of linear algebra theory and data management for index-based mathematical objects across relevant chapters, showcases prototypical object-oriented solutions through clearly presented program examples, and offers numerical results as evidence of effectiveness and precision of algorithms.