Numerical Integration

From KratosWiki
(Difference between revisions)
Jump to: navigation, search
(→Gauss-Legendre Numerical Integration)
(→Numerical Integration for Isoparametric Triangular Domains)
 
(30 intermediate revisions by one user not shown)
Line 27: Line 27:
  
  
−
{| border="1" class="wikitable" style="margin:auto; background:lightgray;"
+
{| border="1" cellpadding="5" cellspacing="0" class="wikitable" style="margin:auto; background:white;"
−
! Number of points, ''p'' !! Points, ''&xi;''<sub>''i'' !! Weights, ''w''<sub>''i''</sub>
+
! Number of points, ''p'' !! Points, ''&plusmn;&xi;''<sub>''i'' !! Weights, ''w''<sub>''i''</sub>
 
|- align="center"
 
|- align="center"
−
| 1 || 0 || 2
+
| <math>1\,</math> || <math>0.0 \,</math> || <math>2.0\,</math>
 
|- align="center"
 
|- align="center"
−
| 2 || <math>\pm\sqrt{1/3}</math> || 1
+
| <math>2\,</math> || <math>\pm\sqrt{1/3}</math> || <math>1.0\,</math>
 
|- align="center"
 
|- align="center"
−
| rowspan="2" | 3 || 0 || <sup>8</sup>⁄<sub>9</sub>
+
| rowspan="2" | <math>3\,</math> || <math>0.0 \,</math> || <math>\frac{8}{9}</math>
 
|- align="center"
 
|- align="center"
−
| <math>\pm\sqrt{3/5}</math> || <sup>5</sup>⁄<sub>9</sub>
+
| <math>\pm\sqrt{3/5}</math> || <math>\frac{5}{9}</math>
 
|- align="center"
 
|- align="center"
−
| rowspan="2" | 4 || <math>\pm\sqrt{\Big( 3 - 2\sqrt{6/5} \Big)/7}</math> || <math>\tfrac{18+\sqrt{30}}{36}</math>
+
| rowspan="2" | <math>4\,</math> || <math>\pm\sqrt{\Big( 3 - 2\sqrt{6/5} \Big)/7}</math> || <math>\tfrac{18+\sqrt{30}}{36}</math>
 
|- align="center"
 
|- align="center"
 
| <math>\pm\sqrt{\Big( 3 + 2\sqrt{6/5} \Big)/7}</math> || <math>\tfrac{18-\sqrt{30}}{36}</math>
 
| <math>\pm\sqrt{\Big( 3 + 2\sqrt{6/5} \Big)/7}</math> || <math>\tfrac{18-\sqrt{30}}{36}</math>
 
|- align="center"
 
|- align="center"
−
| rowspan="3" | 5 || 0 || <sup>128</sup>⁄<sub>225</sub>
+
| rowspan="3" | <math>5\,</math> || <math>0.0 \,</math> || <math>\frac{128}{225}</math>
 
|- align="center"
 
|- align="center"
 
| <math>\pm\tfrac13\sqrt{5-2\sqrt{10/7}}</math> || <math>\tfrac{322+13\sqrt{70}}{900}</math>
 
| <math>\pm\tfrac13\sqrt{5-2\sqrt{10/7}}</math> || <math>\tfrac{322+13\sqrt{70}}{900}</math>
Line 48: Line 48:
 
| <math>\pm\tfrac13\sqrt{5+2\sqrt{10/7}}</math> || <math>\tfrac{322-13\sqrt{70}}{900}</math>
 
| <math>\pm\tfrac13\sqrt{5+2\sqrt{10/7}}</math> || <math>\tfrac{322-13\sqrt{70}}{900}</math>
 
|}
 
|}
 +
 +
 +
or, using numerical values:
 +
 +
 +
{| border="1" cellpadding="5" cellspacing="0" class="wikitable" style="margin:auto; background:white;"
 +
! Number of points, ''p'' !! Points, ''&plusmn;&xi;''<sub>''i'' !! Weights, ''w''<sub>''i''</sub>
 +
|-
 +
| 1 || 0.0 || 2.0
 +
|-
 +
| 2 || 0.5773502692 || 1.0
 +
|-
 +
| rowspan="2" | 3 || 0.0 || 0.8888888889
 +
|-
 +
| 0.774596697 || 0.5555555556
 +
|-
 +
| rowspan="2" | 4 || 0.3399810436 || 0.6521451549
 +
|-
 +
| 0.8611363116 || 0.3478548451
 +
|-
 +
| rowspan="3" | 5 || 0.0 || 0.5688888889
 +
|-
 +
| 0.5384693101 || 0.4786286705
 +
|-
 +
| 0.9061798459 || 0.2369268851
 +
|-
 +
| rowspan="3" | 6 || 0.2386191861 || 0.4679139346
 +
|-
 +
| 0.6612093865 || 0.3607615730
 +
|-
 +
| 0.9324695142 || 0.1713244924
 +
|-
 +
| rowspan="4" | 7 || 0.0 || 0.4179591837
 +
|-
 +
| 0.4058451514 || 0.3818300505
 +
|-
 +
| 0.7415311856 || 0.2797053915
 +
|-
 +
| 0.9491079123 || 0.1294849662
 +
|-
 +
| rowspan="4" | 8 || 0.1834346425 || 0.3626837834
 +
|-
 +
| 0.5255324099 || 0.3137066459
 +
|-
 +
| 0.7966664774 || 0.2223810345
 +
|-
 +
| 0.9602898565 || 0.1012285636
 +
|}
 +
 +
 +
=== Example of a one dimensional integration ===
 +
 +
 +
For the function: <math>f(x)=1+x+x^2+x^3+x^4 \,</math>, the exact integration in [-1,+1] is:
 +
 +
:<math>I=\int_{-1}^{+1} f(x) dx = \left . \left ( x + \frac{x^2}{2} + \frac{x^3}{3} + \frac{x^4}{4} + \frac{x^5}{5}  \right ) \right |_{-1}^{+1} = 2 + 2 \frac{1}{3} + 2 \frac{1}{5} = 3.0666</math>
 +
 +
Numerically:
 +
 +
:'''First order Gauss-Legendre Quadrature:'''
 +
 +
::<math>p=1, x_1=0, W_1=2; \qquad I=W_1 f(x_1)=2</math>
 +
 +
 +
:'''Second order Gauss-Legendre Quadrature:'''
 +
 +
::<math>p=2
 +
\begin{cases}
 +
    x_1 = - 0.57735, & W_1 = 1 \\
 +
    x_2 = + 0.57735, & W_2 = 1
 +
\end{cases} \qquad I=W_1 f(x_1) + W_2 f(x_2) = 0.67464 + 2.21424 = 2.8888</math>
 +
 +
 +
:'''Third order Gauss-Legendre Quadrature:'''
 +
 +
::<math>p=3
 +
\begin{cases}
 +
    x_1 = - 0.77459, & W_1 = 0.5555 \\
 +
    x_2 = - 0.00000, & W_2 = 0.8888 \\
 +
    x_3 = + 0.77459, & W_3 = 0.5555
 +
\end{cases}
 +
</math>
 +
 +
::<math>I=W_1 f(x_1) + W_2 f(x_2) + W_3 f(x_3) = 0.7204·0.5555 + 1.0·0.8888 + 3.19931·0.5555 = 3.0666 \,</math>
 +
 +
 +
:That is the exact value, because for any polynomial function of '''p<sup>th</sup>''' order it is enough to use '''p-1''' integration points.
 +
 +
 +
== Two Dimensional Numerical Integration ==
 +
 +
By using isoparametric formulation we can use natural coordinates to compute any integration. In addition, we can still use the Gauss-Legendre quadrature.
 +
 +
=== Numerical Integration for Isoparametric Triangular Domains ===
 +
 +
 +
A general integral expression form for two dimensional domains can be written in terms of the [[Two-dimensional_Shape_Functions#Areal_Coordinates | '''area coordinates''']] and, therefore, computed by using the Gauss quadrature:
 +
 +
 +
::<math>\int_0^1 \int_0^{1-L_3} f(L_1,L_2,L_3) dL_2 dL_3 = \sum_{p=1}^{n_p} f(L_{1_p},L_{2_p},L_{3_p}) W_p</math>
 +
 +
with:
 +
* <math>n_p \,</math> &nbsp; the number of integration points;
 +
* <math>L_{1_p}, L_{2_p}, L_{3_p} \,</math> &nbsp; the value of the area coordinates;
 +
* <math>W_p \,</math> &nbsp; the weight in the integration point '''''p''''';
 +
 +
The following table and picture shows the integration points and weights for triangles obtained from the Gaussian quadrature ('''precision''' means the degree of polynomial for exact integration):
 +
 +
 +
{| border="1" cellpadding="5" cellspacing="0" class="wikitable" style="margin:auto; background:white;"
 +
! Number of points, ''n'' !! ''precision'' !! Points !! '''L<sub>1</sub>''' !! '''L<sub>2</sub>''' !! '''L<sub>3</sub>''' !! ''W''<sub>''i''</sub>
 +
|- align="center"
 +
| 1 || Linear || ''a'' || 1/3 || 1/3 || 1/3 || 1/2
 +
|- align="center"
 +
| rowspan="3" | 3 || rowspan="3" | Quadratic || ''a'' || 1/2 || 1/2 || 0 || 1/6
 +
|- align="center"
 +
| ''b'' || 0 || 1/2 || 1/2 || 1/6
 +
|- align="center"
 +
| ''c'' || 1/2 || 0 || 1/2 || 1/6
 +
|- align="center"
 +
| rowspan="4" | 4 || rowspan="4" | Cubic || ''a'' || 1/3 || 1/3 || 1/3 || -9/32
 +
|- align="center"
 +
| ''b'' || 0.6 || 0.2 || 0.2 || 25/96
 +
|- align="center"
 +
| ''c'' || 0.2 || 0.6 || 0.2 || 25/96
 +
|- align="center"
 +
| ''d'' || 0.2 || 0.2 || 0.6 || 25/96
 +
|- align="center"
 +
| rowspan="7" | 7 || rowspan="7" | Quartic || ''a'' || 0 || 0 || 1 || 1/40
 +
|- align="center"
 +
| ''b'' || 1/2 || 0 || 1/2 || 1/15
 +
|- align="center"
 +
| ''c'' || 1 || 0 || 0 || 1/40
 +
|- align="center"
 +
| ''d'' || 1/2 || 1/2 || 0 || 1/15
 +
|- align="center"
 +
| ''e'' || 0 || 1 || 0 || 1/40
 +
|- align="center"
 +
| ''f'' || 0 || 1/2 || 1/2 || 1/15
 +
|- align="center"
 +
| ''g'' || 1/3 || 1/3 || 1/3 || 9/40
 +
|}
 +
 +
 +
::[[Image:IntegrationPointsTriangularElement.jpg]]
 +
 +
 +
Note that the weight values has been normalised in order to sum 1/2 to maintain the exact value for the element area.
 +
 +
 +
Therefore:
 +
 +
:<math>A^{(e)} = \int \int_{A^{(e)}} dx dy = \int_0^1 \int_0^{1-\beta} |J^{(e)}| d\alpha d\beta = |J^{(e)}| \sum_p W_p = \frac{|J^{(e)}|}{2}</math>
  
 
== References ==
 
== References ==

Latest revision as of 18:41, 11 November 2009

Numerical integration refers to all the procedures, algorithms and techniques in the numerical analysis to obtain an approximate solution to a definite integral.

That is, how to obtain a numerical value of:

\int_{\lambda_a}^{\lambda_b}\! f(\lambda)\, d\lambda.

where \lambda \, can be a 1D, 2D or 3D domain.


For our interest in the Finite Element Method, the purpose is to describe how the element matrices can be integrated numerically.


Contents

Gauss-Legendre Numerical Integration

To fix the most basic concepts on Numerical Integration, we will focus our description on a one dimensional integration using the Gauss-Legendre quadrature, that is, to solve:

I=\int_{-1}^{+1} f(\xi) d\xi


The Gauss-Legendre quadrature establish that the definite integral of a function can be approximate by using a weighted sum of function values at specified points within the domain of integration. An p-point Gaussian quadrature rule is constructed to yield an exact result for polynomials of degree 2p − 1 or less by a suitable choice of the points   \xi_i \,   and weights   w_i \,   for   i = 1, \cdots p \,.


\int_{-1}^{+1} f(\xi)\,d\xi \approx \sum_{i=1}^p w_i f(\xi_i)


The coordinates and related weights are:


Number of points, p Points, ±ξi Weights, wi
1\, 0.0 \, 2.0\,
2\, \pm\sqrt{1/3} 1.0\,
3\, 0.0 \, \frac{8}{9}
\pm\sqrt{3/5} \frac{5}{9}
4\, \pm\sqrt{\Big( 3 - 2\sqrt{6/5} \Big)/7} \tfrac{18+\sqrt{30}}{36}
\pm\sqrt{\Big( 3 + 2\sqrt{6/5} \Big)/7} \tfrac{18-\sqrt{30}}{36}
5\, 0.0 \, \frac{128}{225}
\pm\tfrac13\sqrt{5-2\sqrt{10/7}} \tfrac{322+13\sqrt{70}}{900}
\pm\tfrac13\sqrt{5+2\sqrt{10/7}} \tfrac{322-13\sqrt{70}}{900}


or, using numerical values:


Number of points, p Points, ±ξi Weights, wi
1 0.0 2.0
2 0.5773502692 1.0
3 0.0 0.8888888889
0.774596697 0.5555555556
4 0.3399810436 0.6521451549
0.8611363116 0.3478548451
5 0.0 0.5688888889
0.5384693101 0.4786286705
0.9061798459 0.2369268851
6 0.2386191861 0.4679139346
0.6612093865 0.3607615730
0.9324695142 0.1713244924
7 0.0 0.4179591837
0.4058451514 0.3818300505
0.7415311856 0.2797053915
0.9491079123 0.1294849662
8 0.1834346425 0.3626837834
0.5255324099 0.3137066459
0.7966664774 0.2223810345
0.9602898565 0.1012285636


Example of a one dimensional integration

For the function: f(x)=1+x+x^2+x^3+x^4 \,, the exact integration in [-1,+1] is:

I=\int_{-1}^{+1} f(x) dx = \left . \left ( x + \frac{x^2}{2} + \frac{x^3}{3} + \frac{x^4}{4} + \frac{x^5}{5} \right ) \right |_{-1}^{+1} = 2 + 2 \frac{1}{3} + 2 \frac{1}{5} = 3.0666

Numerically:

First order Gauss-Legendre Quadrature:
p=1, x_1=0, W_1=2; \qquad I=W_1 f(x_1)=2


Second order Gauss-Legendre Quadrature:
p=2 \begin{cases} x_1 = - 0.57735, & W_1 = 1 \\ x_2 = + 0.57735, & W_2 = 1 \end{cases} \qquad I=W_1 f(x_1) + W_2 f(x_2) = 0.67464 + 2.21424 = 2.8888


Third order Gauss-Legendre Quadrature:
p=3 \begin{cases} x_1 = - 0.77459, & W_1 = 0.5555 \\ x_2 = - 0.00000, & W_2 = 0.8888 \\ x_3 = + 0.77459, & W_3 = 0.5555 \end{cases}
I=W_1 f(x_1) + W_2 f(x_2) + W_3 f(x_3) = 0.7204·0.5555 + 1.0·0.8888 + 3.19931·0.5555 = 3.0666 \,


That is the exact value, because for any polynomial function of pth order it is enough to use p-1 integration points.


Two Dimensional Numerical Integration

By using isoparametric formulation we can use natural coordinates to compute any integration. In addition, we can still use the Gauss-Legendre quadrature.

Numerical Integration for Isoparametric Triangular Domains

A general integral expression form for two dimensional domains can be written in terms of the area coordinates and, therefore, computed by using the Gauss quadrature:


\int_0^1 \int_0^{1-L_3} f(L_1,L_2,L_3) dL_2 dL_3 = \sum_{p=1}^{n_p} f(L_{1_p},L_{2_p},L_{3_p}) W_p

with:

  • n_p \,   the number of integration points;
  • L_{1_p}, L_{2_p}, L_{3_p} \,   the value of the area coordinates;
  • W_p \,   the weight in the integration point p;

The following table and picture shows the integration points and weights for triangles obtained from the Gaussian quadrature (precision means the degree of polynomial for exact integration):


Number of points, n precision Points L1 L2 L3 Wi
1 Linear a 1/3 1/3 1/3 1/2
3 Quadratic a 1/2 1/2 0 1/6
b 0 1/2 1/2 1/6
c 1/2 0 1/2 1/6
4 Cubic a 1/3 1/3 1/3 -9/32
b 0.6 0.2 0.2 25/96
c 0.2 0.6 0.2 25/96
d 0.2 0.2 0.6 25/96
7 Quartic a 0 0 1 1/40
b 1/2 0 1/2 1/15
c 1 0 0 1/40
d 1/2 1/2 0 1/15
e 0 1 0 1/40
f 0 1/2 1/2 1/15
g 1/3 1/3 1/3 9/40


IntegrationPointsTriangularElement.jpg


Note that the weight values has been normalised in order to sum 1/2 to maintain the exact value for the element area.


Therefore:

A^{(e)} = \int \int_{A^{(e)}} dx dy = \int_0^1 \int_0^{1-\beta} |J^{(e)}| d\alpha d\beta = |J^{(e)}| \sum_p W_p = \frac{|J^{(e)}|}{2}

References

  1. Carlos A. Felippa, "A compendium of FEM integration formulas for symbolic work", Engineering Computations, Vol. 21 No. 8, 2004, pp. 867-890, (c) Emerald Group Publishing Limited [1]
  2. Numerical Integration
  3. Gaussian Quadrature

Retrieved from "https://kratos-wiki.cimne.upc.edu/index.php/Numerical_Integration"
Personal tools
Categories
Theory Manual