Transcript ENGINEERING 212 98A

Chapter 22
Eigenvalues
Appendix A
Eigenvalue Problems
Engineering Problems involving vibrations,
elasticity, oscillating systems, etc.,
Determine the eigenvalues for n homogenous linear
equations in n unknowns
[ A]{ x }  { x }
[ A]  [ I ]{ x}  {0}
Non-homogeneous system
homogeneous system
λ : eigenvalue; x : eigenvector
Mathematical Background
 a11  
 a
 21
 A   I  x   a31

 
 an 1
a12
a13

a22  
a23

a32

an 2
a33   


an 3

a1 n   x 1 
a2 n   x 2 
  
a3 n   x 3  

   
 
ann     x n 
0 
0 
 
0 
0 
 
0 
 For nontrivial solutions ==> det( A  I )  0
 Characteristic polynomial: det( ) = fn( )
( 1 ) N N  C N 1 N 1  C N 2 N 2    C 2 2  C1   C0  0
The root of fn( ) = 0 are the solutions for the eigenvalues
Mass-Spring System
Equilibrium positions

m1

m
 2
d 2 x1
  kx1  k ( x 2  x1 )
2
dt
d 2 x2
  k ( x 2  x1 )  kx2
2
dt
Mass-Spring System
2
let x1  X 1 sint , x 2  X 2 sint ;  
Tp

m1



m
2


 2 k
k
2
d 2 x1




X

X2  0

1

k
(

2
x

x
)

0


2
1
m1
 m 1

dt 2

2
d x2
 k X   2 k   2  X  0

k
(
x

2
x
)

0
1
2
1
m
 2
 m
dt 2
2
 2


 Homogeneous system

 2 k / m1   2

  k / m2

  X 1  0 
 
2 
2 k / m2     X 2  0 

 k / m2

 Find the eigenvalues  from det[ ] = 0
Polynomial Method
 m1 = m2 = 40 kg, k = 200 N/m
 Characteristic equation det[ ] = 0
( 2 k / m1   2 )
 k / m2
 k / m2
( 2 k / m2   )
2

10   2
5
5
10   2
  4  20 2  75  0  ( 2  15)( 2  5 )  0
 Two eigenvalues  = 3.873s1 or 2.236 s 1
 Period Tp = 2/ = 1.62 s or 2.81 s
Principal Modes of Vibration
Tp = 1.62 s
X1 = X2
Tp = 2.81 s
X1 = X2
Buckling of Column
 Axially loaded column
 Buckling modes
 M – bending moment
 E – modulus of elasticity
 I – moment of inertia
d 2 y M
 2 
EI
 dx
 M   Py

d2y
P
2


y


p
y
2
EI
dx
y( 0 )  y( L )  0
Curvature:
d2y
dx 2
Buckling of Axially Loaded Column
Eigenvalue problem
d2 y
2


p
y ; y ( 0 )  y( L )  0
2
dx
y  A sin px  B cos px
Buckling loads
 y( 0 )  B  0
 pL  n , n  1, 2, 

 y( L)  A sin pL  0
2 2
n
 EI
P  p 2 EI 
L2
Fundamental mode: n = 1
Pcritical 
 2 EI
L2
Euler formula
Buckling
Modes
n
L
2 2
n
 EI
2
P  p EI 
L2
p
Polynomial Method
ODE
d2y M
2



p
y ; y(0 )  y( L)  0
2
EI
dx
Finite-difference method
y i  1  2 yi  yi  1
2
2 2

p
y

0


y

(
2

h
p ) yi  y i  1  0
i
i 1
2
hi
 2  h2 p 2

 1

0




0
1
0
2  h2 p 2
1
1
2  h2 p 2
0
1
0
0
  y1  0 
 y   
0
  2  0 
  y3   0 
1
0
   


     
 2  h2 p 2   yn  0 

0
Which
Scheme?
Order of
Errors?
Characteristic equation: (2n)th-order polynomial
det  ( 2  h2 p 2 ) n    0
Polynomial Method
One interior node (h = L/2)
pexact 

L
2 2 2
( 2  h 2 p 2 ) y1  0  p 

( a  10 %)
h
L
Two interior nodes (h = L/3)

 2  h2 p 2

 1
2
pexact  ,
L L
 1   y1  0 
2 2 2


(
2

h
p ) 10




2 2
2  h p   y2  0 
3 3 3
ph  1,  3  p  ,
( a  4.5 %, 17.3 %)
L
L
Polynomial Method
Three interior nodes (h = L/4)

 2  h2 p 2

 1

0

2 3
pexact  ,
,
L L L
  y1  0 
1
0
   
2  h2 p 2
 1   y2   0  
1
2  h 2 p 2   y3  0 
( 2  h2 p 2 ) 3  2( 2  h2 p 2 )  0
p
 ph   2 ,  2  2
4 2 2 4 2 4 2 2
,
,
( a  2.6 %, 10.0 %, 21.6%)
L
L
L
Power Method
Power method for finding eigenvalues
1. Start with an initial guess for x
2. Calculate w = Ax
3. Largest value (magnitude) in w is the
estimate of eigenvalue
4. Get next x by rescaling w (to avoid the
computation of very large matrix An )
5. Continue until converged
Power method also gives you eigenvectors
Power Method
Start with initial guess z = x0
 (k1 )  w k( 1 )
w ( 1 )  Az ( 1 )
z( 2 ) 
w(1)

(1)
k

w ( 2 )  Az ( 2 )
z
(3)

w( 2 )

(2)
k

Az ( 1 )

rescaling
(1)
k
 (k2 )  w k( 2 )
k is the
Az ( 2 )
(k2 )
dominant
eigenvalue

If 1   2   3     n , then
1
k
  2
k
  3
k
    n
k
Power Method
1. Initial guess z (1)  x0  1,1,...,1
T
 normalize z by biggest wk 
Calculate Az (2)  w(2)   z (3) ;  normalize z by biggest wk 
2. Calculate Az (1)  w(1)   z (2) ;
3.
... Calculate Az ( k )  w( k )   z ( k 1)
For large number of iterations,  should
converge to the largest eigenvalue
The normalization make the right hand side
converges to  , rather than n
Example: Power Method
Consider
Start with
Az ( 1 )
 2 8 10
A   8 3 4 
10 4 7 
1
 
(1)
z  x 0  1
1
 
Assume all eigenvalues
are equally important,
since we don’t know
which one is dominant
 2 8 10 1  20
0.9524
   




  8 3 4  1   15  21 0.7143


10 4 7  1  21
1
.
0


eigenvalue eigenvector
Example
Current estimate for largest eigenvalue is 21
Rescale w by eigenvalue to get new x
 20 0.9524
w

 1   
x
    15  0.7143
m ax(abs( w ))  21    

21
1
.
0
  

Check Convergence (Norm < tol?)
 2 8 10 0.9524
0.9524   2.3812
Ax  x   8 3 4  0.7143  21 * 0.7143    1.2382
10 4 7   1.0 
 1.0    1.6188
Norm
Ax  x  ( 2.3812) 2  ( 1.2382) 2  ( 1.6188) 2  3.1343
 Update the estimated eigenvector and repeat
Az ( 2 )
 2 8 10 0.9524 17.619

 

  8 3 4  0.7143  13.762
10 4 7  
1.0   19.381
 New estimate for largest eigenvalue is 19.381
Rescale w by eigenvalue to get new x
17.619 0.9091
w
 

 1 
x

 13.762  0.7101
m ax(abs( w ))  19.381 
  1.0 
19
.
381

 

 2 8 10 0.9091
0.9091   0.1203
Ax  x   8 3 4  0.7101  19.381* 0.7101    0.3594
10 4 7   1.0 
 1.0    0.4496
Norm
Ax  x  ( 0.1203) 2 ( 0.3594) 2  ( 0.4496) 2  0.5880
Example
One more iteration
Az ( 3 )
 2 8 10 0.9091 17.499
0.9243

 





  8 3 4  0.7101  13.403  18.9310.7080
 1.0 
10 4 7   1.0   18.931


 2 8 10 0.9243
0.9243   0.0147
Ax  x   8 3 4  0.7080  18.931* 0.7080    0.1153
10 4 7   1.0 
 1.0    0.1440
Norm
Ax  x  ( 0.0147) 2  ( 0.1153) 2  ( 0.1440) 2  0.1851
Convergence criterion -- Norm (or relative error) < tol
Example: Power Method
Az ( 4 )
Az ( 5 )
Az ( 6 )
Az ( 7 )
 2
  8
10
 2
  8
10
8 10 0.9243  17.513
 0.9181

 



3 4  0.7080  13.519  19.0750.7087

4 7   1.0  19.075
1.0 

8 10  0.9181 17.506
0.9206

 



3 4  0.7087   13.471  19.0160.7084


4 7  
1.0  19.016
1
.
0


 2 8 10 0.9206 17.508
0.9196

 



  8 3 4  0.7084  13.490  19.0400.7085


10 4 7  
1.0  19.040
1
.
0


 2 8 10 0.9196 17.507
0.9200

 



  8 3 4  0.7085   13.482  19.0300.7085
 1.0 
10 4 7  
1.0   19.030


Script file: Power_eig.m
Norm  Ax  x
MATLAB
Example:
Power
Method
» A=[2 8 10; 8 3 4; 10 4 7]
A =
2
8
10
8
3
4
10
4
7
» [z,m] = Power_eig(A,100,0.001);
it
m
z(1)
1.0000
21.0000
2.0000
19.3810
3.0000
18.9312
4.0000
19.0753
5.0000
19.0155
6.0000
19.0396
7.0000
19.0299
8.0000
19.0338
9.0000
19.0322
error =
8.3175e-004
» z
z =
0.9199
0.7085
1.0000
» m
m =
19.0322
» x=eig(A)
x =
-7.7013
0.6686
19.0327
z(2)
z(3)
z(4)
0.9524
0.7143
0.9091
0.7101
0.9243
0.7080
0.9181
0.7087
0.9206
0.7084
0.9196
0.7085
0.9200
0.7085
0.9198
0.7085
0.9199
0.7085
eigenvector
eigenvalue
MATLAB function
z(5)
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
MATLAB’s Methods
e = eig(A)
gives eigenvalues of A
 [V, D] = eig(A)
eigenvectors in V(:,k)
eigenvalues = Dii (diagonal matrix D)
 [V, D] = eig(A, B) (more general eigenvalue
problems) (Ax = Bx)
AV = BVD
Inverse Power Method
Power method give the largest eigenvalue
Inverse Power method gives the smallest
*Eigenvalues of B = A-1 are inverse of eigenvalues
of A (i.e.,  = 1/)
So one could use power method on w = Bx to get
largest eigenvalue of B - smallest of A
Calculating B is wasteful - instead use
1
w  Bx  A x
 x  Aw
Inverse Power Method
Basic power method gives the dominant eigenvalue
Inverse power method gives the smallest eigenvalue
Ax  x ; B  A 1
x  A  1 Ax  A  1 x  A  1 x  Bx
 1
Bx  A x    x ;

1
 1/ 
μ is thedominanteigenvaluefor B  A 1
dominant  smallest 
Script file for Inverse Power Method
Use LU_factor and LU_solve
» A=[2 8 10; 8 3 4; 10 4 7]
A =
2
8
10
8
3
4
10
4
7
» max_it=100; tol=0.001;
» [z,m] = InvPower(A,max_it,tol);
L =
1.0000
0
0
4.0000
1.0000
0
5.0000
1.2414
1.0000
U =
2.0000
8.0000
10.0000
0 -29.0000 -36.0000
0
0
1.6897
B =
2
8
10
8
3
4
10
4
7
[B] = [L][U]
A =
2
8
10
8
3
4
10
4
7
1.0000 12.7826 0.3000 1.0000 -0.5333
it =
1
2.0000 0.7123 0.1205 1.0000 -0.8013
it =
» z
z =
[L]
[U]
it =
it =
3.0000
0.6687
0.1167
1.0000 -0.8152
4.0000
0.6686
0.1163
1.0000 -0.8155
4
»
m
»
x
0.1163
1.0000
-0.8155
m
=
0.6686
x=eig(A)
=
-7.7013
0.6686
19.0327
eigenvector
eigenvalue
MATLAB
function
Smallest
eigenvalue
CVEN 302-501
Homework No. 15
Finish the HW but do not hand in. I will
post the solution on the net.