Interval finite element

Applications of the interval parameters to the modeling of uncertainty

Solution of the following equation

a x = b

where a and b are real numbers is equal to x = b a .

Very often exact values of the parameters a and b are unknown.
Let's assume that a ∈ [ 1 , 2 ] = a and b ∈ [ 1 , 4 ] = b . In that case it is necessary to solve the following equation

[ 1 , 2 ] x = [ 1 , 4 ]

There are several definition of the solution set of the equation with the interval parameters.

United solution set

In this approach the solution is the following set

x = { x : a x = b , a ∈ a , b ∈ b } = b a = [ 1 , 4 ] [ 1 , 2 ] = [ 0.5 , 4 ]

This is the most popular solution set of the interval equation and this solution set will be applied in this article.

In the multidimensional case the united solutions set is much more complicated. Solution set of the following system of linear interval equations

[ [ − 4 , − 3 ] [ − 2 , 2 ] [ − 2 , 2 ] [ − 4 , − 3 ] ] [ x 1 x 2 ] = [ [ − 8 , 8 ] [ − 8 , 8 ] ]

is shown on the following picture

∑ ∃ ∃ ( A , b ) = { x : A x = b , A ∈ A , b ∈ b }

Exact solution set is very complicated, because of that in applications it is necessary to find the smallest interval which contain the exact solution set

♢ ( ∑ ∃ ∃ ( A , b ) ) = ♢ { x : A x = b , A ∈ A , b ∈ b }

or simply

♢ ( ∑ ∃ ∃ ( A , b ) ) = [ x _ 1 , x ¯ 1 ] × [ x _ 2 , x ¯ 2 ] × . . . × [ x _ n , x ¯ n ]

where

x _ i = m i n { x i : A x = b , A ∈ A , b ∈ b } ,     x ¯ i = m a x { x i : A x = b , A ∈ A , b ∈ b }

x i ∈ { x i : A x = b , A ∈ A , b ∈ b } = [ x _ i , x ¯ i ]

See also [2]

Parametric solution set of interval linear system

Interval Finite Element Method require the solution of parameter dependent system of equations (usually with symmetric positive definite matrix). Example of the solution set of general parameter dependent system of equations

[ p 1 p 2 p 2 + 1 p 1 ] [ u 1 u 2 ] = [ p 1 + 6 p 2 5.0 2 p 1 − 6 ] ,       f o r     p 1 ∈ [ 2 , 4 ] , p 2 ∈ [ − 2 , 1 ] .

is shown on the picture below (E. Popova, Parametric Solution Set of Interval Linear System [3]).

Algebraic solution

In this approach x is such interval number for which the equation

[ 1 , 2 ] x = [ 1 , 4 ]

is satisfied. In other words, left side of the equation is equal to the right side of the equation. In this particular case the solution is equal to x = [ 1 , 2 ] because

a x = [ 1 , 2 ] [ 1 , 2 ] = [ 1 , 4 ]

If the uncertainty is bigger i.e. a = [ 1 , 4 ] , then x = [ 1 , 1 ] because

a x = [ 1 , 4 ] [ 1 , 1 ] = [ 1 , 4 ]

If the uncertainty is even bigger i.e. a = [ 1 , 8 ] , then the solution doesn't exist. It is really hard to find physical interpretation of the algebraic interval solution set. Because of that in applications usually the united solution set is applied.

The method

Consider PDE with the interval parameters

( 1 )       G ( x , u , p ) = 0

where p = ( p 1 , … , p m ) ∈ p is a vector of parameters which belong to given intervals

p i ∈ [ p _ i , p ¯ i ] = p i , p = p 1 × p 2 × ⋯ × p m .

For example, the heat transfer equation

k x ∂ 2 u ∂ x 2 + k y ∂ 2 u ∂ y 2 + q = 0  for  x ∈ Ω u ( x ) = u ∗ ( x )  for  x ∈ ∂ Ω

where k x , k y are the interval parameters (i.e. k x ∈ k x ,   k y ∈ k y ).

Solution of the equation (1) can be defined in the following way

u ~ ( x ) := { u ( x ) : G ( x , u , p ) = 0 , p ∈ p }

For example, in the case of the heat transfer equation

u ~ ( x ) = { u ( x ) : k x ∂ 2 u ∂ x 2 + k y ∂ 2 u ∂ y 2 + q = 0  for  x ∈ Ω , u ( x ) = u ∗ ( x )  for  x ∈ ∂ Ω , k x ∈ k x ,   k y ∈ k y }

Solution u ~ is very complicated because of that in practice it is more interesting to find the smallest possible interval which contain the exact solution set u ~ .

u ( x ) = ◊ u ~ ( x ) = ◊ { u ( x ) : G ( x , u , p ) = 0 , p ∈ p }

For example, in the case of the heat transfer equation

u ( x ) = ◊ { u ( x ) : k x ∂ 2 u ∂ x 2 + k y ∂ 2 u ∂ y 2 + q = 0  for  x ∈ Ω , u ( x ) = u ∗ ( x )  for  x ∈ ∂ Ω , k x ∈ k x ,   k y ∈ k y }

Finite element method lead to the following parameter dependent system of algebraic equations

K ( p ) u = Q ( p ) ,       p ∈ p

where K is a stiffness matrix and Q is a right hand side.

Interval solution can be defined as a multivalued function

u = ◊ { u : K ( p ) u = Q ( p ) , p ∈ p }

In the simplest case above system can be treat as a system of linear interval equations.

It is also possible to define the interval solution as a solution of the following optimization problem

u _ i = min { u i : K ( p ) u = Q ( p ) , p ∈ p } u ¯ i = max { u i : K ( p ) u = Q ( p ) , p ∈ p }

In multidimensional case the intrval solution can be written as

u = u 1 × ⋯ × u n = [ u _ 1 , u ¯ 1 ] × ⋯ × [ u _ n , u ¯ n ]

History

Ben-Haim Y., Elishakoff I., 1990, Convex Models of Uncertainty in Applied Mechanics. Elsevier Science Publishers, New York

Valliappan S., Pham T.D., 1993, Fuzzy Finite Element Analysis of A Foundation on Elastic Soil Medium. International Journal for Numerical and Analytical Methods in Geomechanics, Vol.17, pp. 771–789

Elishakoff I., Li Y.W., Starnes J.H., 1994, A deterministic method to predict the effect of unknown-but-bounded elastic moduli on the buckling of composite structures. Computer methods in applied mechanics and engineering, Vol.111, pp. 155–167

Valliappan S. Pham T.D., 1995, Elasto-Plastic Finite Element Analysis with Fuzzy Parameters. International Journal for Numerical Methods in Engineering, 38, pp. 531–548

Rao S.S., Sawyer J.P., 1995, Fuzzy Finite Element Approach for the Analysis of Imprecisly Defined Systems. AIAA Journal, Vol.33, No.12, pp. 2364–2370

Köylüoglu H.U., Cakmak A., Nielsen S.R.K., 1995, Interval mapping in structural mechanics. In: Spanos, ed. Computational Stochastic Mechanics. 125-133. Balkema, Rotterdam

Muhanna, R. L. and R. L. Mullen (1995). "Development of Interval Based Methods for Fuzziness in Continuum Mechanics" in Proceedings of the 3rd International Symposium on Uncertainty Modeling and Analysis and Annual Conference of the North American Fuzzy Information Processing Society (ISUMA–NAFIPS'95),IEEE, 705–710

More references can be found here [4]

Interval solution versus probabilistic solution

It is important to know that the interval parameters generate different results than uniformly distributed random variables.

Interval parameter p = [ p _ , p ¯ ] take into account all possible probability distributions (for p ∈ [ p _ , p ¯ ] ).

In order to define the interval parameter it is necessary to know only upper p ¯ and lower bound p _ .

Calculations of probabilistic characteristics require the knowledge of a lot of experimental results.

It is possible to show that the sum of n interval numbers is n times wider than the sum of appropriate normally distributed random variables.

Sum of n interval number p = [ p _ , p ¯ ] is equal to

n p = [ n p _ , n p ¯ ]

Width of that interval is equal to

n p ¯ − n p _ = n ( p ¯ − p _ ) = n Δ p

Let us consider normally distributed random variable X such that

m X = E [ X ] = p ¯ + p _ 2 , σ X = V a r [ X ] = Δ p 6

Sum of n normally distributed random variable is a normally distributed random variable with the following characteristics (see Six Sigma)

E [ n X ] = n p ¯ + p _ 2 , σ n X = n V a r [ X ] = n σ = n Δ p 6

We can assume that the width of the probabilistic result is equal to 6 sigma (compare Six Sigma).

6 σ n X = 6 n Δ p 6 = n Δ p

Now we can compare the width of the interval result and the probabilistic result

w i d t h   o f   n   i n t e r v a l s w i d t h   o f   n   r a n d o m   v a r i a b l e s = n Δ p n Δ p = n

Because of that the results of the interval finite element (or in general worst case analysis) may be overestimated in comparison to the stochastic fem analysis (see also propagation of uncertainty). However, in the case of nonprobabilistic uncertainty it is not possible to apply pure probabilistic methods. Because probabilistic characteristic in that case are not known exactly [Elishakoff 2000].

It is possible to consider random (and fuzzy random variables) with the interval parameters (e.g. with the interval mean, variance etc.). Some researchers use interval (fuzzy) measurements in statistical calculations (e.g. [5]). As a results of such calculations we will get so called imprecise probability.

Imprecise probability is understood in a very wide sense. It is used as a generic term to cover all mathematical models which measure chance or uncertainty without sharp numerical probabilities. It includes both qualitative (comparative probability, partial preference orderings, …) and quantitative modes (interval probabilities, belief functions, upper and lower previsions, …). Imprecise probability models are needed in inference problems where the relevant information is scarce, vague or conflicting, and in decision problems where preferences may also be incomplete [6].

1-dimension example

In the tension-compression problem, the following equation shows the relationship between displacement u and force P :

E A L u = P

where L is length, A is the area of a cross-section, and E is Young's modulus.

If the Young's modulus and force are uncertain, then

E ∈ [ E _ , E ¯ ] , P ∈ [ P _ , P ¯ ]

To find upper and lower bounds of the displacement u , calculate the following partial derivatives:

∂ u ∂ E = − P L E 2 A < 0 ∂ u ∂ P = L E A > 0

Calculate extreme values of the displacement as follows:

u _ = u ( E ¯ , P _ ) = P _ L E ¯ A u ¯ = u ( E _ , P ¯ ) = P ¯ L E _ A

Calculate strain using following formula:

ε = 1 L u

Calculate derivative of the strain using derivative from the displacements:

∂ ε ∂ E = 1 L ∂ u ∂ E = − P E 2 A < 0 ∂ ε ∂ P = 1 L ∂ u ∂ P = 1 E A > 0

Calculate extreme values of the displacement as follows:

ε _ = ε ( E ¯ , P _ ) = P _ E ¯ A ε ¯ = ε ( E _ , P ¯ ) = P ¯ E _ A

It is also possible to calculate extreme values of strain using the displacements

∂ ε ∂ u = 1 L > 0

then

ε _ = ε ( u _ ) = P _ E ¯ A ε ¯ = ε ( u ¯ ) = P ¯ E _ A

The same methodology can be applied to the stress

σ = E ε

then

∂ σ ∂ E = ε + E ∂ ε ∂ E = ε + E 1 L ∂ u ∂ E = P E A − P E A = 0 ∂ σ ∂ P = E ∂ ε ∂ P = E 1 L ∂ u ∂ P = 1 A > 0

and

σ _ = σ ( P _ ) = P _ A σ ¯ = σ ( P ¯ ) = P ¯ A

If we treat stress as a function of strain then

∂ σ ∂ ε = ∂ ∂ ε ( E ε ) = E > 0

then

σ _ = σ ( ε _ ) = E ε _ = P _ A σ ¯ = σ ( ε ¯ ) = E ε ¯ = P ¯ A

Structure is safe if stress σ is smaller than a given value σ 0 i.e.

σ < σ 0

this condition is true if

σ ¯ < σ 0

After calculation we know that this relation is satisfied if

P ¯ A < σ 0

The example is very simple but it shows the applications of the interval parameters in mechanics. Interval FEM use very similar methodology in multidimensional cases [Pownuk 2004].

However, in the multidimensional cases relation between the uncertain parameters and the solution is not always monotone. In that cases more complicated optimization methods have to be applied [7].

Multidimensional example

In the case of tension-compression problem the equilibrium equation has the following form

d d x ( E A d u d x ) + n = 0

where u is displacement, E is Young's modulus, A is an area of cross-section, and n is a distributed load. In order to get unique solution it is necessary to add appropriate boundary conditions e.g.

u ( 0 ) = 0 d u ( 0 ) d x E A = P

If Young's modulus E and n are uncertain then the interval solution can be defined in the following way

u ( x ) = { u ( x ) : d d x ( E A d u d x ) + n = 0 , u ( 0 ) = 0 , d u ( 0 ) d x E A = P , E ∈ [ E _ , E ¯ ] , P ∈ [ P _ , P ¯ ] }

For each FEM element it is possible to multiply the equation by the test function v

∫ 0 L e ( d d x ( E A d u d x ) + n ) v = 0

where x ∈ [ 0 , L ( e ) ] .
After integration by parts we will get the equation in the weak form

∫ 0 L ( e ) E A d u d x d v d x d x = ∫ 0 L ( e ) n v d x

where x ∈ [ 0 , L ( e ) ] .
Let's introduce a set of grid points x 0 , x 1 , . . . , x N e , where N e is a number of elements, and linear shape functions for each FEM element

N 1 ( e ) ( x ) = 1 − 1 − x 0 ( e ) x 1 ( e ) − x 0 ( e ) ,     N 2 ( e ) ( x ) = 1 − x 0 ( e ) x 1 ( e ) − x 0 ( e ) .

where x ∈ [ x 0 ( e ) , x 1 ( e ) ] .
x 1 ( e ) left endpoint of the element, x 1 ( e ) left endpoint of the element number "e". Approximate solution in the "e"-th element is a linear combination of the shape functions

u h ( e ) ( x ) = u 1 e N 1 ( e ) ( x ) + u 2 e N 2 ( e ) ( x ) ,     v h ( e ) ( x ) = u 1 e N 1 ( e ) ( x ) + u 2 e N 2 ( e ) ( x )

After substitution to the weak form of the equation we will get the following system of equations

[ E ( e ) A ( e ) L ( e ) − E ( e ) A ( e ) L ( e ) − E ( e ) A ( e ) L ( e ) E ( e ) A ( e ) L ( e ) ] [ u 1 ( e ) u 2 ( e ) ] = [ ∫ 0 L ( e ) n N 1 ( e ) ( x ) d x ∫ 0 L ( e ) n N 2 ( e ) ( x ) d x ]

or in the matrix form

K ( e ) u ( e ) = Q ( e )

In order to assemble the global stiffness matrix it is necessary to consider an equilibrium equations in each node. After that the equation has the following matrix form

K u = Q

where

K = [ K 11 ( 1 ) K 12 ( 1 ) 0 . . . 0 K 21 ( 1 ) K 22 ( 1 ) + K 11 ( 2 ) K 12 ( 2 ) . . . 0 0 K 21 ( 2 ) K 22 ( 2 ) + K 11 ( 3 ) . . . 0 . . . . . . . . . . . . . . . 0 0 . . . K 22 ( N e − 1 ) + K 11 ( N e ) K 11 ( N e ) 0 0 . . . K 21 ( N e ) K 22 ( N e ) ]

is the global stiffness matrix,

u = [ u 0 u 1 . . . u N e ]

is the solution vector,

Q = [ Q 0 Q 1 . . . Q N e ]

is the right hand side.

In the case of tension-compression problem

K = [ E ( 1 ) A ( 1 ) L ( 1 ) − E ( 1 ) A ( 1 ) L ( 1 ) 0 . . . 0 − E ( 1 ) A ( 1 ) L ( 1 ) E ( 1 ) A ( 1 ) L ( 1 ) + E ( 2 ) A ( 2 ) L ( 2 ) − E ( 2 ) A ( 2 ) L ( 2 ) . . . 0 0 − E ( 2 ) A ( 2 ) L ( 2 ) E ( 2 ) A ( 2 ) L ( 2 ) + E ( 3 ) A ( 3 ) L ( 3 ) . . . 0 . . . . . . . . . . . . . . . 0 0 . . . E ( N e − 1 ) A ( N e − 1 ) L ( N e − 1 ) + E ( N e ) A ( N e ) L ( N e ) − E ( N e ) A ( N e ) L ( N e ) 0 0 . . . − E ( N e ) A ( N e ) L ( N e ) E ( N e ) A ( N e ) L ( N e ) ]

If we neglect the distributed load n

Q = [ R 0 . . . 0 P ]

After taking into account the boundary conditions the stiffness matrix has the following form

K = [ 1 0 0 . . . 0 0 E ( 1 ) A ( 1 ) L ( 1 ) + E ( 2 ) A ( 2 ) L ( 2 ) − E ( 2 ) A ( 2 ) L ( 2 ) . . . 0 0 − E ( 2 ) A ( 2 ) L ( 2 ) E ( 2 ) A ( 2 ) L ( 2 ) + E ( 3 ) A ( 3 ) L ( 3 ) . . . 0 . . . . . . . . . . . . . . . 0 0 . . . E ( e − 1 ) A ( e − 1 ) L ( e − 1 ) + E ( e ) A ( e ) L ( e ) − E ( e ) A ( e ) L ( e ) 0 0 . . . − E ( e ) A ( e ) L ( e ) E ( e ) A ( e ) L ( e ) ] = K ( E , A ) = K ( E ( 1 ) , . . . , E ( N e ) , A ( 1 ) , . . . , A ( N e ) )

Right-hand side has the following form

Q = [ 0 0 . . . 0 P ] = Q ( P )

Let's assume that Young's modulus E , area of cross-section A and the load P are uncertain and belong to some intervals

E ( e ) ∈ [ E _ ( e ) , E ¯ ( e ) ] A ( e ) ∈ [ A _ ( e ) , A ¯ ( e ) ] P ∈ [ P _ , P ¯ ]

The interval solution can be defined calculating the following way

u = ◊ { u : K ( E , A ) u = Q ( P ) , E ( e ) ∈ [ E _ ( e ) , E ¯ ( e ) ] , A ( e ) ∈ [ A _ ( e ) , A ¯ ( e ) ] , P ∈ [ P _ , P ¯ ] }

Calculation of the interval vector u is in general NP-hard, however in specific cases it is possible to calculate the solution which can be used in many engineering applications.

The results of the calculations are the interval displacements

u i ∈ [ u _ i , u ¯ i ]

Let's assume that the displacements in the column have to be smaller than some given value (due to safety).

u i < u i m a x

The uncertain system is safe if the interval solution satisfy all safety conditions.

In this particular case

u i < u i m a x ,       u i ∈ [ u _ i , u ¯ i ]

or simple

u ¯ i < u i m a x

In postprocessing it is possible to calculate the interval stress, the interval strain and the interval limit state functions and use these values in the design process.

The interval finite element method can be applied to the solution of problems in which there is not enough information to create reliable probabilistic characteristic of the structures [Elishakoff 2000]. Interval finite element method can be also applied in the theory of imprecise probability.

Endpoints combination method

It is possible to solve the equation K ( p ) u ( p ) = Q ( p ) for all possible combinations of endpoints of the interval p ^ .
The list of all vertices of the interval p ^ can be written as L = { p 1 ∗ , . . . , p n ∗ } .
Upper and lower bound of the solution can be calculated in the following way

u _ i = m i n { u i ( p k ∗ ) : K ( p k ∗ ) u ( p k ∗ ) = Q ( p k ∗ ) , p k ∗ ∈ L } u ¯ i = m a x { u i ( p k ∗ ) : K ( p k ∗ ) u ( p k ∗ ) = Q ( p k ∗ ) , p k ∗ ∈ L }

Endpoints combination method gives solution which is usually exact; unfortunately the method has exponential computational complexity and cannot be applied to the problems with many interval parameters [Neumaier 1990].

Taylor expansion method

The function u = u ( p ) can be expanded by using Taylor series. In the simplest case the Taylor series use only linear approximation

u i ( p ) ≈ u i ( p 0 ) + ∑ j ∂ u ( p 0 ) ∂ p j Δ p j

Upper and lower bound of the solution can be calculated by using the following formula

u _ i ≈ u i ( p 0 ) − | ∑ j ∂ u ( p 0 ) ∂ p j | Δ p j u ¯ i ≈ u i ( p 0 ) + | ∑ j ∂ u ( p 0 ) ∂ p j | Δ p j

The method is very efficient however it is not very accurate.
In order to improve accuracy it is possible to apply higher order Taylor expansion [Pownuk 2004].
This approach can be also applied in the interval finite difference method and the interval boundary element method.

Gradient method

If the sign of the derivatives ∂ u i ∂ p j is constant then the functions u i = u i ( p ) is monotone and the exact solution can be calculated very fast.

if ∂ u i ∂ p j ≥ 0 then p i m i n = p _ i ,   p i m a x = p ¯ i if ∂ u i ∂ p j < 0 then p i m i n = p ¯ i ,   p i m a x = p _ i

Extreme values of the solution can be calculated in the following way

u _ i = u i ( p m i n ) ,   u ¯ i = u i ( p m a x )

In many structural engineering applications the method gives exact solution.
If the solution is not monotone the solution is usually reasonable. In order to improve accuracy of the method it is possible to apply monotonicity tests and higher order sensitivity analysis. The method can be applied to the solution of linear and nonlinear problems of computational mechanics [Pownuk 2004]. Applications of sensitivity analysis method to the solution of civil engineering problems can be found in the following paper [M.V. Rama Rao, A. Pownuk and I. Skalna 2008].
This approach can be also applied in the interval finite difference method and the interval boundary element method.

Element by element method

Muhanna and Mullen applied element by element formulation to the solution of finite element equation with the interval parameters [Muhanna, Mullen 2001]. Using that method it is possible to get the solution with guaranteed accuracy in the case of truss and frame structures.

Perturbation methods

The solution u = u ( p ) stiffness matrix K = K ( p ) and the load vector Q = Q ( p ) can be expanded by using perturbation theory. Perturbation theory lead to the approximate value of the interval solution [Qiu, Elishakoff 1998]. The method is very efficient and can be applied to large problems of computational mechanics.

Response surface method

It is possible to approximate the solution u = u ( p ) by using response surface. Then it is possible to use the response surface to the get the interval solution [Akpan 2000]. Using response surface method it is possible to solve very complex problem of computational mechanics [Beer 2008].

Pure interval methods

Several authors tried to apply pure interval methods to the solution of finite element problems with the interval parameters. In some cases it is possible to get very interesting results e.g. [Popova, Iankov, Bonev 2008]. However, in general the method generates very overestimated results [Kulpa, Pownuk, Skalna 1998].

Parametric interval systems

[Popova 2001] and [Skalna 2006] introduced the methods for the solution of the system of linear equations in which the coefficients are linear combinations of interval parameters. In this case it is possible to get very accurate solution of the interval equations with guaranteed accuracy.