```python
import theme
theme.load_style()
```




<style>
@import url('http://fonts.googleapis.com/css?family=Source+Code+Pro');
@import url('http://fonts.googleapis.com/css?family=Vollkorn');
@import url('http://fonts.googleapis.com/css?family=Arimo');
@import url('http://fonts.googleapis.com/css?family=Arvo:400,700,400italic')

#notebook li { /* More space between bullet points */
    margin-top:0.8em;
}

/* draw border around running cells */
div.cell.border-box-sizing.code_cell.running {
    border: 3px solid #111;
}

div.cell.code_cell {
    margin-top: 1em;
    margin-left: auto;
    margin-right: auto;
    width: 800px;
    border: 1px solid rgba(200, 200, 200, 0.8);
    border-radius: 8px; /* rounded borders */
}

div.text_cell_render{
    font-family: 'Arvo' sans-serif;
    /* font-family: 'Arimo',verdana,arial,sans-serif; */
    line-height: 130%;
    font-size: 115%;
    margin-left:0%;
    margin-right:auto;
}

.text_cell_render h1 {
    font-family: 'Open sans',verdana,arial,sans-serif;
    font-weight: 300;
    font-style: bold;
    font-size: 30pt;
    line-height: 100%;
    color:#c76c0c;
    margin-top: 0.75em;
    margin-bottom: 0.05em;
    display: block;
}

.text_cell_render h2 {
    font-family: 'Open sans',verdana,arial,sans-serif;
    font-weight: 200;
    font-size: 20pt;
    line-height: 100%;
    color:#c76c0c;
    margin-top: 0.75em;
    margin-bottom: 0.05em;
    display: block;
}

.text_cell_render h3 {
    font-family: 'Open sans',verdana,arial,sans-serif;
    font-weight: 200;
    font-size: 16pt;
    line-height: 100%;
    color:#d77c0c;
    margin-top: 0.75em;
    margin-bottom: 0.05em;
    display: block;
}

.text_cell_render h4 {
    font-family: 'Open sans',verdana,arial,sans-serif;
    font-weight: 200;
    font-size: 14pt;
    font-style: italic;
    line-height: 100%;
    color:#d77c0c;
    margin-top: 0.75em;
    margin-bottom: 0.05em;
    display: block;
}

.text_cell_render h5 {
    font-family: 'Open sans',verdana,arial,sans-serif;
    font-weight: 200;
    font-size: 12pt;
    color: #1d3b84;
    margin-top: 0.75em;
    margin-bottom: 0.05em;
}

.rendered_html code{
    background-color: transparent;
}

ul{
    margin: 2em;
}

ul li{
    padding-left: 0.5em;
    margin-bottom: 0.5em;
    margin-top: 0.5em;
}

ul li li{
    padding-left: 0.2em;
    margin-bottom: 0.2em;
    margin-top: 0.2em;
}

ol{
    margin: 2em;
}

ol li{
    padding-left: 0.5em;
    margin-bottom: 0.5em;
    margin-top: 0.5em;
}

ul li{
    padding-left: 0.5em;
    margin-bottom: 0.5em;
    margin-top: 0.2em;
}

a:link{
    font-weight: bold;
    color:#447adb;
}

a:visited{
    font-weight: bold;
    color: #1d3b84;
}

a:hover{
    font-weight: bold;
    color: #1d3b84;
}

a:focus{
    font-weight: bold;
    color:#447adb;
}

a:active{
    font-weight: bold;
    color:#447adb;
}

.rendered_html :link {
    text-decoration: underline;
}

.rendered_html :hover {
    text-decoration: none;
}

.rendered_html :visited {
    text-decoration: none;
}

.rendered_html :focus {
    text-decoration: none;
}

.rendered_html :active {
    text-decoration: none;
}

.warning{
    color: rgb( 240, 20, 20 )
}

hr {
    color: #f3f3f3;
    background-color: #f3f3f3;
    height: 1px;
}

blockquote{
    display:block;
    background: #fcfcfc;
    border-left: 5px solid #c76c0c;
    font-family: 'Open sans',verdana,arial,sans-serif;
    width:680px;
    padding: 10px 10px 10px 10px;
}

blockquote p {
    margin-bottom: 0;
    line-height: 125%;
    font-size: 100%;
}

table#mytable,
table#mytable td,
table#mytable tr
{
    border: none !important;
}

table#mytable-1,
table#mytable-1 td,
table#mytable-1 tr
{
    border: none !important;
    align: center;
    cellspacing: 10;
}

.msg {
    width: 90%;
    margin-left: auto;
    margin-right: auto;
    margin-top: 10px;
    margin-bottom: 10px;
    padding: 14px 14px 14px 14px;
    color: #c09853;
    text-shadow: 0 1px 0 rgba(255, 255, 255, 0.5);
    background-color: #fcf8e3;
    border: 1px solid #fbeed5;
    -webkit-border-radius: 4px;
    -moz-border-radius: 4px;
    border-radius: 4px;
}

.msg-heading {
    color: inherit;
}

.msg .close {
    position: relative;
    top: -2px;
    right: -21px;
    line-height: 18px;
}

.msg-success {
    color: #468847;
    background-color: #dff0d8;
    border-color: #d6e9c6;
}

.msg-danger,
.msg-error {
    color: #b94a48;
    background-color: #f2dede;
    border-color: #eed3d7;
}

.msg-info {
    color: #3a87ad;
    background-color: #d9edf7;
    border-color: #bce8f1;
}

.boxed {
    color: #b94a48;
    background-color: #f2dede;
    border-color: #eed3d7;
    margin-left: auto;
    margin-right: auto;
    padding: 14px 14px 14px 14px;
}

.w-300 { width: 300px; }
.w-400 { width: 400px; }
.w-500 { width: 500px; }

</style>




# Lesson 12: Numerical Integration

<a rel="license" href="http://creativecommons.org/licenses/by/4.0/"></a>

This lecture by Tim Fuller is licensed under the
<a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.  All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT).

<a id='top'></a>
# Topics

- [Numerical Integration](#num_int)  
- [Gauss Legendre Quadrature](#gl_quad)
- [Natural Coordinates](#nat_coords)
  - [Approximation of Geometry](#appx_geom)
  - [\*Parametric Forms](#parametric_forms)
- [Examples](#examples)
  - [Example 1: Element Stiffness Integration](#ex-1)
  - [Example 2: Reduced Integration](#ex-2)

# <a id='num_int'> </a> Numerical Integration[](#top)

Consider the element stiffness

$$
K_{ij}^{(e)} = \int_{x_1^{(e)}}^{x_2^{(e)}} a^{(e)} \phi'_i \phi'_j dx
$$

In the event that $a_e$ is a complicated function of $x$, the exact evaluation of the stiffness may not be possible.  In such cases, numerical integration is necessary.  Likewise, with other integrals encountered in the FEM.

Numerical evaluation of integrals is called numerical quadrature and typically involves approximating the integrand by a polynomial of sufficient degree.

For example, consider the integral

$$
J[f] = \int_{x_A}^{x_B} f(x) dx
$$

We approximate $f(x)$ by

$$
f(x) \approx \sum_i^n f(x_i) \phi_i(x)
$$

where the $x_i$ are the $i^{\rm th}$ point of the interval $[x_A, x_B]$ and $\phi_i(x_i)$ are polynomials of degree $n-1$.  Substitution of the approximate solution for $f$ in to the integral equation gives

$$
\begin{align}
J[f] &= \int_{x_A}^{x_B} f(x) dx \\
 &\approx \sum_{i=1}^n f(x_i) \int_{x_A}^{x_B} \phi_i(x) dx \\
& = f(x_1) \int_{x_A}^{x_B} \phi_1(x) dx + \ldots + f(x_N) \int_{x_A}^{x_B} \phi_n(x) dx
\end{align}
$$

Suppose we choose Lagrange linear interpolation polynomials for the approximation of $f$ with $n=2$.  Then

$$
J[f] \approx f(x_A) \int_{x_A}^{x_B} \frac{x_B - x}{h} dx + f(x_B) \int_{x_A}^{x_B} \frac{x - x_A}{h}  dx = \frac{1}{2}h\left(f(x_A) + f(x_B)\right)
$$

Graphically, this represents the area of a trapezoid used to approximate the area under the function $f(x)$



This is known as the trapzoidal rule.  If we use quadratic Lagrange interpolation of $f(x)$ we obtain

$$
J[f] \approx \frac{1}{6}h \left(f(x_A) + 4 f(x_A + \frac{1}{2}h) + f(x_B)\right)
$$

The previous examples demonstrate the basic idea behind numerical quadrature.  In the trapezoidal type rules, we choose the $x_i$ and the weights shake out naturally.  It can be shown that this form of numerical integration is exact for polynomials of up to degree $n-1$.  The trapezoidal rules also demonstrate a more general form of numerical quadrature:

$$
\begin{align}
J[f] &= \int_{x_A}^{x_B} f(x) w(x) dx \\
J_n[f] &= \int_{x_A}^{x_B}w(x) \sum_{i=1}^n f(x_i) L_i(x) 
= \sum_{i=1}^n w_i f(x_i)
\end{align}
$$

with

$$
w_i = \int_{x_A}^{x_B}L_i(x)w(x) dx, \quad 
L_i^n(x)=\prod_{\stackrel{j=1}{j \neq i}}^n \frac{x-x_j}{x_j-x_i}
$$

where the $x_i$ are called quadrature points, $w_i$ are quadrature weights, and $L_i(x)$ are $n^{\text{th}}$ order Lagrange polynomials.  The goal of numerical quadrature is to choose the $x_i$ and the $w_i$ such that $J_n$ is exact for polynomials up to degree $2n-1$.

# <a id='gl_quad'></a> Gauss-Legendre Quadrature[](#top)

So, how do we choose the $x_i$ and $w_i$?  That question is the subject of numerical analysis, and the exact details are beyond the scope of this class.  Here, we seek to determine the quadrature points and weights for small $n$ and appeal to reason and mathematical proofs to furnish the results for all $n$.  To that end, we consider only the special case that $w(x)=1$, $a=-1$, $b=1$, and $f(x) \in P_n$, where $P_n$ is the space of polynomials of degree $n$.  This quadrature rule is known as Gauss-Legendre quadrature.  Other quadrature rules are distinguished by their choice of $w$, $a$, and $b$.

So, how do we go about finding the $x_i$ and $w_i$?  Suppose we want to determine the $w_i$ and $x_i$ for $n=2$:

$$
\begin{align}
J_2[f] &= \sum_{i=1}^2 f(x_i) \int_{-1}^{1} L_i(x) = w_1 f(x_1) + w_2 f(x_2)
\end{align}
$$

Evidently, to solve for the unkown $w_i$ and $x_i$, we need 4 equations.  Let us form them by considering the first four polynomials $f=1, x, x^2, x^3$:

$$
\begin{align}
w_1 \cdot 1 + w_2 \cdot 1 = \int_{-1}^1 dx = 2 &\qquad
w_1 \cdot x_1 + w_2 \cdot x_2 = \int_{-1}^1 x dx = 0 \\
w_1 \cdot x_1^2 + w_2 \cdot x_2^2 = \int_{-1}^1 x^2 dx = \frac{2}{3} &\qquad
w_1 \cdot x_1^3 + w_2 \cdot x_2^3 = \int_{-1}^1 x^3 dx = 0
\end{align}
$$

After some algebra and using the methods of undetermined coefficients to solve the nonlinear system of equations, we find that

$$
w_1 = w_2 = 1, \qquad x_1 = -x_2 = -\frac{\sqrt{3}}{3}
$$

and

$$
J_2[f] = f\left(-\frac{\sqrt{3}}{3}\right) + f\left(\frac{\sqrt{3}}{3}\right)
$$

This is called a 2 point Gauss-Legendre quadrature rule.  **The Gauss-Legendre quadrature rule gives exact results if $f(x)$ is a polynomial of degree $2n-1$.**

How can we find the $w_i$ and $x_i$ for higher order Gauss-Legendre quadrature?  We can go through the same process as with $n=2$, but it is slow and tedious.  Instead, using properties of orthogonal polynomials, the $x_i$ and $w_i$ can be derived for many $w_x$ and domains of integration $[a, b]$.  For Gauss-Legendre quadrature, it turns out that the $x_i$ are the roots of the $n^{\text{th}}$ order Legendre polynomial.  The $w_i$ can then be found from knowledge of the $x_i$ and properties of the Legendre polynomials.  The details are beyond the scope of this class, but tabulated values of the $x_i$ and $w_i$ can be found in Table 4.1 of the text.

# <a id='nat_coords'></a> Natural Coordinates[](#top)

To make use of the results from the previous section, the domain of integration for an integral of interest must be $[-1,1]$.  For integration of element functions, this implies that we must perform a change of coordinates from the element domain $[x_1^{(e)},x_2^{(e)}]$ to $[-1,1]$.  The transformed domain $[-1, 1]$ is referred to as the **natural coordinates**.  The actual element domain will be referred to as the **physical coordinates**.  To distinguish natural coordinates from the physical coordinates, the natural coordinates will be labeled $\xi$.

Mathematically, we seek a transformation such that 

$$
\xi = \begin{cases}
-1 & \text{for } x = x_1^{(e)} \\
1 & \text{for } x = x_2^{(e)} \\
\end{cases}
$$

Then the transformation between $x$ and $\xi$ can be represented by

$$
x = a + b\xi
$$

which is known as a "stretch" transformation.  We determine the constants $a$ and $b$ be requiring that

$$
\begin{align}
x_1^{(e)} &= a + b(-1) \\
x_2^{(e)} &= a + b(1)
\end{align}
$$

Solving for $a$ and $b$ we get

$$
a = b = \frac{1}{2} h^{(e)}
$$

Then, the transformation becomes

$$
x = x_1^{(e)} + \frac{1}{2}h^{(e)} \left( 1 + \xi \right)
$$

where $x_1^{(e)}$ denotes the left end node of the element and $h^{(e)}$ is the element length.  

The value of the natural coordinate always lies in $[-1,1]$ with origin at the element center, as shown:



## <a id='appx_geom'></a> Approximation of Geometry

Recall that the element stiffness and force arrays (for second order boundary value problems) are given by

$$
k_{ij} = \int_{x_A}^{x_B} a\phi'_i\phi'_j + c\phi_i\phi_j dx
$$
$$
f_i = \int_{x_A}^{x_B} q\phi_i + P\delta(x - x_i) dx
$$

These element matrices involve derivatives of the interpolation functions with respect to the physical coordinates $x$.  Thus, we need to find a transformation of the form

$$
x = \chi(\xi)
$$

so that the integrals can be cast in terms of $\xi$ on $[-1, 1]$.  The transformation $\chi$ is assumed to be 1 to 1.  An example of one such transformation is

$$
\chi(\xi) = x_1^{(e)} + \frac{1}{2}h^{(e)}\left(1+\xi\right)
$$

which is the same as our original transformation.  In this transformation $\chi(\xi)$ is a linear function of $\xi$.  It will map a straight line to a straight line.  If $\chi$ were non-linear, it would map a straight line to a curved line.

It is natural to think of approximating the geometry in the same way as we approximated a dependent variable.  In other words, the transformation $x=\chi(\xi)$ can be written as

$$
x = \chi(\xi) = \sum_{i=1}^m x_i^{(e)} \phi_i^{(e)}(\xi)
$$

where $x_i^{(e)}$ is the physical coordinate of the $i^{\rm th}$ node of element $(e)$.  $\phi_i^{(e)}(\xi)$ are interpolating polynomials of degree $m-1$. This equation represents the shape of the element and the $\phi_i^{(e)}(\xi)$ are called the **shape functions**, aften denoted as $\{N(\xi)\}$.

Physically, this mapping maps a geometric shape from $\xi$ to $x$.  In other words, for any given value of $\xi$ we can determine the corresponding value of $x$.  When the element is a straight line, the preceding transformation is the same as that for the linear case.

<div class='msg msg-info'>

The transformation is not used to change the actual geometry of the element - the transformation is used to write

$$
\int_{x_1^{(e)}}^{x_2^{(e)}}f(x)dx = \int_{-1}^{1}f(\xi)d\xi
$$

so that Gauss quadrature can be used to approximate the integral.

</div>

To complete the transformation, we need $dx$ in terms of $d\xi$.  $dx$ in the global coordinates is related to $d\xi$ by

$$
dx = \frac{\partial x}{\partial \xi} d\xi = J^{(e)} d\xi
$$

where $J^{(e)}$ is called the "Jacobian" of the transformation.  Using our transformation, we have

$$
J^{(e)} = \frac{\partial x}{\partial \xi} =
\frac{\partial }{\partial \xi} \sum_{i=1}^{m} x_i^{(e)}\frac{\partial \phi_i^{(e)}(\xi)}{\partial \xi}
$$

## <a id='parametric_forms'></a> \*Parametric Formulations

Recall that our approximate solution is of the form

$$
u_N(x) = \sum_{j=1}^{n} u_j \phi_j(x)
$$

In general, the degree of the transformation $m$ does not have to equal the degree of the approximation $n$.  In other words, $\phi_i(x) \neq \phi_i(\xi)$.

Two elements can be used in the finite element analysis: one for the approximation of the geometry and the other for the interpolation of $u$.  Finite element formulations are classified as

- Subparametric: m < n
- Isoparametric: m = n
- Superparametric: m > n

Examples:

- Subparametric: In the Euler-Bernouli beam element the transverse deflection is approximated by the cubic Hermite polynomials and straight beams are approximated by linear polynomials.
- Isoparametric: The same element is used to approximate the geometry and dependent variable.
- Superparametric: The geometry is approximated by higher order polynomial than the dependent variable.  Rarely used in practice.

# <a id='examples'></a> Examples[](#top)

## <a id='ex-1'></a> Example 1: Element Stiffness

We now piece together what we have done so far.  Consider the FE stiffness matrix in global coordinates

$$
k_{ij}^{(e)} = \int_{x_A}^{x_B} a_e \frac{d\phi_i^{(e)}}{dx} \frac{d\phi_j^{(e)}}{dx} dx
$$

To transform the the integral from $x$ to $\xi$, we note that

$$
a(x) = a(\chi(\xi)), \quad dx = Jd\xi
$$

$$
\frac{d\phi}{dx} = \frac{d\phi}{d\xi}\frac{d\xi}{dx} = \frac{1}{J}\frac{d\phi}{d\xi}
$$

Therefore, the stiffness can be written as

$$\begin{align}
k_{ij}^{(e)} &= \int_{-1}^{1} a^{(e)}(\chi(\xi)) 
                \frac{1}{J^{(e)}}\frac{d\phi_i^{(e)}}{d\xi}\frac{1}{J^{(e)}}\frac{d\phi_j^{(e)}}{d\xi}
                J d\xi \\
             &= \int_{-1}^{1} \frac{a^{(e)}(\chi(\xi))}{J^{(e)}} \frac{d\phi_i^{(e)}}{d\xi}\frac{d\phi_j^{(e)}}{d\xi}d\xi \\
             &= \sum_{I=1}^n Q_{ij}^{(e)}\left(\xi_I\right)w_I
\end{align}
$$

where

$$
Q_{ij}^{(e)} = \frac{a^{(e)}}{J^{(e)}} \frac{d\phi_i^{(e)}}{d\xi}\frac{d\phi_j^{(e)}}{d\xi}
$$

For linear interpolation elements and an isoparametric mapping, the shape functions are

$$
\{\phi(\xi)\}^{(e)}
= \frac{1}{2}
\begin{Bmatrix}
1 - \xi \\
1 + \xi
\end{Bmatrix}
$$

and

$$
J = \sum_{i=1}^m x_i^{(e)} \frac{d\phi_i^{(e)}}{d\xi} = \frac{1}{2}\left(x_1(-1) + x_2(1)\right) = \frac{1}{2}h^{(e)}
$$

For elementwise constant values of $a$ and $h$

$$\begin{align}
k_{ij}^{(e)} &= Q_{ij}^{(e)}(\xi_1) w_1 + Q_{ij}^{(e)}(\xi_2) w_2 \\
       &= \frac{2a^{(e)}}{h^{(e)}}\left(2\frac{d\phi_i^{(e)}}{d\xi}\frac{d\phi_j^{(e)}}{d\xi}\right)
\end{align}
$$

from which

$$ k_{11}^{(e)} = \frac{2a^{(e)}}{h^{(e)}} \left(2(-.5)(-.5)\right) = \frac{a}{h}^{(e)} $$
$$ k_{12}^{(e)} = \frac{2a^{(e)}}{h^{(e)}} \left(2(-.5)(.5)\right) = -\frac{a}{h}^{(e)} $$
$$ k_{21}^{(e)} = \frac{2a^{(e)}}{h^{(e)}} \left(2(.5)(-.5)\right) = -\frac{a}{h}^{(e)} $$
$$ k_{22}^{(e)} = \frac{2a^{(e)}}{h^{(e)}} \left(2(.5)(.5)\right) = \frac{a}{h}^{(e)} $$

and

$$
[k]^{(e)} = \frac{a^{(e)}}{h^{(e)}} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}
$$

As noted before, the Jacobian matrix will be the same - $J^{(e)} = h^{(e)} / 2$ - when the element is straight, even when the coordinate transformation is not linear.  However, when the element is curved, the Jacobian is a function of $\xi$ for transformations other than linear.

## <a id='ex-2'></a> Example 2: Reduced Integration

It is often beneficial to reduce the number of integration points that are evaluated.  In the case of the element stiffness with constant properties, single point Gauss-Legendre quadrature can be used to find $k_{ij}^{(e)}$ exactly:

$$
\begin{align}
k_{ij} &= Q_{ij}^{(e)}(\xi_1)w_1 \\
&= \frac{2a^{(e)}}{h^{(e)}}\left(\frac{d\phi_i^{(e)}}{d\xi}\frac{d\phi_j^{(e)}}{d\xi}\right) 2 \\
& = \frac{a^{(e)}}{h^{(e)}} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}
\end{align}
$$

where we used the fact that $w_1=2$ for single point integration.

It turns out that in higher dimensions, using reduced integration can have a dramatic influence on the finite element solution - both good and bad!
