Transcript MANE 4240 & CIVL 4240 Introduction to Finite Elements Prof. Suvranu De Numerical Integration in 1D.

MANE 4240 & CIVL 4240
Introduction to Finite Elements
Prof. Suvranu De
Numerical Integration in 1D
Reading assignment:
Lecture notes, Logan 10.4
Summary:
• Newton-Cotes Integration Schemes
•Gaussian quadrature
Axially loaded elastic bar
y
A(x) = cross section at x
b(x) = body force distribution
(force per unit length)
x E(x) = Young’s modulus
F
x
x=L
x=0
x1
1
For each element
x2
2
Element stiffness matrix
x2
k   B EB Adx
T
x1
dNi ( x)
where Bi 
dx
x2
kij   Bi EB j Adx d1x
x1
d2x
Only for a linear finite element

x2
x1
B EB Adx 
1
T
 x 2  x1 
2
 1 1 x2
 1 1  x1 AEdx 



x2
x1
AEdx
 x  x 
1
2
1
2
 1 1
 1 1 


Element nodal load vector
x2
f b   N b dx
x1
T
x2
f b i   N i b dx
x1
Question: How do we compute these integrals using a computer?
Any integral from x1 to x2 can be transformed to the following
integral on (-1, 1)
1
I   f ( ) d
1
Use the following change of variables
1
1 
x
x1 
x2
2
2
Goal: Obtain a good approximate value of this integral
1. Newton-Cotes Schemes (trapezoidal rule, Simpson’s rule, etc)
2. Gauss Integration Schemes
NOTE: Integration schemes in 1D are referred to as “quadrature
rules”
Trapezoidal rule: Approximate the function f by a straight line
g that passes through the end points and integrate the straight
line
g
f1
f1
f

-1
g ( ) 
1
1 
f (1) 
f (1)
2
2
1
1
1
1
1
I   f ( ) d   g ( ) d  f (1)  f (1)
•Requires the function f(x) to be evaluated at 2 points
(-1, 1)
• Constants and linear functions are exactly integrated
• Not good for quadratic and higher order polynomials
How can I make this better?
Simpson’s rule: Approximate the function f by a parabola g
that passes through the end points and through f(0) and integrate
the parabola
f1
g
f
f1
f(0

-1
g ( ) 
1
 (  1)
 (1   )
f (1)  (1   )(1   ) f (0) 
f (1)
2
2
1
4
1
I   f ( ) d   g ( ) d  f (1)  f (0)  f (1)
1
1
3
3
3
1
1
•Requires the function f(x) to be evaluated at 3 points
(-1,0, 1)
• Constants, linear functions and parabolas are exactly
integrated
• Not good for cubic and higher order polynomials
How to generalize this formula?
Notice that both the integration formulas had the general form
M
I   f ( ) d  Wi f ( i )
1
1
i 1
Weight
Integration point
Trapezoidal rule:
W1  1
M=2
1  1
W2  1  2  1
Simpson’s rule: W1  1 / 3 1  1
W2  4 / 3  2  0
M=3
W2  1 / 3
2  1
Accurate for polynomial of
degree at most 1 (=M-1)
Accurate for polynomial of
degree at most 2 (=M-1)
Generalization of these two integration rules: Newton-Cotes
• 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
value of this polynomial will be equal to the value of the function
at these M-1 points)
• Integrate this polynomial to obtain an approximate value of the
integral
f1
f
f1
g

-1
1
With ‘M’ points we may integrate a polynomial of degree ‘M-1’
exactly.
Is this the best we can do ?
With ‘M’ integration points and ‘M’ weights, I should be able to
integrate a polynomial of degree 2M-1 exactly!!
Gauss integration rule
See table 10-1 (p 405) of Logan
Gauss quadrature
M
I   f ( ) d  Wi f ( i )
1
1
i 1
Weight
Integration point
How can we choose the integration points and weights to exactly
integrate a polynomial of degree 2M-1?
Remember that now we do not know, a priori, the location of the
integration points.
Example: M=1 (Midpoint qudrature)
1
I   f ( ) d  W1 f (1 )
1
How can we choose W1 and 1 so that we may integrate a
(2M-1=1) linear polynomial exactly?
f ( )  a0  a1
1

1
f ( ) d  2a0
But we want
1

1
f ( ) d  W1 f (1 )  a0W1  a0W1
Hence, we obtain the identity
2a0  a0W1  a1W11
For this to hold for arbitrary a0 and a1 we need to satisfy 2
conditions
Condition1 : W1  2
Condition2 : W11  0
i.e.,
W1  2; 1  0
For M=1
1
I   f ( ) d  2 f (0)
1
f
f0
g

-1
1
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)
Example: M=2
1
I   f ( ) d  W1 f (1 )  W2 f ( 2 )
1
How can we choose W1 ,W2 1 and 2 so that we may integrate a
polynomial of degree (2M-1=4-1=3) exactly?
f ( )  a0  a1  a2 2  a3 3
2
1 f ( ) d  2a0  3 a2
1
But we want
1

1
f ( ) d  W1 f (1 )  W2 f ( 2 )

 
 a0 W1  W2   a1 W11  W2 2   a2 W11  W2 2  a3 W11  W2 2
2
2
3
3

Hence, we obtain the 4 conditions to determine the 4 unknowns
(W1 ,W2 1 and 2 )
Condition1 : W1  W2  2
Condition2 : W11  W2 2  0
Condition3 : W11  W2 2
2
Condition4 : W11  W2 2
3
2
3
2

3
0
Check that the following is the solution
W1  W2  1
1  
1
3
; 2 
1
3
1
For M=2
I   f ( ) d  f (
1
1
f (
1
3
)
f(
3
1
3
) f(
1
3
)
)
f
-1
*
*

1
• Only two evaluations of f is required.
• Scheme is accurate for polynomials of degree at most 3
(compare with Simpson’s rule)
Exercise: Derive the 6 conditions required to find the
integration points and weights for a 3-point Gauss quadrature
rule
Newton-Cotes
1. ‘M’ integration points are
necessary to exactly integrate
a polynomial of degree ‘M-1’
2. More expensive
Gauss quadrature
1. ‘M’ integration points are
necessary to exactly integrate a
polynomial of degree ‘2M-1’
2. Less expensive
3. Exponential convergence,
error proportional to  1  2M


2
M


Example
1
I   f( )d where f( )   3   2
1
Exact integration
I
2
3
Integrate and check!
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?
Gauss quadrature:
1
I   f( )d
1
 f(
2
3
1
3
)  f(
1
3
)
Exact answer!
Comparison of Gauss quadrature and Newton-Cotes for the
1
integral
I   cos(2x) dx
1
Newton-Cotes
Gauss quadrature
In FEM we ALWAYS use Gauss quadrature
Linear Element
2
1
  1
 1
Stiffness matrix
1  1 1 1
k   B EB Ad  
AEd 


 1

1
4  1 1 
 1
T


1  1 1
1 AEd 4 1 1 
1
Nodal load vector
1
f b   N b d
1
T
1
f b i   Ni b d
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
Quadratic Element
1
Nodal shape functions
N1 ( ) 

2
N 2 ( )  1   2 
N 3 ( ) 

2
  1
  1
2
3
  0  1
You should be able to derive these!
  1
Stiffness matrix
1
1
k   B EB Ad  AE  B B d
1
T
1
T
Assuming E and A are constants
d N  dN1 dN 2 dN3   1
B

    2  1  2
d  d d d   2
1
 2  1
2

1
1
k   B EB Ad  AE  B B d
T
1
T
1


   1/ 2 2
2   1/ 2 
 2  1/ 4 
1 

 AE   2   1/ 2 
4 2
2   1/ 2   d

1 
2
2
2   1/ 2    1/ 2  
   1/ 4




Need to exactly integrate quadratic terms.
Hence we need a 2-point Gauss quadrature scheme..Why?