IMPACT OF MATRIX INVERSION ON THE COMPLEXITY OF THE FINITE ELEMENT METHOD

Institute of Civil and Geoingeenering, Poznan University of Life Sciences, ul. Piątkowska 94, 60-649 Poznan, Poland, email msybis@up.poznan.pl, ORCID 0000-0002-0032-6313 Institute of Civil and Geoingeenering, Poznan University of Life Sciences, ul. Piątkowska 94, 60-649 Poznan, Poland, email asw72@up.poznan.pl, ORCID 0000-0002-7202-6678 Institute of Civil and Geoingeenering, Poznan University of Life Sciences, ul. Piątkowska 94, 60-649 Poznan, Poland, email agraczyk@up.poznan.pl, ORCID 0000-0002-1187-9087


Introduction
The development of a wide construction market and a desire to design innovative architectural building constructions has resulted in the need to create complex numerical models of objects having increasingly higher computational complexity.
There are a lot of numerical methods designed to provide an approximate solutions to the considered problem including finite element method (FEM), boundary element method (BEM) or finite difference method (FDM).
The article is organized as follows.The first section describes the main characteristics of Finite Element Methods for the beam element.The second point concerns the computational complexity of the FEM.In the next point, an estimate of the complexity of the individual stages of the FEM is presented.In the fourth section the basic methods used to solve set of linear equations are presented.The next point is dedicated to the presentation and discussion on the complexity of the implemented algorithms, depending on the implemented matrix inversion method and the number of nodes used in calculations.Summary of the article is included in section six.

Purpose
The purpose of this work is to show that choosing a proper method for solving the set of equations can improve the calculation time (reduce the complexity) by a few levels of magnitude.

Methodology
FEM for the beam element: Finite element method is one of the most popular tools used in engineering calculations.One of the basic assumptions of this method is the abandonment of the analytical division of the area for the division into a finite number of elements [12].The calculation takes place only at specific points, called nodes, and the results for the rest of the area are approximated using, so-called, shape functions on the basis of the results obtained for each node [10], [2], [4].
The basic steps in the FEM for the beam element can be described as follows [7], [8]: − In the first step, we divide the beam in the given number of finite elements and select nodes.We accept both the local and the global coordinate system, − Number of components of the displacement vector q e need to be set, − Shape functions N need to be assumed: [ ] , , , Shape functions for the beam construction can be assumed as a third degree polynomial written in the classical form as: ( ) ( ) ( ) ( ) or in the matrix form: The stiffness matrices K e are determined in each node Calculation of substitute nodal loads matrix Calculation of the B matrix that describes the deflection at each point of the element with the following equation: where L is the differentiation operator described as follows: Substituting equations ( 6) and ( 10) to (9) gives: The stiffness and nodal load matrices are aggregated to the global system with application of the boundary conditions.
We solve the equation for the balance of the structure (12) identifying unknown deflections u: Computational complexity: A key issue regarding the performance of all algorithms, including FEM [7], is the computational complexity or, in other words, computing requirements.We call computational requirements all the necessary arithmetic operations that need to be performed during calculations.
The concept of computational complexity is also associated with the issue of the resources availability.Resources can be defined in the form of time or memory [8], [9].
− Time complexity -the measurement of the time complexity is performed by the calculation of the number of basic operations that depend on the size of the input data.Measuring the actual clock time can vary depending on the implementation of the algorithm, computing machine or a compiler used.
− Memory complexity -the memory comp-lexity expresses the amount of memory used, expressed as the number of memory cells in a function of the size of input data or expressed in bits or bytes, to perform specific computing.
The creators of complexity theory are Juris Hartmanis and Richard Stearns [3].It can be said that the larger the scale of the problem and thus the larger size of the input data, the more resources it will need to perform computing operations [12].
A common way to compare the complexity of various algorithms is to determine their asymptotic growth rate.It is a measure of the time requirements with increasing size of the input data.It can be said that it describes how quickly time requirements are changing.
To describe the asymptotic growth rate socalled O notation is used.It can be said that a function f has a complexity of O(g(n)) if there is a positive constant value c and non-negative value N satisfying the condition: In the remainder of this article the above notation will be used to determine the computational complexity of different stages of FEM algorithm.
Estimating the complexity of the FEM algorithm: In this section we will investigate the complexity of the different stages of the FEM algorithm and we will assess its impact on the complexity of the entire algorithm.a) Dividing the Ω area into finite elements.
In the case of beam element, this step has a very low complexity, which can be estimated as O(n).
b) Adoption of shape functions vector e N and degrees of freedom e q .The complexity of this stage is constant and is equal to O(1).It is due to the fact that regardless the size of the problem (e.g.number of nodes), this step is performed only once.a) Calculation of elementary stiffness matrices e K and nodal loads vectors e Q .This step, similar like the previous one, has the complexity of the order of O(1).This is because the number of selected nodes does not affect the complexity of calculating the elementary matrices e K and e Q .b) Aggregation of the stiffness matrix K and the loads vector Q for the entire object.
The complexity of this stage is linear and thus can be noted as O(n).This is related to the need to add the e K matrix to the aggregated matrix K n-times.b) Boundary conditions.Because taking into account boundary conditions takes place only at the ends of the beam it not depend on the number of nodes.It is therefore assumed that this step has a complexity of the order of O(1).a) Solving Kq Q = equation On the basis of the literature it can be stated that the complexity of this stage is in the order of O(n 3 ).

Methods of solving equations set:
As shown in section 4 the most complex part of each of the algorithms is a step designated to solve the set of equations.This problem is even more complex in the case of finite element methods because this method has to generally deal with matrices of considerable sizes, which entails the need of conducting significant number of mathematical operations.Below, collection of the most popular methods for solving linear sets of equations (and matrix inversion), which can be found in the literature is presented.
The method of algebraic complements: Application of this method requires the appointment of an inverse matrix which subsequently need to be multiplied by the vector of the right side values.The algorithm for inverse matrix calculation can be presented in a following steps [3]: − Calculation of the determinant of the matrix A, − Calculation of the algebra complements of all elements of a matrix, − Transposition of the matrix containing algebraic complements, − Determination of the inverse matrix.These steps can be represented by the following formula: where detA is the determinant of a matrix A and the Gauss-Jordan method: This method is another method whose immediate goal is to determine the inverse matrix (further step in the form of multiplication of inverse matrix and right values vector need to be performed).Gaussian elimination is also called matrix method of the attached identity matrix.
To find the inverse of a matrix A the following set of equations need to be solved: Taking into account that ⎦ is obtained.In order to get the inverse of matrix A the sub-matrix A in the [ ] | A I need to be converted to unitary sub-matrix with the use of elementary operations on the rows.
Gauss elimination: This method was developed by a famous German mathematician Carl Friedrich Gauss [6].Its purpose is to bring the matrix to form a stepwise matrix (Lower or Upper triangular matrix) using elementary operations on the rows.It should be also checked the existence of solution by using the Kronecker-Capelli theorem.
Cholesky decomposition: Every positively definite matrix A can be written as [11,5,13] where: L -the lower-triangular matrix with positive values on the main diagonal, referred to as the «square root» of positive definite matrix A.
In equation ( 17), the unknowns are the elements of the matrix L: Using the formulas for calculating the matrix product it is possible to describe the relations that allow to calculate individual elements of the matrix L .As a result: Expression ( 22) is so-called the Schur complement.Due to the positive value of the expression for positive definite matrix, decomposition algorithm can be applied to the matrix diminished by first row and first column.
LU decomposition: The method of LU decomposition is to produce a matrix A as a product of two triangular matrices: lower L and upper U with the addition of zero elements, above and below, respectively, the main diagonal of the matrix [1].This formally can be written as: where: Analyzed set of equations can be described as: Solving the equation ( 34) is achieved by solving two sets of equations: In this way, we obtain a searched vector x.

Findings
This section provides an analysis of the computational requirements of different matrix inversion algorithms for finite element method.
The calculations were performed for the fixedend beam of the length equal to l x = 6.5 m, loaded with a uniformly distributed load q equal to 15 kN / m, which diagram is shown in Fig. 1.For each method, the calculations of deflection of the beam, which is presented in Fig. 2, divided into 5, 50, 100 and 200 nodes have been performed.The process of calculation of the equation ( 12) has been carried out using the Gauss method, Cholesky decomposition and LU decomposition.In addition, each method has been further optimized.This optimization consisted in the use of some knowledge of matrices.For example, it has been taken into account that the matrices are symmetrical and contain a large number of coefficients equal to 0, which meant that some multiplications and additions can be omitted when performing calculations.The obtained results have been shown in the latter part of this section.Table 1 presents the number of mathematical operations performed during calculation of the deflection of the beam divided into 5 nodes, depending on the implemented method for matrix inversion (solving equation set).
Second column shows the assignment operations necessary to be performed.Assignment operations are the operations of saving the calculated values in the memory of the computer (after the value is calculated it needs to be stored in the memory).The third column shows the number of addition operations and in subsequent columns, number of multiplications, divisions and square root calculations are presented.Last column presents the total number of required operations (sum of the operation of the preceding columns).Hereafter the name basic methods has been used regarding the Gauss, LU decomposition and Cholesky decomposition methods while the term optimized methods has been used regarding the optimized versions of these algorithms.
Comparing the Gauss method to the optimized Gauss method it can be noticed that the number of arithmetical operations performed is reduced by about 31%.For the LU decomposition method and its optimized version complexity reduction is at the level of 43%.In the case of optimized Cholesky decomposition method compared to its non optimized version the profits is equal to approximately 13%.
While taking into account the non optimized versions of the algorithms the number of operations performed by a Cholesky method is the smallest.However, while taking the optimized versions into account more preferably is to use LU and Gauss method.
For 50 nodes, we can observe a significant reduction in the number of mathematical operations required to perform for optimized algorithms compared to basic methods.The obtained results are shown in Table 2. Reduction of the number of arithmetic operations in the Gauss method and LU decomposition compared to their non optimized methods is about 99.5%.The difference between these two optimized methods is however low and is equal to about 2.5% in favor of the LU method.Basic Cholesky method gives results more than 17% better (in terms of required complexity) compared to other basic methods.Optimization of this method gives a gain of about 96%.However, the situation is changing if optimized methods are compared because other methods are more than 80% less complex.For 100 nodes (see Table 3) application of the optimized methods result in a reduction in the number of arithmetic operations performed in relation to the basic methods, even by three orders of magnitude, giving a gain of 99.9%.It is also worth noting that optimized LU decomposition method requires 10 % less operations than the optimized Gauss method and 92% fewer than the Cholesky method.By analyzing the amount of necessary operations to be performed for 200 nodes (Table 4) it can be noticed that the difference between the optimized versions and the basic algorithms is in the order of four orders of magnitude (profit is more than 99.99%) for the Gauss method and LU decomposition.
Profit for the Cholesky decomposition between the basic version and its optimized version reaches the value of 99%.
Similar as for 100 nodes the difference between the optimal versions of LU and Gauss is still around 10% (in favor of the LU decomposition method).Optimized Cholesky method offers the highest complexity compared with other optimized algorithms where number of required operations is more than 95% lower.
The basic version of Cholesky method is about 16% less complex compared to other methods considered.

Originality and practical value
The main achievement of this work is that it shows the impact of the used methodology on the complexity of solving the problem (or equivalently, time needed to obtain results).The difference between the best (the less complex) and the worst (the most complex) is in the row of few orders of magnitude.This result shows that choosing wrong methodology may enlarge time needed to perform calculation significantly.

Conclusions
In the article the impact of the method used for calculating the equations set on the number of arithmetic operations performed is investigated.Deflection of the beam when divided into 5, 50, 100 and 200 nodes using the finite element method has been calculated.To solve the set of equations the Gauss method, LU and Cholesky decomposition methods were implemented.In addition, for each method an optimized version of the algorithm was also investigated, what was intended to further reduce the number of computational operations performed.It was found that in order to maximally reduce the number of necessary to conduct operations at first the advantage of all known properties of the resulting matrix should be taken.Afterwards, the selection of a method that aims at the highest possible reduction of complexity in solving the system of equations should be made.
In the case of optimized methods number of operations can be reduced, by several orders of magnitude.It is also important to choose a proper method because in the basic versions it is preferably to choose Cholesky method while the worst choice is the LU decomposition.However, the situation changes for the optimized methods.The best results are observed for LU decomposition method giving a gain of 3-10% compared to the Gaussian method, and more than 90% compared to the for the Cholesky method.
It was also found that the highest differences between the complexity of the methods are observed if the beam is divided into a larger number of nodes.

11 D , 12 D
, …, nn D represent the algebraic complement values of consecutive elements of matrix A [12].
15) where: A -output matrix, B -inverse matrix, Iidentity matrix.It is necessary for both submatrices [ ] | A I to multiply them by matrix B, to get the matrix [ ]