Reading assignment:
MANE 4240 & CIVL 4240
Introduction to Finite Elements Lecture notes, Logan 10.4
Summary:
Prof. Suvranu De • Newton-Cotes Integration Schemes
•Gaussian quadrature
Numerical Integration in 1D
Only for a linear finite element
Axially loaded elastic bar
y A(x) = cross section at x
b(x) = body force distribution
∫x1
x2
B EB Adx =
T 1
( x2 − x1 )
2
⎡ 1 −1⎤ x2
⎢ ⎥ ∫ AEdx =
⎣−1 1 ⎦ x1
( ∫ AEdx) ( x −1x ) ⎡⎢⎣−11 −11⎤⎥⎦
x2
x1
2 1
2
F (force per unit length)
x E(x) = Young’s modulus
x Element nodal load vector
x=0 x=L x2
x2 f b i = ∫ N i b dx
x2
x1 f b = ∫ N b dx
T
For each element 1 2 x1 x1
Element stiffness matrix
x2 d2x
x2
k = ∫ B EB Adx
T kij = ∫ Bi EBj Adx d1x
x1 x1
Question: How do we compute these integrals using a computer?
dN (x)
where Bi = i
dx
1
Any integral from x1 to x2 can be transformed to the following Trapezoidal rule: Approximate the function f(ξ) by a straight line
integral on (-1, 1) g(ξ) that passes through the end points and integrate the straight
line
1 g(ξ) f(1)
I = ∫ f (ξ ) dξ
−1
f(−1)
Use the following change of variables f(ξ)
1− ξ 1+ ξ
x= x1 + x2 ξ
2 2 -1 1
1−ξ 1+ ξ
Goal: Obtain a good approximate value of this integral g(ξ) = f (−1) + f (1)
1. Newton-Cotes Schemes (trapezoidal rule, Simpson’s rule, etc) 2 2
1 1
2. Gauss Integration Schemes I = ∫ f (ξ ) dξ ≈ ∫ g(ξ ) dξ = f (1) + f (−1)
−1 −1
NOTE: Integration schemes in 1D are referred to as “quadrature
rules”
Simpson’s rule: Approximate the function f(ξ) by a parabola g(ξ)
•Requires the function f(x) to be evaluated at 2 points that passes through the end points and through f(0) and integrate
the parabola f(1)
(-1, 1) g(ξ)
f(ξ)
• Constants and linear functions are exactly integrated
• Not good for quadratic and higher order polynomials f(−1)
f(0)
ξ
-1 1
How can I make this better?
ξ(ξ −1) ξ(1+ξ)
g(ξ) = f (−1) + (1−ξ)(1+ξ) f (0) + f (1)
2 2
1 1 1 4 1
I = ∫ f (ξ ) dξ ≈ ∫ g(ξ ) dξ = f (1) + f (0) + f (−1)
−1 −1 3 3 3
2
Notice that both the integration formulas had the general form
•Requires the function f(x) to be evaluated at 3 points M
I = ∫ f (ξ ) dξ ≈ ∑Wi f (ξi )
1
(-1,0, 1) −1
i =1
• Constants, linear functions and parabolas are exactly
integrated Weight Integration point
• Not good for cubic and higher order polynomials Trapezoidal rule: Accurate for polynomial of
W1 = 1 ξ1 = −1
M=2 degree at most 1 (=M-1)
How to generalize this formula? W2 = 1 ξ2 = 1
Simpson’s rule: W1 = 1/ 3 ξ1 = −1 Accurate for polynomial of
M=3 W2 = 4 / 3 ξ2 = 0 degree at most 2 (=M-1)
W2 = 1/ 3 ξ2 = 1
Generalization of these two integration rules: Newton-Cotes With ‘M’ points we may integrate a polynomial of degree ‘M-1’
exactly.
• Divide the interval (-1,1) into M-1 equal intervals using M points
• Pass a polynomial of degree M-1 through these M points (the Is this the best we can do ?
value of this polynomial will be equal to the value of the function
at these M-1 points) With ‘M’ integration points and ‘M’ weights, I should be able to
• Integrate this polynomial to obtain an approximate value of the integrate a polynomial of degree 2M-1 exactly!!
integral Gauss integration rule
f(1)
f(ξ)
f(−1)
g(ξ) See table 10-1 (p 405) of Logan
ξ
-1 1
3
Gauss quadrature
Example: M=1 (Midpoint qudrature)
M
I = ∫ f (ξ ) dξ ≈ ∑Wi f (ξi )
1 1
−1
I = ∫ f (ξ ) dξ ≈ W1 f (ξ1 )
−1
i =1
Weight Integration point How can we choose W1 and ξ1 so that we may integrate a
(2M-1=1) linear polynomial exactly?
How can we choose the integration points and weights to exactly f (ξ) = a0 + a1ξ
integrate a polynomial of degree 2M-1?
1
Remember that now we do not know, a priori, the location of the ∫ −1
f (ξ) dξ = 2a0
integration points. But we want
1
∫
−1
f (ξ) dξ = W1 f (ξ1) = a0W1 + a0Wξ1
1
Hence, we obtain the identity For M=1 I = ∫ f (ξ ) dξ ≈ 2 f (0)
−1
2a0 = a0W1 + a1W1ξ1
For this to hold for arbitrary a0 and a1 we need to satisfy 2 f(ξ)
conditions f(0)
Condition 1 : W 1 = 2 g(ξ)
Condition 2 : W 1ξ 1 = 0
ξ
-1 1
i.e., W 1 = 2; ξ 1 = 0 Midpoint quadrature rule:
• Only one evaluation of f(ξ) is required at the midpoint of the
interval.
• Scheme is accurate for constants and linear polynomials
(compare with Trapezoidal rule)
4
Example: M=2 Hence, we obtain the 4 conditions to determine the 4 unknowns
1
(W1 ,W2 ξ1 and ξ2 )
I = ∫ f (ξ ) dξ ≈ W1 f (ξ1 ) + W2 f (ξ2 )
−1 Condition 1 : W 1 + W 2 = 2
Condition 2 : W 1ξ 1 + W 2 ξ 2 = 0
How can we choose W1 ,W2 ξ1 and ξ2 so that we may integrate a
2
polynomial of degree (2M-1=4-1=3) exactly? Condition 3 : W 1ξ 1 + W 2 ξ 2 =
2 2
3
f (ξ) = a0 + a1ξ + a2ξ + a3ξ
2 3
Condition 4 : W 1ξ 1 + W 2 ξ 2 =0
3 3
2
∫
1
f (ξ) dξ = 2a0 + a2
−1 3 Check that the following is the solution
But we want W1 = W 2 = 1
1 1 1
∫
−1
f (ξ) dξ = W1 f (ξ1) +W2 f (ξ2 ) ξ1 = −
3
; ξ2 =
3
(
= a0 (W1 +W2 ) + a1(W1ξ1 +W2ξ2 ) + a2 W1ξ1 +W2ξ2 + a3 W1ξ1 +W2ξ2
2 2
) ( 3 3
)
1 1 1
For M=2 I = ∫ f (ξ ) dξ ≈ f (− )+ f( )
−1
3 3
Exercise: Derive the 6 conditions required to find the
1 1 integration points and weights for a 3-point Gauss quadrature
f (− ) f( )
3 3 rule
f(ξ)
Newton-Cotes Gauss quadrature
1. ‘M’ integration points are 1. ‘M’ integration points are
* * ξ
-1 necessary to exactly integrate necessary to exactly integrate a
1
a polynomial of degree ‘M-1’ polynomial of degree ‘2M-1’
2. More expensive 2. Less expensive
• Only two evaluations of f(ξ) is required. 3. Exponential convergence,
• Scheme is accurate for polynomials of degree at most 3 error proportional to ⎛ 1 ⎞ 2 M
(compare with Simpson’s rule) ⎜ ⎟
⎝ 2M ⎠
5
Example Gauss quadrature:
1
I = ∫ f(ξ)dξ where f(ξ) = ξ3 + ξ2 I = ∫ f(ξ)dξ
1
−1 −1
Exact integration 1 1
= f(- ) + f( )
3 3
2
I= Integrate and check! 2
3 = Exact answer!
3
Newton-Cotes
To exactly integrate this I need a 4-point Newton-Cotes
formula. Why?
Gauss
To exactly integrate this I need a 2-point Gauss formula.
Why?
In FEM we ALWAYS use Gauss quadrature
Comparison of Gauss quadrature and Newton-Cotes for the
integral Linear Element
1 2
I = ∫ cos(2πx) dx 1
−1
ξ = −1 ξ =1
Stiffness matrix
Newton-Cotes
k =∫
ξ =1
ξ =−1
T 1 ⎡ 1 −1⎤ 1
4 ⎣−1 1 ⎦ ∫−1
B EB Adξ = ⎢ ⎥ AEdξ = ( ∫ AEdξ ) 14 ⎡⎢⎣−11 −11⎤⎥⎦
1
−1
Nodal load vector
1
f b i = ∫ N i b dξ
1
f b = ∫ N b dξ
T
Gauss quadrature −1 −1
Usually a 2-point Gauss integration is used. Note that if A, E and b
are complex functions of x, they will not be accurately integrated
6
Quadratic Element
k = ∫ B EB Adξ = AE∫ B B dξ
1 T 1 T
Nodal shape functions 1 2 3 −1 −1
ξ ξ = −1 ξ =0 ξ =1
⎡ (ξ −1/ 2)
⎢
2
−2ξ (ξ −1/ 2) ( )
ξ 2 −1/ 4 ⎤
⎥
N1 (ξ ) = (ξ − 1) 1
= AE ∫ −2ξ (ξ −1/ 2)
⎢ 4ξ 2
−2ξ (ξ +1/ 2) ⎥ dξ
2
−1 ⎢ ⎥
N 2 (ξ ) = (1 − ξ 2 ) You should be able to derive these!
( )
⎢ ξ 2 −1/ 4
⎣
−2ξ (ξ +1/ 2) (ξ +1/ 2) ⎥
2
⎦
ξ
N 3 (ξ ) = (ξ + 1) Need to exactly integrate quadratic terms.
2 Hence we need a 2-point Gauss quadrature scheme..Why?
Stiffness matrix
1 1
k = ∫ B EB Adξ = AE∫ B B dξ
T T
Assuming E and A are constants
−1 −1
dN ⎡ dN1 dN2 dN3 ⎤ ⎡ 1
( 2ξ +1)⎤⎥
1
B= =⎢ ⎥ = ( 2ξ −1) − 2ξ
dξ ⎣ dξ dξ dξ ⎦ ⎢⎣ 2 2 ⎦
