II[IIIgmnrgnflnmnlqmnII
new“

A}. Univtfli‘)’ g :-

 

 

This is to certify that the
thesis entitled

DESIGN SENSITIVITY IN DYNAMICAL SYSTEMS

presented by
JOSEPH EUGENE WHITESELL, JR.

has been accepted towards fulfillment
of the requirements for

PH.D. degreein MECHANICAL ENGINEERING I

W

f
Major professor

 

 

Date 8-8-80

 

0-7639

 

 

ovanoug FINES:
25¢ per day per in: ,
hng'mue Liamv mains:

Place in book return to rénové

charge from circulation records .I I '

 

cc-

DESIGN SENSITIVITY IN DYNAMICAL SYSTEMS

By
Joseph Eugene Nhitesell, Jr.

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

1980

ABSTRACT
DESIGN SENSITIVITY IN DYNAMICAL SYSTEMS
By
Joseph Eugene Whitesell, Jr.

Modal analysis is a standard tool used in the design of
dynamical systems. By numerically determining certain eigen-
values and eigenvectors associated with a linear structural model,
its vibratory behavior can be investigated. The difficulty with
this technique is that it may be costly in computer time if
several redesigns (and re-evaluations of the eigenvalues and
eigenvectors) are needed to find a suitable design.

The thesis presents a technique to represent selected eigen-
values and eigenvectors by a Taylor series in a design variable.
Since, with the technique, the Taylor coefficients can be computed
efficiently and to arbitrary order, the Taylor series can be used
to accurately estimate the consequences of a design change. The
technique, therefore, has the potential to substantially reduce
the number of calculations necessary to re-evaluate the eigensystem
after a design change.

Each Taylor coefficient is determined by evaluating a recursive
formula which is derived through an application of generalized in-
verse theory. If the eigenvector is found through an inverse itera-
tion, the technique is especially efficient. Most of computational
effort spent to find the eigenvector can then be re-applied to find
the Taylor coefficients. Two representative examples illustrate the

power of the method for practical design problems.

DEDICATION

To my father and mother.

ii

ACKNOWLEDGEMENTS

There are several people who deserve acknowledgement for
helping me during my studies and through the preparation of
this dissertation. The advice and friendship of Dr. James
Bernard, my major professor, has been invaluable. His quick
mathematical insight and his engineering intuition as well as
his high personal values and his good humor have made working
with him stimulating and enjoyable. It was Dr. Albert Andry
who first stimulated my interest in vibration problems and in
the topic of this thesis. His ability as a lecturer and his good
advice and friendship were major influences during the past three years.
Both Drs. Andry and Bernard were cited last year by Michigan State
University for excellence-in-teaching. Dr. Robert Schlueter, who was my
advisor during my Master's degree studies, and a member of my thesis
committee, has helped me numerous times during my studies here and I am
greatly indebted to him. The West Coast Branch of my thesis committee
consists of Herb, Margie, Marta, Erik, Evan and Lorin Rauch who have
been good friends and they have enriched my life in many ways.

I want to thank especially Dr. J.S. Frame for introducing me to
linear algebra through his extraordinary teaching ability. Attending

his course was a rare privilege.

For the difficult job of typing this thesis I thank Jan
Swift of the Department of Mechanical Engineering. I appreciate
the programming help of Mr. Mark Zykin of the Case Center for
Computer-Aided Design. I have much enjoyed working with Mark.
Finally, I would like to thank Dr. John Brighton, Chairman of the
Department of Mechanical Engineering, for his consistant support

and leadership and for his advice.

iv

Chapter I -
l.l

l.2

Chapter 2 -
2.l

2.2

Chapter 3 -
3.1

3.2

3.3

3.4

3.5

Chapter 4 -
4.l

4.2

4.3

Chapter 5 -
5.1

5.2

.3

Chapter 6 -
6.1

6.2

6.3
Bibliography

TABLE OF CONTENTS

Introduction

Some Historical Notes
Thesis Overview

Mathematical Preliminaries

Concepts from Linear Algebra
Concepts from Functional Analysis

Eigensystem Differentiability

Preliminary Differentiability Results

The Reduction Process

Eigenvector Differentiability

Differentiation with Respect to Several
Parameters

Conclusions

Differentiation of Eigensystems

Differentiation of Fully Analytic Matrices

Eigensystem Derivatives with Generalized
Inverse

Numerical Approach

Design Sensitivity in Finite Element Models
of Dynamical Structures

The Finite Element Method
Dynamical Problems
Examples

Discussion and Further Study
Discussion

Methods
Suggestions for Further Study

(JON

l9
25
25
28
32

37
42

43
43
49
59
64
64
7o
78
94
94
95
96

97

LIST OF TABLES

Table Page
5.3-6 The Effect of a Design Change on an Eigenvalue 85

vi

Figure
5.l-ll

5.3-7
5.3—8
5.3-9
5.3-10
5.3-ll
5.3-12

LIST OF FIGURES

The Triangulation of a Domain
Triangular Element

Eigenvalue Dependence on Design
Baseline Design

First Design Change

Second Design Change

Third Design Change

Fourth Design Change

vii

Page
68

80
86
89
90
91
92
93

Chapter 1

Introduction

 

The numerical determination of static and dynamic responses
for linear structural models has become a routine procedure in
engineering design [l,2]. Typically a sequence of design and
redesign is followed until certain features such as weight, natural
frequency, mode shape, and deflection and stress magnitudes attain
suitable values. At each stage the designer changes the design
and performs a numerical analysis to evaluate the consequences
of the changes.

This approach relies heavily on human judgement to modify
the design at each stage and since the appropriate design modi-
fications may not be obvious, several redesigns and reanalyses
are often necessary. Alternative methods to improve this pro-
cedure have been proposed. In one method, the full structural
model is replaced by a Taylor series which approximates its
behavior in a neighborhood of the current design [3]. Combined
with interactive computer graphics hardware, this method allows
a designer to preview the consequences of a proposed design while
avoiding a costly reanalysis. Another approach is the optimal
design method [3], which combines nonlinear programming methods
with structural design techniques. The procedure identifies

a cost functional with certain important design features and

attempts to minimize it by changing the design. When the method
is successful, a sequence of numerically determined designs con-
verge to a feasible design for which the selected cost functional
attains a local minimum.

A major difficulty arises when these methods are applied to
large problems since the determination of the system's state
variables and especially their sensitivity values may require
an excessive amount of computation. However, the potential
benefits of these methods provide a strong incentive to develop

improved numerical methods.

l.l Some Historical Notes

 

When the design involves dynamical measures such as natural
frequencies or mode shapes, efficient methods for determining the
sensitivities of eigenvalues and eigenvectors to design changes
become important. Although methods for computing the first
derivative of an eigenvalue with respect to some system parameter
have been known since the nineteenth century work of Jacobi [4],
eigenvector derivatives have a much shorter history. The first
work in this area was done by Fox and Kapoor [5], who presented
two methods for determining an eigenvector derivative for the
symmetric eigenvalue problem (K - AM)u = 0. Unfortunately both
methods are computationally inefficient. The first method,
although it requires only knowledge of the specified eigenvalue
and eigenvector, involves the multiplication of an (n + l)xn

matrix by its transpose and the subsequent solution of a fully

p0pulated symmetric system. The matrix multiplication is a
lengthy computation and usually leads to a loss of any sparseness
possessed by the factor matrices. The second method avoids these
problems, but it expresses the derivative of an eigenvector in
terms of a complete set of eigenvalues and eigenvectors. This
formulation has only theoretical value when the eigensystem is
large since it is difficult to extract a full set of eigenvec-
tors. In refs. [6-8] the work of Fox and Kapoor was extended to
non-symmetric systems but no improvement in computational ef-
ficiency was made.

These computational difficulties were first discussed by
Nelson [9]. As a remedy, he proposed a technique by which the
rank (n - l) matrix L — AI is modified by zeroing certain of its
entries. The modified matrix, which describes a system of
equations which must be solved to determine the eigenvector deri-
vative, is well conditioned and retains the sparseness of the
original system. This is a significant improvement since the
problems of mathematical physics and engineering often involve
sparse matrices. In this thesis further improvement in computing
eigensystem sensitivities are presented. Two representative
examples illustrate the power of the method for practical design

problems.

l.2 Thesis Overview

 

Notation and preliminary mathematical concepts are presented

in Chapter II. In Chapter III a survey of eigensystem differen-

tiability results are presented. The presentation is based on

the literature of the perturbation theory for linear operators.
After two example problems the chapter closes with a discussion
of eigensystem dependence on several parameters.

Chapter IV is concerned with methods for differentiating
eigenvalues and eigenvectors. The chapter states some basic
results, followed by an application of generalized inverse theory
which leads to further improvements in computing eigensystem
derivatives. A recursive numerical method is then described for
efficiently computing eigensystems derivatives of arbitrary
order. Each eigenvalue and eigenvector derivative is determined
by solving a sparse triangular system and involves no more than
0(n2) multiplications.

Chapter V begins with a statement of the design sensitivity
problem in finite element models. The methods of Chapter IV are
then extended to the generalized eigenvalue problem (K - AM)u = 0.
The chapter ends with two examples which involve finite element
formulations.

In the first Example (5.3-l), an eigenvalue problem associated
with a single element plane elastic vibration problem is considered.
After introducing a design variable to the problem, a Taylor
series representing an eigenvalue as a function of the design
variable is determined with the technique presented in this
thesis. The results are compared with direct evaluations of the

eigenvalue at various values of the design variable.

In the second example (5.3-7) an eigenvalue which depends on
the boundary shape of a fixed-fixed vibrating plate is estimated
using a first order approximation. The estimated eigenvalue is
compared with a direct evaluation of the eigenvalue for several
design modifications.

Chapter VI presents concluding remarks and suggestions for

future work.

Chapter 2

Mathematical Preliminaries

 

This chapter contains a summary of definitions and theorems
drawn from linear algebra and functional analysis which are used
in the following chapters. The reader may find this material in
many sources such as refs. [10-14].

In this thesis an n dimensional column vector is denoted by
x. Its transpose to a row vector is denoted by xT. The con-
jugation of xT is denoted by x*. To avoid confusion between
the ith vector in an indexed sequence of vectors and the ith
component in a particular vector, 5i denotes an indexed vector
and xi denotes the ith component of the vector x. The jth com-
ponent of the vector 5i is denoted by xij' An nxm matrix is
denoted by A, B, C etc. The entry in the ith row and jth column
of A is denoted by aij' The ith matrix in an indexed sequence of
matrices is denoted by A1. The entry in the jth row and kth
column of the matrix A1 is denoted by aijk'

The real (complex) field is denoted by 32(3). The real
(complex) vector space of dimension n is denoted by 3?"(@n).

The space of real (complex) nxm matrices is denoted by ‘gynxm

(gnxm) .

2.l Concepts from Linear Algebra

 

Definition (2.l-l): A linear space X is a set of elements called

 

vectors which is closed under addition (if XeX and yeX then
x+y = 2eX) and under multiplication by a real or complex

scalar a (if XeX then ax = yeX).

Definition (2.l-2): Let X be a linear space over a field F and let
V be a set of k vectors {54’52""’5k} of X. Then if for
some set of k scalars {c],c2,...,ck} from F, not all zero,

the linear combination

C154 + c252 + ... + Ckék = 0,

 

the vectors {§I’§Z""’§k} are said to be linearly dependent

over F. If instead

c154 + c252 + ... + ckxk = 0

only if each ci=0 then the vectors {54’52"°”5k} are

linearly independent.

 

Definition (2.1-3): A set of vectors {54’52’°"’5k’°"} is said to
span X if every vector XeX can be expressed as a linear

combination of the set {54’52’°"’5k""}'

Definition (2.l-4): A set of vectors {54’52""’§k} of a linear

space V is said to be a Hamel basis if and only if

 

i) the set spans V and

ii) the vectors {54’52""’5k} are linearly independent.

Definitions (2.l-l) through (2.l-4) are concerned with
only the algebraic properties of a linear space (algebraic
linear space). If these notions are combined with the topo-
logical notions of length, distance and convergence, a space
with both algebraic and topological properties results (topo-

logical linear space).

Definition (2.l-5): A norm on a linear space X is a positive real

valued function

 

 

 

 

with the following properties
1) IIXII f 0 if x f O, lIx|| = 0 otherwise
1'1')||9><II=I9I°I|><|L 96F

iii) ||§1+§2||1|l51||+|l52l|

Definition (2.l-6): An inner product is a scalar function,

 

(-,-), of two elements x,yeX such that

X142 + x3) = (51.52) + (51.253)

54: 052) = aIéqgéz), aeF

x,x) > 0 if x f 0, (x,x) = 0 otherwise.

Definition (2.l-7): Two vectors x and y are orthogonal with

 

respect to an inner product (~,-) if

(x,y) = 0.

Definition (2.l-8): A set of vectors 54 forms an orthonormal

 

§§t_with respect to an inner product (.,.) if

Definition (2.l-9): Two sets of vectors 5i and yfi form a bi;

orthonormal set with respect to an inner product (-,-) if

 

(‘é-l :Xj) : 51i-

Definition (2.l-lO): The transpose of a matrix A is the matrix

T- =
A - B such that aij bji'
matrix A is the matrix A* = B such that aij = bji'

The Hermitian transpose of a

 

Definition (2.l-ll): A matrix A is called symmetric (Hermitian)

 

if A = AT (A = A*) or skew-symmetric (skew-Hermitian) if

A = -AT (A = -A*).

Definition (2.l-l2): A matrix A is upper triangular if aij = O

for i > j, lower triangular if aij = D for i < j and diagonal

 

. = 0 for i f j.

if aiJ

Definition (2.l-l3): The matrix A whose elements are all zero with

the exception that a1. = l is given the special symbol, A = e..

J 1J
and is called the ij_matrix unit.

 

Definition (2.l-l4): The product of an mxk matrix A and an kxn

matrix B is the matrix C = AB such that

k
Cij = Z aipbpr
p=l

Definition (2.l-l5): For positive integers, the powers of an nxn

matrix A=A1 are defined by

Definition (2.l-l6): A matrix A is idempotent if A2 = A.

 

Definition (2.l-l7): If the negative integral powers of the nxn

matrix A exist, they are defined by

A.“ = (A-I)n

IO

where A'1 is called the inverse of A and

I

If A" exists the matrix A is called non-singular otherwise

 

it is called singular.

Definition (2.l-l8): The rank r of an mxn matrix A is the maximum

number of linearly independent columns of A.

Definition (2.l-l9): A generalized inverse of an mxn matrix A

of rank r is any nxm matrix AI of rank r such that

 

AAIA = A.

I

The mxm matrix AAI and the nxn matrix A A are each idem-

potent matrices [10].

Theorem (2.l-20): Any solution X of the linear system

AX

H
.<

where A is an mxn matrix of rank r
X is an nxk matrix

Y is an mxk matrix

ll

may be expressed as

X = A Y + A Z

where A satisfies AA = 0 and Z is arbitrary [10].

O 0

Definition (2.1-21): The scalars A and corresponding vectors u

which satisfy the equation

Lu = Au

are called the eigenvalues and corresponding right eigenvectors

 

 

of the matrix L. The eigenvectors of LT are called the left

eigenvectors of L. The subspace consisting of the origin

 

and all the right (left) eigenvectors corresponding to A
is called the right (left) eigenspace of L corresponding

to A.

Theorem (2.1-22): The eigenvalues of a real symmetric matrix

are real.

Theorem (2.1-23): Corresponding to any real symmetric matrix

is a set of orthonormal eigenvectors.

Definition (2.1-24): A matrix A is similar to a matrix B if there

exists a matrix S such that

12

B = 5' AS.

Theorem (2.1-25): If A and B are similar then they have the same

eigenvalues.

Definition (2.1-26): A matrix A is called diagonable if it is
l

 

similar to a diagonal matrix such that 3‘ AS = A, where the

diagonal matrix A has the eigenvalues A1 of A as its entries
and where the ith column of S and the ith row of S'1 are

respective right and left eigenvectors corresponding to A1.

Definition (2.1-27): A matrix A such that AA* = A*A is called a

normal matrix.

Definition (2.1-28): A matrix S is called orthogonal (unitary) if

sTs = I (s*s = I).

 

Theorem (2.1-29): Every real symmetric matrix A may be diagonalized

by an orthogonal matrix S

S AS = A

where the diagonal matrix A has the eigenvalues, A1 of A,

as its entries.

13

Definition (2.1-30): If f(A) is a polynomial such that

m m-l
f(A) fox + f1A + ... + f A + fm

m-1

then

f(A) Am'I + ... + f

m
fOA + f1 m-lA + fmI

is the corresponding matrix polynomial.

 

Theorem (2.1-31): If f(A) is an analytic function defined on a
simply connected domain Dg; i? and if A is diagonable matrix

such that

-1 l

S AS = A and A = SAS-
then
f(A) = Sf(A)S
where f(A) = :E:f(i.)e .
Theorem (2.1-32): If f(A) and A are as defined in Theorem (2.1-31)

then

14

where A, = SeiiS']. The matrix EH is called the ith rank

one constituent idempotent matrix and AiAj = 5iin and

n
2 A]. = I. Furthermore A]. = g]. y_1.' where u]. is the ith
i=1

column of S (right eigenvector) and v: is the ith row of S"1

(left eigenvector).
Alternatively, if A has 5 distinct eigenvalues A, with

multiplicity mi then

S
f(A) = 2 mpg,

i=1
where
mi mi
51: Z Sens-1:2 Aj’ jg{j|ii = A3,}.
3 J

The matrix 5i is called the ith rank mi constituent idempotent

 

5
matrix and AiAj = 61.in and Z A. = I.
i=1

Theorem (2.1-33): (Interpolation formula) The ith rank mi con-

situent idempotent matrix Ai corresponding to a diagonable

matrix A with 5 distinct eigenvalues is given by

(A-A.I)

S
#‘I

LAC—J.

15

Definition (2.1-34): The matrix RZ = R(z) = (A-zI)'] is called the

resolvent matrix of A.

Theorem (2.1-35): (Residue Theorem) Let F be a closed curve in
the complex plane enclosing eigenvalues A1 of a matrix A,

i = 1,k. Then

 

k
_ 1
Z A,- - - 21”.] R(z)dz.
i=1 r
Theorem (2.1-36): (Dunford-Taylor Integral) Let f(z) be an analytic
function defined on a simply connected domain D g g and

let P be a closed curve in ‘6’ enclosing all of the eigen-

values Ai of a matrix A. Then

 

if R is defined for 21 and 22.

Definition (2.1-38): The matrix AA = A - AI is called the charac-

teristic matrix of A. The equation det(AA) = 0 is called

 

the characteristic eqpation of A. The notation det(-)

 

denotes the determinant of a matrix [10].

16

Theorem (2.1-39): (Cayley-Hamilton) Every matrix satisfies its

characteristic equation.

Definition (2.1-40): A function u(x): €+€n is called a

vector-valued function.

Definition (2.1-41): A function A(x): 952+ gm" is called a

matrix-valued function.

Definition (2.1-42): Suppose D is a simply connected domain in if ,
i‘I‘

L(x) is an nxn matrix-valued function of x e D and {Ai(x) 1

is a set of eigenvalues of L(x). Then if Ai(x) + A(x])
as x + x1 where A(x]) is an eigenvalue of L(x1) with multi-

plicity m and x1 2 D, the set {Ai(x)}? is called the

A-gY‘OUE.

Now let P be a closed curve in i? enclosing a A-group of

eigenvalues. Then

TIT

P(x) = --217 J{.R(z]x)dz where R(z]x) = (L(x) - 21)"1
P

is called the total projector1 for the A-group [14].

 

 

1Note that P(xo) is identical to the rank m constituent idempotent
matrix corresponding to A(x0). Note also the connection to Theorem

(2.1-35) which shows that the total projector is the sum of the in-
dividual consituent idempotent matrices of the A-group. That P(x) is
also idempotent follows from Theorems (2.1-32) and (2.1-35).

17

Theorem (2.1-43): (Product Rule)
i) If A(x) and B(x) are differentiable matrix-valued functions

such that C(x) = A(x)B(x) then

0.0.
xlci
ll
3:
0.0..
XW
+
9—0.
><3>
CD

ii) If u(x) is a differentiable vector-valued function such

that v(x) = A(x)u(x) then

where gig-ord a—-indicates the differentiation of the individual

entries Cij or Vi with respect to x.

Theorem (2.1-44):'If A(x) is differentiable and n is a positive

integer

n

WEE: -1 dA An- i

i=1

and if A(x) is also non-singular

dA'" _ -n 6A"
3;- - - A ai‘ A

Definition (2.1-45): An eigenvalue A of a matrix A is called simple
if it is not repeated. An eigenvalue A of multiplicity m is

called semi-simple if its eigenspace has dimension m. All the

18

-T(I~
' TF'I. ..

eigenvalues of a diagonable matrix are semi-simple.

2.2 Concepts from Functional Analysis

 

Definition (2.2-l): Let {xn} be an infinite sequence of elements
in a normed space X such that lim||xn - x|| = 0. Then {xn}

n-mo

is said to converge (strongly) to x.

Definition (2.2-2): A sequence {xn} in a normed space is called

a Cauchy sequence if for any 6 > 0 there is an N(e) such that
||xm - xnll < e for every m, n > N.
Definition (2.2-3): A normed space X is said to be complete if it

contains the limit point of every Cauchy sequence in X. A

complete normed linear space is called a Banach space.

 

Definition (2.2-4): An inner product space is a linear space X

 

on which (x,y) is defined for each pair of elements x,y in
X. A complete inner product space is called a Hilbert
space. (Since every inner product generates a norm, any

Hilbert space is also a Banach Space.)

19

Definition (2.2-5): A linear operator L is a mapping between

 

linear spaces such that
i) The domain 99(L) of L is a linear space and the range
159(L) lies in a linear space over the same field F

ii) For all x,y c @(L) and aeF

L(x + y) = Lx + Ly
L(aX) = aLx.

Definition (2.2-6): The Null Space of L denoted./V(L) is the set

 

of all x c @(L) such that Lx = O.

The following standard mapping notation will be used. A

function f which maps from a set A into a set B, is denoted by f:A + B

l. A function f:A + B is surjective (onto) if and only if each

 

b e B is an image of some element of A.

2. A function f:A + B is injective (one to one) if for each b
e 39(f) there is exactly one a e A such that b = f(a).

3. A function f:A + B is bijective (one to one and onto) if and
only if it is subjective ggg_injective (that is if and only

if every b e B is the unique image of some a e A.)

Definition (2.2-7): A function f:A + B is called a homeOmorphism

 

if it is bijective and if both f and f-]:B + A are continuous.

20

I‘OI'J l i‘ |-: : '..‘

Definition (2.2-8): Let X and Y be normed spaces and L:£3(L) + Y
a linear operator, where @(L)§ X. The operator is said to
be bounded if there is a positive real number c such that for all

Xe@(L),
IILXII :C ||X||o

Theorem (2.2-9): If a normed space X is finite dimensional, then

every linear operator on X is bounded.

Theorem (2.2-10): If L is a linear operator such that the mapping
L: £Z(L) + Y is injective then there exists the mapping

-1

L :9I’(L) + @(L) which maps each y e 92(L) onto x 59(L)

for which Lx = y.

Theorem (2.2-ll): Let X and Y be linear spaces and L:QD(L) + Y.
Then
a) L']:.§?(L) —> @(L) exists if and only .A/(L) is empty.
b) L"1 is a linear operator.

Definition (2.2-12): Let X be a linear space, a function f defined

for each u e X is a linear functional on X if

 

f[au + 8v] = af[u] + Bf[v]

for all u,v of X and a,8 e F.

21

Definition (2.2-13): The linear space consisting of all the bounded

linear functionals defined on X is called the dual space of

 

X and is denoted by X'.

Theorem (2.2-l4): Let X and Y be Hilbert spaces with inner
products (-,-)x and (-,-)y respectively and L:X + Y a linear
operator. Then corresponding to L there exists a unique

operator L*:Y + X called the adjoint of L which satisfies
* =
(L v,U)X (v,LU)y

for every u e X and v e Y. If L* = L then L is self adjoint.

Definition (2.2-15): A linear self-adjoint operator L is called

positive definite (u,Lu) > 0 unless u 0.
Theorem (2.2-l6): Let A be a non—zero scalar. Let X be a Hilbert
space with inner product (-,-)x and L:X + X a bounded linear

operator. Then if

Lx - Ax = y (I)
Lx - Ax = 0 (2)
L*f - Af = g (3)
L*f - Af = 0 (4)

22

i) (1) has a solution x for every y e X if and only if (2)
has only the trivial solution x = 0.
ii) (3) has a solution f for every 9 e X if and only if (4)
has only the trivial solution f = 0.
iii) (1) has a solution x (is normally solvable) if and only
if (f,y)x = 0 for all solutions of (4).
iv) (3) has a solution f (is normally solvable) if and only

if (g,x)X = 0 for all solutions x of (2)-

Definition (2.2-l7): Let X and Y be Banach spaces. Then a function
f:X + Y is called Frechet differentiable at x e X if there

exists a unique linear operator f'(x):X + Y such that

Hm llflx + h) - f(x) - f'(x)hll . 0.
llh||+O Ilhll

Definition (2.2-18): Let X be a Hilbert space with inner product
(-,-)X and then if f:X-+3? is a Frechet differentiable
real valued function there exists a unique vector in X

denoted by vf(x) called the gradient of f at x such that

f'(x)h = (h,vf(x))x.

Definition (2.2-l9): Let X be a Hilbert space with inner product

(-,-)X and f:X-+3? , then the directional derivative of f at
x is defined by

23

 

provided that this limit exists. If this limit exists for

any direction h then f is called Gateaux differentiable at

 

XeX.

Remark: If f is Frechet differentiable then Df(x,h) exists for

any direction h and
Df(x,h) = '(x)h = (h,vf(x))x.

Theorem (2.2-20): (Projection Theorem [15]) Let X be a Hilbert
space with inner product (-,-)X, let ||-||X be the norm
generated by (-,-)X and let M be a closed subspace of X.
Then corresponding to any x e X there is a unique vector

m e M such that

||v - m*||X §_||v - m||X for all m e M. (2.2-21)

Furthermore (x - m*,m) = 0 is a necessary and sufficient
condition that m* e M is the unique vector which satisfies

(2.2-21).

24

Chapter 3

Eigensystem Differentiability

 

In this chapter various results related to the differentia-

bility of eigenvalues and eigenvectors are presented.

3.1 Preliminary Differentiability Results

 

Suppose L(x) is an nxn matrix-valued function (See 2.1-4) whose
elements are analytic functions of x e D _C_ 93’. It is natural to
ask if and when L(x) has eigenvalues and eigenvectors which are
also analytic functions of x. This question, which has been
studied extensively by Rellich, Kato and others [15-17], un-
fortunately does not have a succinct answer, since the problem
inherits all of the complications associated with repeated eigen-
values and defective eigenspaces from the underlying eigenproblem.
Nevertheless it is vital to recognize whether the eigenvalue problem
behaves smoothly if eigenvalue and eigenvector derivatives are to
be used in computations.

The following theorem gives an important sufficient con-
dition for eigenvalue differentiability. Its proof may be found

in cited literature.

Theorem (3.1-l): Let L(x) be an nxn matrix-valued function of x e 3? ,

which is analytic in some neighborhood of x = 0. Then

25

a) If A(O) is a non-repeated eigenvalue of L(O),
A(x) is an analytic function in some neighborhood
of x = O and there exist left and right eigenvectors
v(x) and u(x) which are analytic vector-valued func-

tions of x in some neighborhood of x = 0.

b) If A(O) has multiplicity m and {Ai(x)}m

1 is the A-group,

Ai(x) may be expanded in a Puiseaux series [18] as

1 2
Ai(x) = A(O) + ulxp + nzxp + ... , p is an integer :_m

Proof: [14, 18 and 19].

It is easy to show that the conditions of Theorem (3.1-l) (a)
are not necessary by constructing a matrix which has repeated

analytic eigenvalues and analytic eigenvectors.
Example (3.1-2): Let
x) (3.I-3)

where the Ai(x) are analytic scalar functions not necessarily
distinct for x e D g; Q and the _u_1.(x) are analytic vector-
valued functions of x e D such that the matrix S(x), whose

ith column is g4(x), is invertible. Then S'I(x) is also

26

analytic [l4], and we may take its ith row as v4(x).

Theorem (2.1-32) shows that the Ai(x) are analytic eigen-
values of L(x) and the vi(x) and ui(x) are analytic left and
right eigenvectors of L(x), even if Ai(x1) = Aj(x]), i f j
for some isolated value x1 (exceptional point). So L(x) has

the desired properties.

The matrix L(x) in example (3.1-2), although it has repeated
eigenvalues, is diagonable due to the manner in which it was
constructed [10]. That an analytic matrix is diagonable however,
may not be taken as a necessary or sufficient condition for

eigenvalue analyticity as the next example shows.

Example (3.1-4): The matrix

[0 x O
L(x) = 0 0 x
x 0 l

.. J

 

 

is diagonable for x3 f 4/27 yet its eigenvalues cannot be
written as a Taylor series in x at x = D [14] whereas the

matrix

27

L(x) =

is not diagonable for any x 5'8’ but its eigenvalue x is

clearly analytic.

3.2 The Reduction Process

 

To gain further insight into the behavior of eigenvalue

differentiability the reduction process of Kato [14] is useful.

 

Suppose L(x) is an analytic matrix-valued function with the ex-

pansion

L(x) = j[: xl L(l)(0). (3.2-l)

If A(x) is an analytic eigenvalue of L(x) with multiplicity m

then a particular branch of A(x) may be expanded as

A(x) = A(O) + xA(])(O) + x2A(2)(O) + ... + ... . (3.2-2)
Kato has demonstrated that under some circumstances the sequence
of coefficients A(1)(O) may be identified with certain eigen—

values of a sequence of matrices, L(i)(x), which are described

below. The following theorem is central to the process.

28

Theorem (3.2-3): Let L(x) be an analytic matrix with expansion
L(x) = Z x‘LI‘)(0) (3.2-4)
i=0

and let A(O) = A be a semi-simple eigenvalue of L(D) = L
of multiplicity m with total projector P(D) = P. If Agl)
is an eigenvalue of L(]) = PL(I)P in the subspace 99(P)

with corresponding constituent idempotent matrix Pg1) then

L(x) has exactly m§1) = dim Pgl) repeated eigenvalues of the

form A + XA§]) + 0(x).

Proof: From the first resolvent equation Theorem (2.1-37) we find

that
(L(x) - AI) R(z,x) = I + (z - A) R(z,x) . (3.2-5)

Suppose that r is a closed curve in ‘8’ which encloses

the A-group eigenvalues of L(x) then

 

(L(x) - AI) P(x) = - I “/7 (z - A)R(z,x)dz (3.2-6)
P

an

follows from applying theorem (2.1-35) to (3.2-5) and re-
calling Definition (2.l-42). Since R(z,x) is analytic in x

[14] (3.2-6) can be written as

29

(L(x) - AI) P(x) = jg: xl£(l)(0) (3.2—7)

where

 

L(1)(0) = — 1"}f (z - A) R(i)(z,0) dz . (3.2-8)

Since A is a semi-simple eigenvalue both sides of (3.2-7) must

vanish at x = D so L(O)(O) = 0 and

“,8 (L(x) )2 AUPIX) = L(1)(O)
X—>

Now if uj(x) is an eigenvector of Aj(x) where Aj(x) is a member

of the A-group then

 

 

lim (L(x) - All P(x)u.(x) = lim j .
X+0 X J X+0 X J

which shows that A§I)(O) is an eigenvalue of L(])(O). To complete
the proof we must evaluate (3.2-8) for i = 1.

From Theorem (2.1-44)
R(I)(z,0) = - R(z,0) L(])(O) R(z,0)
we may expand R(z,0) as

30

R(z,0) = (z — A)-]P(O) + 5(2)

where S(z) has no singularities so

—1

Amen) = [(2 - A>“P<0) + S(zn Wm) [(2 - A) 9(0) + 5(2)].

Integrating (z - A)R(I)(z,0) along r shows that

This theorem can be re-applied to L(])(x) if its eigenvalues
are also semi-simple to derive a higher order expansion for

A(x). In this case we have
A(x) = A+ Xigll + x2 g2) + 0(x2), k = i,m§‘) (3.2-5)

where the Aji are the repeated eigenvalues of

~(2) _ (l) (2) (l) (1)(1)+(l) (l)
L Pj L Pj - Pj L LAL Pj (3.2-6)

in the subspace 3?(P§])). The matrix L: is defined by

n
-l . . .
L: = - ii: (A - A11) Li if L 15 diagonable. (3.2-7)
i=1

The reduction process then consists of repeated applications of

Theorem (3.1-8) to the matrices L(i).

31

The reduction process may also be used to generate expan-
sions for the constituent idempotent matrices, so if it can be
shown that each L(i) has semi-simple eigenvalues then the process
may be continued indefinitely to show that the eigenvalues of
L(x) are analytic and to establish the existance of analytic
eigenvectors (See Theorem (4.1-4)).

In the next section, the question of eigenvector differentia-
bility is discussed through examples involving the interpolation

formula Theorem (2.1-33) for constituent idempotent matrices.

3.3 Eigenvector Differentiability

 

Eigenvector differentiability is more complicated than eigen-
value differentiability. This is due, in part, to the possibility
that the dimension of the eigenspace associated with a particular
eigenvalue may change abruptly or disappear entirely. In this
section we shall use the interpolation formula of Theorem (2.1-33)
to study the behavior of eigenvectors corresponding to analytic
semi-simple eigenvalues in two examples. In the first example
the matrix L(x) has no exceptional point (See Example (3.2-2)) in

its domain D.

Example (3.3-l): Suppose the eigenvalues of L(x) can be represented
by the s analytic functions A](x), A2(x),...,AS(x) where Ai(x)
is assumed to have multiplicity mi and where L(x) is analytic

and diagonable for x e D.

32

Since L(x) is diagonable, we may express the rank mi con-
stituent idempotent matrices of L by the interpolation
formula of Theorem 2.1-33.

For the ith constituent idempotent matrix, Li’ we have

 

L.(x) = j , x e 9 (3.3-2)
“ LI (L(x)-i (x))

If there is no exceptional point in D, the Lj(x) are
analytic on D. Any column of L4(x) which is non-zero on D
represents an analytic eigenvector gg(x) corresponding to
Ai(x) (See Theorem (2.132)). We now assume, without loss of
generality, that x = 0 is in D. Then if g4(0) is an eigen-
vector of L(D) corresponding to Ai(0), Li(x)gg(0) is an
analytic eigenvector corresponding to Ai(x) in some neigh-

borhood of x = 0. This is true since (1)

S
(A.(0)-A
L.(D) .(0) = ‘
‘1 3‘ U (A1.(O)-A
.fl.

 

.(0
J u.(O) =u.(0)7‘0
3‘0

LIL].

and the continuity of L4(x) imply that g4(x) is non-zero in
a neighborhood of x = 0 and since (2) by the Cayley-Hamilton

theorem, Li(x)g4(0) is in the null space of L(x) - Ai(x)I.

33

Now suppose that x = 0 is an exceptional point in D such
that Ak(0) = Ai(0). Then L4(x) is not continuous at x = 0. It may
have a removable discontinuity at x = 0 however.1 In the fol-
lowing example we study the behavior of L4(x) in a deleted

neighborhood of an exceptional point.

Example (3.3-3): Suppose L(x) is an nxn matrix which is dia-
gonable for x e D with n eigenvalues Ai(x) which are distinct

analytic functions on D. If A 0) = Ai(0) for x = 0 e D

k(
then x = 0 is an exceptional point in D since by hypothesis
Ak(x) and Ai(x) are distinct analytic functions. We will now

investigate the behavior of

 

S
(L )-A.( I
Li(x)y_(0) = n (X 3 X) ) 2(0) (3-3-4)
i=1 (A1(X)-Aj(X))
Jfl

in a deleted neighborhood of x = 0, where g(0) is any
eigenvector taken from the two dimensional eigenspace cor-
responding to Ai(0) = Ak(0).

Both the numerator and denominator of L(x)u(0)

become zero at x = D so

 

lim L1.(x)y_(0) = lim d" dx 3(0) (3.3—5)
X+O X+O

9&3IX) dA (x)

de - dxk

 

1For example if L(x) is a normal matrix on x e D then L§(x) = Li(x) [10]
with ||Li(x)|| = l implying that Li(x) has no singularities on D [14].

34

after applying L'Hopital's rule. If by hypothesis, Ai(x),
Ak(x) have different slopes at x = O, the denominator of (3.3-4)
is non-zero at O and the limit (3.3-5) exists.1 However,
(3.3-5) is an eigenvector only if its numerator is also non-
zero.

To study the behavior of the numerator of (3.3—5) let
P = P(O) be the total projector corresponding to Ai(0)
= Ak(0). Then (3.3-15) may be written

 

lim Li(x)g(0) = lim d ( )- ( ) 3(0) (3.3-6)
+0 +0 . d
x x 3%] x _ Hgk x

since g(0) = Pg(0).

Now suppose that g(0) is an eigenvector of

 

” _ dL(O)

. . . dA.(0) 2
corresponding to its eigenvalue ail (See Theorem (2.l-21)).
Then

d d I d o
P 3%”) Pu(0) = aé-l(0)u(0) = 3%”) Pu(0) (3.3-7)

So the numerator of (3.3-6) is non-zero since otherwise

 

1Note also that the matrix L. (x ) has only rank one in the limit as
x+D since A. (x ) is non- repeated in the deleted neighborhood of x = 0.
An eigenvector of L may be selected from the eigenspace of A. (0)
since the range of P is precisely that subspace.

35

ppm) 139(0) = %%(0) P_u_(0) (3.3-8)
and
P gh—(O) Pg(0) = 3353(0) Pgm) (3.3-9)

would contradict (3.3-7).
If g(0) is an eigenvector of L corresponding to

(0) a similar analysis shows that the numerator of

AQQ
boxy

.1-16) does vanish. So if g(0) contains a component

which is an eigenvector corresponding to %§i(0), (3.3-5) is

an analytic eigenvector passing through the exceptional point
x = 0 which corresponds to Ai(x). An analytic eigenvector cor-

responding to Ak(x) can be constructed similarly.

As the final t0pic of this chapter, we consider the question
of differentiability of eigensystems which are functions of
several parameters. As we shall see, the Frechet differentia-
bility of such eigensystems may not be assumed for repeated

eigenvalues even if they are analytic in a single variable.

36

3.4 Differentiation with Respect to Several Parameters

 

The preceding sections concerned differentiability of eigen-
systems with respect to a single variable. The situation is more
complicated when several variables are involved.

Consider the following example [14]:

Let

L(x],x2) = , (x],x2) e 3.2. (3.4-l)

Even though L(x],x2) is symmetric and Frechet differentiable,

2 2)l/2

its eigenvalues A1 2 = :(x1 + x2 are only Gateaux differentiable

when x1 = x2 = 0 [20].
Using Definition (2.2-l9) the directional derivative of

A1 when x1 = x2 = O is determined as

 

((X 'I' th )2 + (X + th )2)I/2 I (X2 + X2)]/2
I I 2 2 l 2
D(A1 2.11) - Ilm t
’ t+0
- :_(h( + h§)I/2. where n = (h,,n2). (3.4-2)

Since D(A1 2,h) exists for all h, A1 2 are Gateaux differentiable,

but they are not Frechet differentiable since D(A h) is not

1,2’
linear in h.

37

This complication has serious consequences if these deri-
vatives are to be used in non-linear programming routines which
use gradients to determine a descent direction. Since the defini-
tion of the gradient of a function depends crucially on the
linearity of the Frechet differential. Recalling Definition
(2.2-20), the Frechet differential of f may be represented in

terms of the gradient of f as

D(f,h) = (Vf,h). (3.4-3)

where (-,-) is the inner product. It is this equivalence which
provides the justification for using the negative gradient as the
direction of steepest descent since if |Ivf|| = l and ||h|| = l
D(f,h) = (vf,h) attains its minimum value for h* = - vf. When
the Gateaux derivative is non-linear in the direction h, the
descent direction is not so easily determined [21, 27].

We may also use Theorem (3.2-3) to compute the directional
derivative for A1 or A2. Since L(x],x2) is Frechet differentiable

we may represent D(L,h) in terms of its gradient VL as

38

Then

D(L,h) =

 

to which Theorem (3.2-3) may be applied since L(t) is an ana-
lytic symmetric function of t. The directional derivatives of

the eigenvalues A1 2 at x1 = x2 = O are the eigenvalues of

L = E_D (L,h) P = P %%-P = D(L,h)
SO
_ 2 2 1/2
D(A12,h) - :(h1 + h2) as before.

Remark (3.4-4): If D(L,h) is a symmetric matrix and if B is any matrix
representing an orthonormal basis for 99(P) then the directional
derivatives are also the eigenvalues of

T“ T

8 L8 = B PD(L,h)PB = BT

D(L,h)B (3 4-5)

which is a formulation given by Haug and Rousselet [20] which

may be easier to apply than Theorem (3.2-3).

The non-linear dependance of the directional derivatives of

eigenvalues can perhaps best understood by writing L as

39

where although the matrix L is linear in hi its eigenvalues cannot
expected to be, unless .9P(P) has dimension one (if the eigen-
value is simple).

The rank one constituent idempotent matrices and eigenvectors
or L(x1,x2) are even 'less' regular than the eigenvalues. To show
this we construct an expansion of L(x],x2) in terms of its con-

stituent idempotent matrices Li(x],x2) (See Theorem (2.l-33).)

L(x],x2) = A](x],x2) L](x],x2) + A2(x],x2) L2(x],x2) (3.4-6)

 

 

where
L - A2 L - A]
L1 = [ifjriigl , L2 = A2 _ A] . (3.4—7)
or
L (x x ) 1 X1 + (XI+X2)1/2 X2
1 1’ 2 g(x$ + x§)I/2 x2 -X] + (x$+x§)]/2
(3.4—8)
1 Xi ' (XI+X§)]/2 X2
L2(x1.X2) 2(X§ + x91/2 x2 _x] _ (x§+x§)i/2
(3.4-9)

4O

However, the idempotents Li(x],x2) possess no limit as

(x],x2)+(0,0) since

“3...:

limo L](O,x2)
2 1 1

whereas

lim L1(x1,0) = (3.4-ll)
x +0
1 O 0

and similarly

N|-—*

limo L2(0,x2)
2 -l 1

whereas

lim L,(x ,0)
X1+O 2 1

This illustrates that although in this example the rank gpg_

constituent idempotent matrices and therefore the eigenvectors

 

can be continued smoothly through the origin along the x] or
x2 axis, their total limits do not exist at the origin. Con—
sequently, the rank one constituent idempotent matrices and eigen-
vectors at not even continuous at x1 = x2 = 0. They do have

directional derivatives since L(t) is symmetric and analytic

in t (see Theorem (4.1-4)).

41

(3.4-10)

(3.4-12)

(3.4-13)

3.5 Conclusions

 

The theorems of this chapter provide sufficient conditions
needed to justify the differentiation of eigenvalues and eigen-
vectors in a number of situations important to applications. In
the next chapter formulations for the actual eigensystem deri-

vatives are developed.

42

Chapter 4

Differentiation of Eigensystems
As the remarks of the previous chapter illustrate, the
question of eigenvalue and eigenvector differentiability in its
full generality is difficult. Some matrices however, have

properties which lead to relatively simple eigensystem deriva-

tives.

4.1 Differentiation of Fully Analytic Matrices

 

Definition (4.1-1): Let the diagonable matrix-valued function L(x)

have the decomposition

' n
L(x) = ZAi(x)ui(x)vT(x) = z Ai(x)L1.(x) x c 0 gig (4.1-2)

-— -n
i=1 i=1

where the g4(x) and 34(x) belong to an analytic biorthogonal set of

eigenvectors such that 14(x)uj(x)=6ij’ Ai(x) are corresponding

analytic eigenvalues and the Li(x) are rank one constituent

idempotent matrices. Then we say, in this thesis, that L(x) is

fully analytic on D.

 

Remark (4.1-3): Theorem (3.1-l) gives sufficient conditions for
full analyticity of L(x) on D if L(x) has no repeated eigen-

values for any x e D.

43

The following theorems form a basis for a study of eigensystem
derivatives of fully analytic matrices and methods to compute
them. The first of these (Theorem (4.1-4)) states a sufficient
condition for full analyticity which is important in applications.
The theorem was originally stated by Rellich [16] but the proof
here follows Kato [14]. (See also [14] for extensions to normal

matrices.)

Theorem (4.1-4) [14, 16]: If L(x) is a Hermitian analytic matrix-
valued function for x e D C_i .9? then L(x) is fully analytic

on D.

Proof: If L(x) is Hermitian on D then the L(i) in Theorem (3.2-2)
are Hermitian and therefore diagonable. The reduction

process can then be continued indefinitely.D

The next two theorems (Theorems (4.1-5 and 6)) present various
algebraic properties of eigensystems derivatives. Part (e) is
essentially the result given by Jacobi [4] although the derivation

is different.

Theorem (4.1-5): Let L(x) be an nxn fully analytic matrix-valued

function for x e D. Then

ll Tl
dL _ dL. dA.
a> o‘z'ziiaiz”z;fi”i
i=1

44

C) ggflj: 48:5
d) I %%-gj = (AJ Ai)!I ggj = (A1 - Aj) %%4 95’ i f j
9’ VI%ui=%i

where dependence on x is understood.

Proof: To prove part (a) apply the product rule of differentiation

to equation (4.1-2). Part (b) follows after differentiating

L? = Li' For part (c), differentiate v§(x)uj(x) = 6ij' Part

(d) and (e) are proved by forming the product
T dL _ dA T T dL
-"-i 632% ‘ 2314‘ liLkij + 24k ii as?" 13'
and applying (b) and (C).D

Theorem (4.1-6): Let L(x) be an nxn analytic real matrix-valued

function for x e D g_ 8? such that LT(x) = L(x). Then

l’l ll
dL _ dL. dA.
9) 37-2 A1371+Za§iLi
i=1 i=1

45

a?” 93 7 24 dx
d) TdL =(A-A)uT99—=(A A)-d—u-Iu
—i 6% J i—i dx 1 j dx —j

e) 1i 2179—1 =33?

Proof: Similar to (4.1-5) after noting that Theorem (4.1-4)
guarantees that L(x) is fully analytic.o

Remark (4.1-7): In part (c) if i = j we have g; 3&4 = 0, showing

that 3i is orthogonal to ggi. This results from the constant
norm condition imposed by gg(x)gj(x) = 6ij' Note also that
L(x):9P:>3ann.

Further insight into the behavior of eigensystem derivatives
may be obtained by studying the effect of changing an individual
element of a matrix. In the following theorem we will assume that
L(O) is a constant matrix such that L(x) = L(O) + eij x is fully

analytic for x e D.

Theorem (4.1-8): Let L(x) = L(O) + eijx be a fully analytic
matrix-valued function for x e D. Let Ai(x) be an analytic

eigenvalue of L(x), let 34(x) and g (x) belong to an analytic

l

46

set of left and right eigenvectors such that g}(x)pj(x) = 6ij

and let Li(x) = gi(x) 31(x) be the corresponding consti-

tuent idempotent matrix. Then

T dL = dA ___ =
a) it 3734. a‘xi‘ lki Bio Lka‘i
T dL _ T du =
b) is $541; ' (At ' As) is aszt lei y—tj

Proof: To prove part a) apply Theorem (3.2-4) e) with g%(x) = eij
and note that xki Ekj is the jith entry in Lk' For part
b) apply Theorem (3.2-4) d). Part c) follows from part a)
and Theorem (2.1-32). Since the sum of eigenvalue deriva-
tives is a constant function of x, the sum of second deri-

vatives with respect to x is zero proving part d).a

Remark (4.1-9): Since the non-zero columns of the idempotent matrices
Li are eigenvectors corresponding to A1, part (a) above es-
tablishes a strong relationship between an eigenvalue's deri-

vatives and its eigenvectors.

47

The following Theorem (4.1-10) states a formulation, given
by Fox and Kapoor [5] (See also [12]), which expresses the deri-

vative of an eigenvector of a symmetric matrix as a linear com-
bination of eigenvectors.

Theorem (4.1-l0) [5]: Let L(x) be an analytic real matrix-valued

function for x e 05; 3? such that LT(x) = L(x) Then

n
Elfin) = Z] CjIijIX)

(4.1-ll)
J

with

< ) 0 1: J

C. X =

J (xxx) - A360)“ ujm $359460
= gJT-(X) 334m 1 7‘J

. d . . . . .
Proof: Since agfl is a vector in an n dimenSional space, we can

express it as a linear combination of n orthonormal eigen-
vectors of L(x)

Tl
%%i(x) = 2 ck(x)9k(x)
k=1

—j ij j a;— = 0 (See Remark (4.1-6),
we have

48

TdU-_ T : Tg_U_'= : °
Ej ail — cj Ed Ed cj and g4 dxi Ci 0,

applying Theorem (4.1-5) (d) completes the proof.D
Although this formulation has theoretical value it is not

useful when a complete set of eigenvectors is not available [9].

However, we may rewrite (4.1-ll) as

n
du. _ -l T dL
6)? ‘ '2. (ii ‘ ii) -”-i 1.13741
i=1
n
= ' IE: [Iii ‘ *i)-]Lj] gF'Ei = ‘ L: 8%”!i (4'1'12)
i=1
The matrix

T]
L: = - jg; [(Aj - Ai)-]Lj]
j=l

as we shall see, in Section 4.3, can be easily computed and will lead

to an efficient means to compute ggi.

4.2 Eigensystem Derivatives with Generalized Inverse

If L(x) is a diagonable fully analytic nxn matrix-valued

function in some neighborhood of x = 0, then we may write

49

[L(x) — A(x)I]u(x) = (2 x1(L(1) - x(1)1))(2 x1u(1)) = o

by expanding L(x), A(x) and u(x) in a Taylor series about x = 0.

Collecting terms in x1

m i
Z1 x Z (L (111 - A(1'11I)u(~11= o. (4.2—2)
i=0 j=0

Each coefficient of x1 must vanish separately so

i c o 1 o c c

2(LI1'31- A(1'311)u(11= o i = O,l,2,... (4.2-3)
i=0

Ol‘

2LI111uI =2A11‘1)u( i

O,l,2,... (4.2-4)

After subtracting the i = j term from (4.2-3) the eigensystem

derivatives are given recursively by

(L - AI)u = o, i = o (4.2—5)
. 1-] . o o o c
(L - AI)u(1) = - (L(1'11- A(1’J)I)u(31, i = l,2,3,...
i=0
(4.2-6)

50

Since we have assumed that L(x) is fully analytic, these
equations must be normally solvable (See Theorem (2.2-16)) for

each i. Hence if v is any vector such that v c.AV(LT - AI), then
1:] c c o o o
vT :E: (L(l-J) _ A(i-J)I)U(J) z 0 (4.2-7)
J=0

follows from Theorem (2.2-l6). Solving for A111 results in

A(1) _ vT L(1) u

 

_ __________ 4.2-8)
VTU (
and i-l 1_]
VT<2 L(i-j)u(.i) _ A(i-i>u(i))
A(i) = J=0 T 3‘1 , i = 2,3,4,...
V U

where v is any vector such that vTu f 0 and v euV(LT - AI) and

u = u(O) lies on u(x), an analytic eigenvector.

We now turn our attention to solving equations (4.2-6) for

u(1).

Since the matrix LA = L - AI has rank n - m we use a

generalized inverse L; and apply Theorem (2.1-20) to find that

i-l
”(1): _ Li Z (L(i-j) _ A(1-j)I)U(j) + Z, (4.2_‘|0)
j=0

51

where z 530(LA) and Li is any rank n - m matrix satisfying

I _
LALALA - LA'

We shall have immediate use for the following:

Lemma (4.2-ll): Let A, AR and AL be nxn matrices satisfying

AA = 0 and A A = 0. Then if AI

R L
of A,

is any generalized inverse

(4.2-12)
is also a generalized inverse of A.

Proof. Multiplying (4.2-12) on both sides by A results in

#

AA A = A(I - A 1

)A1(I - A )A = AA A = A.D

R L

Now suppose P(O) = P is the total projector associated with

A then

since L is diagonable. From Lemma (4.2-ll) if L: is any generalized

inverse of LA then

>2
A
H
l

'U
v
A

43

P0

. -l3)

52

is also a generalized inverse of Lx. Furthermore, since Pu = u

T T

and v P = v , L: satisfies

Liu = 0 (4.2-14)
and
VTL: = o . (4.2-15)
If
1-] o u u o 0
5(1) z _ L: :2: (L(l-J) A(l-J)I)UIJ), 6(0) _ u
j=0
and
u(1) = 9(1) + c.u, i = l,2,3,... (4.2-16)

then Theorem (2.1-20) guarantees that u(11 is a solution of (4.2-6)
for any constants c1 since ciu euV(LA).

We may use the properties (4.2-14 and 15) to simplify
equations (4.2-8 and 9).

If c1 = 0 for all i, then equations (4.2-8 and 9) become

- T (4.2-17)

53

where u is the value of an analytic eigenvector at x = 0 and v is
any vector satisfying v euW(LT - AI) and vTu # 0.
If the eigenvector derivatives are to be used explicitly,

then the constants Ci should be selected to satisfy some normal-

ization conditions such as

x) u(x) = 1 (4.2-18)
and
(u(x))* u(x) = 1. (4.2-19)
If we assume (4.2-l8), we obtain
v u = v u + civ u (4.2-20)

also

 

T

follows from (4.2.16) since v L: = 0.

Now suppose u is also normalized so that (4.2-l9) holds.

Then

(u + xu(1) + x2u(2) + ...)*(u + xu(1) + x2u12) + ... ) - l

54

must be true and the coefficients of x1, i :_1, must be zero so

2 (u(1'J))*u(J) _ O, 1:]
j=0
or
1-1
(u(0))*u(1) + (u(1))*U(O) = _ (”(1-3))*UIJ)
j=1

Equating real parts of both sides

Re[u*u(1)] = - %-Re[:E: (u(1-1))*u(j)]

u*u - u*fl + ciu*u
and

Ci = u*u(1) _ u“1(1).
Then,

Re[c]] = - Re[u*fi(1)]
and

55

(4.2-21)

i-1
Re[ci] = - %—Re[ :2: (u(1'j))*u(j)] - Re[u*d(1)], i 3_2 (4.2-22)
j=1

must be satisfied for (4.2-18 and 19) to hold.

Remark (4.2-23): If L(x) is a normal matrix then u*u(1) = 0 since

u* = vT. If L(x) is a real symmetric matrix then uT = VT,
c1 = 0
and 111 > ( >
_ l (i-J T J
C, - - 2—[ (u ) u 1. i 3_2

To determine u in the event that A = A(D) has multiplicity m
in (4.2-l) we note that since v in (4.2-13) may be selected as any
vector in the m-dimensional null space of LT(0) - A(0)I and since
1(

L(x) is fully analytic, JV(L 0) - A(0)I) is spanned by the values

at x = O of m independent analytic vectors, which we may take as the

rows of the m x n matrix VT(0) = VT. Then from (4.2-6)
1.-.! c o 0 o o
vT (L(1'11 - A(1'111)u(11 = 0 (4.2-23)
i=0
0?“
i-l . . . i-l
v1 IE: LII-J)u(J) z VT :5: A(i-J)U(J) (4.2_24)
i=0 i=0

56

Similarly we form an n x m matrix U(0) = U where the columns

of U are derived from the analytic right eigenvectors of A and

i-l i-l
v1 2 L(i-J’)U(.i) = VT 2 Um A(i-i) z A(i) (4.2-25)
3:0 1:0
In particular, if i = 1,
VTL(1)U = VTPL(1)PU = v1u A(11 = A(11 (4.2-26)

where A(11 is an m x m diagonal matrix consisting of the m first
derivatives of A and the columns of V and U are normalized so that
VTU = I. Equation (4.2-26) shows that V and U must diagonalize
the matrices L(11 and L111 = PL(1)P. Furthermore, if A(l) has
distinct values on its diagonal then V and U are unique. We now

turn to an example.

Example (4.2-27): Let

The eigenvalues of L(x) are A12 = l :_x and 112) = :_l for all
XlEig 1

Analytic eigenvectors of L(x) may be found from the non-zero
columns of the rank one constituent idempotent matrices L1 and

L2. Using Theorem (2.1-33) we arrive at

57

 

 

 

 

J. I.- l -11
12 2 () 2 2
L (x) = L x =
1 41 1_ 2 _i 1
L2" 2 L 2' '2
Then
[./2 ./2
() 7 (1 7
U x = u x = '
52 b7...

 

 

 

 

are the analytic unit eigenvectors of L(x) and they are the unique

unit vectors that form the columns in the matrix S which diagonalizes

 

 

 

 

 

 

 

 

 

./2 /2 - i ./2 ./2 i
/2 /2 ./2 ./2
L2 ”'21 .1 0. f2- '7] -0 4.

Furthermore if v is any vector such that vTu1 # 0 or vTuZ # 0

and v e .4/(L1 - AI) then
0 l o l
vT u1 = l or vT u2 = -1
l 0 l i
T T

58

Since if vT = (x,y) then

 

 

 

 

 

 

(x,y) 0 l 13-1
./2 ./2
/2 x ——-+ y —-
l 0 7 2 2
d: :1
/2 ./2 ./2 ./2
7+y7 X7+Y7
OY‘
(x,y) 0 l %- / /"
_/2 -x —%—+ y —%—
l 0 ‘2'
= =.."
/2 1% £_ :2
X z'y 2 X 2 Y 2

and we have agreement with (4.2-l7) for i = 1. Note that it is
essential that u (or v) lie on the trajectory of an analytic eigen-

vector for this formulation to work.

4.3 Numerical Approach

In numerical work the algorithm (inverse iteration) [22],
(L - AI)ui+] = V, (4.3-1)

Vi = ui/lluilloo (4-3-2)

59

is often used to determine an eigenvector u corresponding to the

eigenvalue A where A = A + 8. Solving (4.3-1) we have

= (L - AI)'1v.

1 (4.3-3)

ui+i
In practice, the matrix (L - AI) is given an LU factorization and
the iteration (4.3-l) is carried out as a series of back sub-
stitution stages. Since each back substitution stage only in-
volves solving a triangular linear system the LU decomposition
stage is the major effort in the computation [23].

These practical considerations provide a strong incentive
to study the structure of (L - AI)"1 with the hope that it can
be utilized in constructing LI.

If L = L(0) is a diagonable matrix with 5 distinct eigen-

values then

5
~ -1 _ ~ -1
J'=1
s
” -1 _ -1 ” -l
(L - Ail) - - 6 L1 i- )E: (AJ - A1) Lj (4.3-4)
J=1
J'i‘i

follows from Theorem (2.1-32).

Now consider

6O

. -1 . -1
113(1‘ Ll)“- - 411) (1' L1) = 2:31:62 (1' £14.,“ " Li)

S
~ -1
+ (Aj — x,) [j]. (4.3-5)
1

LIL].

';i

The first term of the right side of (4.3-5) vanishes since

(I - L4) Li (I — Li) = 0 is a constant and we have
L1 = (1 — L ) (L — A 1)‘1 (1 — L )
Ai -i i -4
s
-l
= - . - . . 4.3-6
IE: (AJ A,) La ( )
J=1
if]
Furthermore, Li is a generalized inverse of LA = L - AiI since
i i

(L - AiI) (1 - Li)(L — xi(1)'1(1 - Lfi)(L - x11) = L - AiI (4.3-7)

Finally since PiIO)’ the total projector of the Ai - group, is

equal to Li we have (dropping the subscript i),

L: = (1 - P)(L - AI)'1(I - P) (4.3-8)

From Lemma (4.2-11) however,

(I - P)L (I - P) = L: (4.3-9)

I
A

61

is also a generalized inverse of LA but since I - P is idempotent

L: = Li. Therefore we write
i _ -1
LA - (I - P)(L - AI) (I - P) (4.3-10)

and note that L: satisfies conditions (4.2-14) and (4.2-15), namely

This form of L: is especially convenient since (L - AI)’1 may
be computed previously to determine the eigenvector u and since
the idempotent matrix I - P is easily multiplied by a vector

particularly if P has only rank one, i.e.,
(I - uvT)b = b - (va)u (4.3-ll)

OY‘

61(1 - uvT) = 61 - (bTu)vT. (4.3-12)

In practice, instead of directly inverting (L - AI)'1, an LU
decomposition of L - AI is formed.
The following has proved to be an efficient procedure if

A is non-repeated.

62

Form an LU decomposition of L - AI.
Use the LU decomposition in an inverse iteration scheme to

find an eigenvector u.
T) dL

(I - uv ax'“ .

Solve for z in LLz

Solve for y in Luy
1)

Z .

y = y - (va)u + c1u

iii-(I

-UV

where Lu and LL are upper and lower triangular matrices computed

in the LU factorization of L - AI.

If A is a repeated eigenvalue with multiplicity m then it

may be necessary to resort to a technique such as the method of

diagonalizing PL(1)P given in Section 4.2 to find an analytic

eigenvector u, before proceeding with,

3.

. dL
Solve for z in LLz (I - P) 32‘”

Solve for y in Luy

II
N

0.9..
X:

= (I - P)y + clu.

63

Chapter 5

Design Sensitivity in Finite Element

 

Models of Dynamical Structures

 

The finite element method is used extensively for solving
boundary value problems involving partial differential operators.
Its success depends largely on its ability to deal effectively
with complicated domains and boundary conditions. Its major
disadvantage is that it often requires large amounts of computer
time to find solutions. It makes sense then to develop efficient
methods to update existing solutions when small design changes
are made in order to avoid a complete re-analysis.

In this chapter, methods for computing a linear dynamical
system's sensitivity to design changes are discussed and two
illustrative examples are presented. The methods described are
directly applicable to solutions generated by the finite element
method [1, 24-26] or similar methods. To begin we shall outline
how the finite element method can be derived for a static homo-

geneous boundary value problem in the plane.

5.1 The Finite Element Method

 

Consider the boundary value problem
Lu = f in 0,

(5.1-l)
u = 0 on an

64

where L: H1(Q) + L2(Q) is the 2nd order partial differential

operator:

= _ 3 9E. - 2.. 9!.
Lu Sin ax] 3y [q 8y] + cu,

1 : 'flz 9—1.1:
H (0) {u c L2(Q). ax uX and By uy c L2(Q)}. p9Q:C

and f e L2(Q), and 9 is a connected, closed, and bounded domain in
922 with boundary 39.
Now suppose that we define the following inner product

on H1(n) to be

(u,v)Hl =j;(uxvx + uyvy + uv) (5.1-2)

for any u, v e H1(n), and consider the quantity

A

(LU.v) = (u,v)L

with v = O on 39.

If (u,v)L is integrated by parts over the domain 0, an

application of Green's identities yields

(u,v)L .1; (puxvx + quyvy + cuv) + J£;(-puxv + quyv).

.1; (puxvx + quVy + cuv). (5.1-4)

65

It can be shown that [25] if p,q and c are strictly positive and

continuous then positive constants a and B exist such that
OLLIUNHUZ 1 [(u.u)L]1/2 : 8[(U.U)]1/2 (5.1-s)

and so (u,u)L = 0 if and only if u = 0.
Since (u,v)L = (v,u)L (symmetric),

and
(oLu1 + BU1,V)L = a(u],v)L + B(u2,v) (linear)

we have established that (u,v)L is an inner product on
H;(n) = {u e H1(n): u = 0 on an}.
Since (u,v)L is an inner product on H;(Q) we may now apply

the Projection Theorem (2.2-20) to find a numerical solution for u.

To find an approximate solution to 5.1.1 let M C H; be
spanned by the Hamel basis 8 = 1¢l’¢2""’¢m1’ then
m
u = :5: 21¢, (5.1-8)

is a vector in M C Hg. We seek G such that
||LG - f|| 5_ ||Lfi - f|| , for all 8 e M

By Theorem (2.2-20), 6 must satisfy

66

(u. 4.), = (i.ej)

J
Tl

(12::1 2,4,. oJlL = (1243-)
ll

12;] (4,. 63.)sz = (f,¢j) for all oj e B (5.1-9)

which we write in matrix form as

Kz

II
H‘

(5.1-10)

where kij = (¢i’¢j)L'
The functions ¢i e H$(Q) may be constructed by triangulating
the domain a [25] as shown in Figure (5.1-ll). We then require

¢i(x,y) to satisfy

where "j = (xj’yj) is a vertex in the triangulation of n. This

choice of the 1i will tend to cause K to have a banded structure

67

since K is populated by inner products (¢i’ ¢j1L which are non-
zero only if "i and nj are sufficiently close in the triangulation.

Since the matrix K was generated by using the Projection
Theorem (2.2-20) with the inner product (-,-)L on the space Hg(9),
we would, in general expect K to change if the domain 9 were

to change. Thus, it is natural to consider the variable finite ele-

ment problem
K(9)z(9) = f(9) (5.1-12)

where 9 is now considered to be a variable domain. Suppose

T(9) = 9', T: 99:2 +592 (5.1-13)
is a mapping operating on the domain 9 such that

9' = 9(6) = 9(0) + on, o e D (5.1-l4)
where A<:992 is a connected, closed and bounded set whose points

are such that every point in A corresponds to some point in
9(0) (2 392 where the addition operation in (5.1-l3) is vector
addition of the points in 9(0) to the corresponding points in A
and where T is a topological homeomorphism for a e 0. Then if

9(0) is given a triangulation, we may construct a matrix K(a)

such that

69

and the variable problem K(a)x(a) = f(a) results. We now define

the directional derivative of K(a) as

 

DK(9(0),A) = lim “9(0) + 0‘4) ' W401) (5.1-15)

a
9+0

(If this limit exists for any 9 as defined above then K(a) is

Gateaux differentiable and DK(9(O),A) is its Gateaux derivative.)

5.2 Dynamical Problems

 

The distributed parameter forced vibration problem (without

damping) can be written as

2
L1 Lg-(xgc) + L2u(x,t) = f(x,t) (5.2-1)
3t
B(u) = g(x,t) on 39x1 (5.2-2)
C(u) = r(x,0) in 9 (5.2-3)

where L1 and L2 are self adjoint positive definite operators on 9,
u is taken from a function space V(9xT) defined on the product
space 9xT derived from the spacial domain 9 (k = 2,3) and T

is the temporal domain [O,t]. In addition, boundary conditions

70

(5.2-2) defined on a9xT, where 39 is the spatial boundary of 9,
may be given along with initial conditions (5.2-3).

When the finite element method is used to obtain an approxi-
mate solution to (5.2-1, 2 and 3), it is usual to arrive at a

lumped parameter model of the form

M'z'(t) + Kz(t) = f(t) (5.2-4)

where M is an n x n matrix

K is an n x n matrix

z(t) is an n-vector consisting of n = km time varying
coefficients of the spline basis functions 11 e B related to the
triangulation of 9
and f(t) is an n vector resulting from projecting f(x,t) onto
the finite subspace spanned by B.

We are mainly concerned with the variable problem

M(e)'z'(a,t) + K(o()z(a,t) = f(a,t) (5.2-5)

and the related generalized eigenvalue problem [28]

[K(a) — A(a) M(a)]u(a) = 0 (5.2-6)

where 6 accounts for a changing spacial domain

0(a) = 9(0) + 9A

71

This problem can be related to the symmetric eigenvalue problem

in standard form through the use of the following transformation [8],

u = M‘1/2v, (5.2-7)

1/2

where M' is the positive definite square root of M'1 and

dependence on a is understood. Then

M'I/ZKM_1/2v = Av (5-2-8)
Since W”2 is positive definite and since K and M are analytic
functions of a e D, W”2 and therefore M'1/2KM'1/2 are analytic

on D. Thus (5.2-6) is equivalent to the variable eigenvalue
problem in standard form of the previous section and the eigen-
values Ai(a) and particular corresponding eigenvectors vi(a) i = l, n,
are analytic for a e D.

It is not convenient to determine the derivatives of Ai(a)
and Vi(a) from (5.2-8) since the transformation M'1/2(a) may be
difficult to compute and the transformed matrix in (5.2-8) is

usually difficult to differentiate. Instead (5.2-6) is used.

After differentiating and rearranging (5.2-6) we arrive at

due- 91:. 9’1-“ -
[K-AMIE; [dd Au a;M]u (5.2 9)

72

Solving for gE-with a generalized inverse Di, we obtain

du _ I dK dM dA

d'oT"DA[doT' *T‘aa‘mu”
where

DA=(K-AM)
and

DI-(K MI .

A - - A ) satisfy

0 DID =9

and z e uV(K - AM). We may assume that analytic eigen-

vectors exist such that

1 .i' 1.1
we have
'T 1/2 1/2 = T _
ui M M ”j ”i M ”j - aij‘

Differentiating (5.2-12) gives

73

(5.2-10)

(5.2-11)

(5.2-12)

(5.2-13)

(5.2-14)

Now suppose u satisfies (5.2-12). Then from Lemma (4.2-ll)

o: = (1 - uuTM)D:(I - MuuT) (5.2-15)
is also a generalized inverse of DA since
DA uuTM = MuuT DA = 0. (5.2-16)
Substituting 0: in (5.2-10) with z = cu gives
= - D:[%§-- Agg]u + cu (5.2-l7)

Multiplying (5.2-l7) by uT

M and applying (5.2-12) shows that
_ l d
C - - §1 U a“; U (5.2-18)
Again we must determine 0: to complete the formulation. The
approach is similar to the one used in the standard eigenvalue

problem.

Suppose S is a matrix such that

ST[K - AM]S = A - AI (5.2-l9)

74

Then we may expand K - AM as [29]

 

n
K - AM = S(A - AI)ST = jg: (A1 - A)Gi (5.2—20)
i=1
where
U.U.T
_ i 1
G1 - T
u.Mu
l
and
GiMGj = GijGi‘ (5.2-22)

T1
(K - AM)'1 = IE: (A, - i)'1e. (5.2-23)

If PG is the sum of the m matrices Gk associated with eigenvalue

A of multiplicity m then

D: = (I - PGM)(K - AM)‘1(I - MPG) (5.2-24)
is a generalized inverse of K - AM such that

DIMu = 0 (5.2-25)

uTMD: = 0. (5.2-26)

75

If u is determined through the inverse iteration algorithm [26]

(K - AM)u = Mv

i+1 i
(5.2-27)

<
ll

uj/IIU-iII0°

then the algorithm given in (4.3) is essentially unchanged.
The eigenvalue derivative may be found by multiplying (5.2-9)
by vT, v e 34(K - AM) and solving for dA/da

dx v1 (§§-- x§§)u T
a...: T , where v Mu f 0. (5-2-28)
“ v Mu

 

The system (5.1-4) can be modified to include dissipative
terms by adding the term C2 to the left side [28]. We then have the

damped system
n(a)2(u,t) + Ci(a,t) + Kz(o,t) = f(a,t) (5.2-29)
and the related generalized eigenvalue problem

a)M(a) + A(a)C(a) + K(a))u(a) = 0 (5.2-30)

where the positive semi-definite matrix C accounts for viscous

damping.

76

The 'damped' eigensystem of equation (5.2-30) is usually
more difficult to solve than the corresponding system (5.2-6)
without damping. However, under special circumstances the eigen-
vectors of the damped system are also eigenvectors of the cor-
responding undamped system. This is true, for example, if C is

some linear combination of K and M, i.e.,
C = 8K + yM, 8,y e32 (5.2-31)

It has been shown that [30]

KM'1c = CM'1K (5.2-32)

is a necessary and sufficient condition for the eigensystem of
(5.2-30) to be uncoupled when transformed to the modal coordinates
of the corresponding eigensystem of (5.2-6). If (5.2-32) holds
for a e D 999 then the eigenvectors of the damped system do
not depend on C(u) and their derivatives may be computed directly
from (5.2-l7).

The eigenvalue derivatives may be determined by differentiating

(5.2-30) and solving for dA/do to obtain

T 2 dM dC dK
(A ETC-1.1+ 1138-1113;)”

vTCu + 2AvTMu

V

 

(5.2-33)

0.10.
9 >2
ll

77

Remark: Some caution must be taken when using this formula since
the eigenvalues are not differentiable for particular values of
the matrix C. For example, suppose u is an eigenvector of (5.2-30)
then

AZUTMU + AuTCu + uTKu = 0 (5.2-34)

is satisfied by

 

-uTCu :_[(uTCu)2 — 4uTMu uTKu]1/2
A = T ’ (5.2-35)
2u Mu
If
(uTCu)2 = 4uTMu uTKu (5.2-36)

then the mode corresponding to u is said to be critically damped [31]
and the eigenvalue A is not analytic. This can be verified by

attempting to differentiate (5.2-35).

5.3 Examples

 

This section is devoted to two design sensitivity examples
involving finite element vibration models.

In the first example we will use the methods presented in
Chapters 4 and 5 to compute a Taylor series for the lowest natural
frequency of a triangular plane elastic element. The design

variation consists of a changing node (vertex) position. The

78

second example involves a fixed-fixed plate assembled from several

plane elastic triangular elements which is undergoing boundary

shape variations.

Example (5.3-l):

tqrf

Figure (5.3-2).

the stiffness matrix is given as a function of x by

K(x) =
4 0 2(x-l)
4 0
x2-2x+2
symmetric

 

-2 —2(x+1) 2
2(x+l) 0 -2(x+1)
-(x-1) -x2 x-l
x -2x+2 x+l -x
x2+2x+2 -(x+l)
x2+2x+2

79

 

_1

Consider the plane elastic element shown in

Following standard procedures and notation [1,2],

(5.3—3)

 

 

 

(XZ’yZ) (x3’y3)
("190) (190)

Figure 5.3-2 Triangular Element (shown in three design positions)

80

where

hld

 

t is the plate thickness,

 

b1 0
_ l
‘11 b1
a] = x3 - x1 2
a2 = x] - x3 x1 - 1
a3 = x2 - x1 -x1 - 1
2A = b1x1 + b2x2 + b3x3 = 2
V = At.
(1) z d_K _
K dx -
0 0 2 0
0 0 2
2x-l -1
2x+2
symmetric

81

and

y2'5’3
y3'yi

=3’1‘3’2

-2x

2x+2

 

 

(5.3-4)

 

 

I0 0 0 0 o 0
o o o o o
(2) _ 1 de _ t 2 o -2 0
K 779
dx 2 0 -2
2 o
symmetric ' 2
K(11=o i>3

For the mass matrix we use

which is a lumped mass formulation [2].

Since the area A is a constant

If the plate thickness t is set equal to 3, M = I and
D (x) = K(x) - A(x)I .

A

The smallest non-zero eigenvalue of K(O) is determined as

82

(5.3-5)

and its corresponding unit eigenvector is

 

 

We may now use (4.2-16 and 17) to determine a Taylor series for A(x).

For L: = D: we use (4.2-13), i.e.

o: = ([1 - uu1]) (K - i1)'1 (1 - uu1)

to find that

 

 

83

The first 10 Taylor series coefficients to single precision accuracy

are

A(11 = 0.0

A(Z) = -.500000
A(31 = 0.0

A(4) = .138890
A15) = 0.0

A(6) = .001543
A(7) = 0.0

A18) = -.052641
A(91 = 0.0
AI10): .061548

With these coefficients, A(x) may be approximated as

n
A(x) = 2 x1AI1).
i=0

In Table (5.3e6a) a Taylor series representing A(x) is evaluated
for x = 0 to x = 1.5 and for n = 2, 4, 6, 8, 10 and 20. In Table
(5.3-6b) a direct numerical solution for A(x) is presented for
the same values of x. Figure (5.3-7) presents these data
graphically.

A comparison between the values of A(x) given by the Taylor
series and the values of A(x) obtained by direct evaluation show

excellent agreement even when n = 2.

84

.o=_m>cmmwm mo
cowpmzpm>m pomcwo

mpm—mw.o
mkppom.o
mmmmmm.o
mvmmoo.F
Kmmmoo.F
mmwmfip.—
mmoomF.F
cmommm._
__mmmm.P
Nm¢mmm.P
wvmmwm._
mmmmm¢.~
mmpom¢.F
mmmom¢.P
mpomm¢.~
ooooom._

Ax?

 

3V

om.P
o¢.~
om.F
om.—
op.P
00.?
om.o
ow.o
o~.o
om.

om.o
o¢.o
om.o
om.o
o—.o
00.0

x

mvmmmp.o¢
momom¢.o—

mvamo.m
omomo¢._
okump.p
mwowmp.~
omomnp.—
mmxmmm.~
mrmowm._
Fmvmmm._
mqmmmm._
Nmmmm¢.P
mmpmm¢._
mmmow¢._
mpommv.F
OOOOOm.—

8h:

monmmm.m
commmo.m
mommmv.—
m¢mm-.—
mkmuvy._
mmmm¢—.F
mexmm—.F
cmommm.—
NmNNmN.F
mmmmmm.~
mmmmwm.P
emmmm¢.F
mmpom¢.—
mmmom¢.F
mpomm¢._
ooooom._

OP":

n_¢mmN.Ou

NmemN.o
«mmmmm.o
memowm.o
vammm.o
pmmnmo.p
vwmeo_.P
poemNN.F
«memmm.p
mm—nmm.—
mmvmmm._
mmmmm¢.—
mmpom¢._
mmmow¢.F
mfiommv.y
ooooom.—

m n :

mim.m mpnmp

Fonmmo.F
ompmoo._
mmpmmo.P
mommno.F
owopop._
Nm¢o¢P.—
ewmomF.F
mommmm.~
wmmmmm._
mmowmm.~
¢ommmm._
mommmv.F
mNFom¢.P
mmmow¢.P
vpomm¢.—
ooooom.P

m h c

mmpwno.F
mmmmmo.~
ommpmo.~
ooowoo.P
mvmwmo.F
mwwmm~.F
mNme~.P
mmmmmm.~
memwwm.—
ooowmm.—
omommm._
mmmmmv.F
mmpmmv.F
mmmow¢._
e—omme.p
ooooom.P

a n :

.mscmp c mo mmwcmm Lopxmh cpwz mmpms_pmm mzfim>cmm_m

momemm.o
ammo—m.o
mmmvmo.o
oooowm.o
ooommw.o
oooooo.~
ooommo._
oooom—.F
ooommm.~
oooomm.~
ooommm.—
oooo~v.~
QQOmm¢._
oooom¢._
ooomm¢._
ooooom.~

N n c

3

om.—
o¢.F

om
om

F
—

op._
oo.P
om.o
om.o
om.o
on o

om

o

o¢.o
om.o
om.o
o_.o
oo.o

.m

85

:mwmwo co mocmvcmamo m:_m>cmmwm mim.m mczmvu

 

 

 

m.F o.— m. m.i O.Fi m.—i
— LI
(1+ ” _ a
w H:
N N:
COWPGJFQ>U
pooewo
0.7... ‘
OUCQVHC .
OoNIlI
OF H :
ON n :
o.mI:i

 

86

The efficiency of the technique may make design sensitivity
calculations of eigenvalues and eigenvectors feasible in many
practical situations. For example, if an inverse iteration method
is used to determine the eigenvector at x = 0 then subsequent
determinations of each of the Taylor coefficients A(11 and u(11
involve little extra effort since the generalized inverse D: is
easily determined from 4.2-13 (see Section 4.3). (The method is
especially efficient if the matrix K(x) is only linear or quad~
ratic in x since in these cases most of the matrices K(11 = L111
in Equations (4.2-16 and 17) are null.) Since the Taylor coef-
ficients can be obtained so economically, there is a strong in-
centive to represent the design with a Taylor series. Then
subsequent evaluations of the eigenvalue and eigenvector as the
design is modified can be carried out with substantially less
computational effort.

The next example involves a finite element model assembled
from thirty two plane elastic elements similar to the element of
example (5.3-l). The example illustrates the results of a design
sensitivity calculation for an eigenvalue depending on a vari-

able boundary shape.

Example (5.3-7): Figures (5.3-8) through (5.3-12) illustrate the
sensitivity on the lowest eigenvalue of a dynamic finite model
to changes of the boundary shape. Figure (5.3-8a) shows the
baseline design. The lowest eigenvalue for the baseline

design is calculated as 14.44. (Figure (5.3-8b) illustrates

87

the mode shape corresponding to this eigenvalue.) In
Figures (5.3-9) through (5.3-12) the shape of the lower edge
of the model is changed. In each case an estimate of the
new eigenvalue is compared with the result of a finite
element run for the new configuration. The estimates are
based on a single first order directional derivative of the
eigenvalue with respect to some direction h (see Definition
(2.2-10)) which is related to the variable boundary. The
direction h is determined by a gradient projection method as
a direction (in the design space) which causes the eigen-
value to increase without increasing the volume (mass) of
the plate. (For a comprehensive discussion of this and
similar techniques see [3].) Each subsequent variation in
the design is determined by moving along the direction h an
increased distance.

A comparison between a direct evaluation of A and the
first order estimate of A shows good agreement for moderate
changes in the design but for extreme design changes the

estimate is substantially higher than the direct evaluation.

88

 

cmwmmo wcwpwmmm mum.m wczmvm

omqq.¢_ n msfim>zmmwm op mcwucoammccou mamzm one: ADV

4223:8491

’>

moo_a coxwa-ooxea to eo_pe_=memeee Aev

 

89

 

 

 

 

 

 

 

 

 

 

mmcmsu cmwmmo pmcwd m-m.m mczmwm

 

 

 

 

 

 

 

 

 

 

mnmm.o~ o=_e>comem oopsanu

oe¢o.o~ o=_e>eomem eoeeeeemm

9O

 

 

mmcmno cmwmmo ucoomm opim.m mesmwd

 

 

 

 

 

 

 

emem.e~ o=Pe>=om2m eooaaeoo

eeee.m~ o=_e>eom_m eooee_omu

91

 

 

omemeo ememoo eeeee __-m.m seamen

 

 

 

 

 

 

vomm.vm
m¢¢m.~m

mspm>cwmwm umpaqsou

mzpm>cmmwm uwumawumm

92

 

 

mucosa :mwmwo sensed Npum.m mczmwm

 

 

 

umpp.ee ospe>eomem 66336268

_mem.em o=_n>comem emceewomm

93

Chapter 6

Discussion and Further Study

6.1 Discussion

 

Modal analysis has become a standard tool used in the design of
dynamical systems. By numerically determining certain eigenvalues
and eigenvectors associated with a linear structural model a designer
can investigate how a particular structural design will perform
dynamically. The eigenvalues are related to the natural frequencies
of the structure and the eigenvectors indicate the associated
vibratory modes.

If the designer is concerned with dynamical performance,
several redesigns and re-evaluations of the eigenvalues and
eigenvectors may be necessary before a suitable design is found.
Since this procedure is costly in computer time, there has been
an incentive to find improved methods.

This thesis presented a technique to represent an eigenvalue
and eigenvector by a Taylor series in a design variable. Since
the Taylor coefficients can be determined to arbitrary order,
excellent estimates may be easily calculated for design changes
which do not exceed the radius of convergence of the series.

Since the coefficients of the Taylor series can be computed with
substantially less effort than recomputing the eigenvalue or

eigenvector, the procedure has the potential to reduce the number

94

of calculations necessary to carry out the design of dynamical
structures. There is also the potential to use the method in the
implementation of optimal design algorithms.

Two representative example problems were presented. In one
of them, twenty Taylor coefficients of a specified eigenvalue
were determined after a single finite element run. These coef-
ficients were used to represent the eigenvalue as a power series
in the design variable. Comparison between separate finite

element runs and the Taylor series showed excellent agreement.

6.2 Methods

 

Each Taylor coefficient is determined by repeatedly evaluating
a simple recursive formula which is derived through an application
of generalized inverse theory. The procedure is simplified by
using a generalized inverse matrix specially selected to annihilate
certain terms in the formulation. The generalized inverse itself
need not be formed explicitly if the eigenvector is determined
through the inverse iteration method, with the result that each
Taylor series coefficient is computed with no more than 0(n2)
multiplications.

The technique was extended to compute derivatives for the
eigenvalues and eigenvectors of the generalized eigenvalue
problem (K a AM)u = 0 and there was some discussion of the damped

2

eigenvalue problem (MA! + AC + K)u = 0.

95

6.3 Suggestions for Further Study

 

There are several questions involving eigensystem derivatives
which deserve further study.

Necessary conditions for eigenvalue differentiability are in
short supply. It may be possible to use Equation (4.2-25) to
construct a sequence of necessary conditions for the full analy-
ticity of diagonable matrices similar to the sufficient condition
given by the reduction process. (Equation (4.2-26) amounts to such
a necessary condition.)

A convenient computational method for determining the radius
of convergence for Taylor series representing eigenvalues and
eigenvectors would be useful in applications. Lower bounds have
been given [14] which unfortunately often severely underestimate
the actual convergence radius.

Most of the results of this thesis can be extended to eigen-
value problems involving linear operators on Hilbert spaces. Ex-

cellent references for continued study are Refs [10-14 and 19].

96

BIBLIOGRAPHY

[l]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

BIBLIOGRAPHY

Zienkiewicz, O.C., The Finite Element Method, McGraw-Hill, New
York, 1978.

Desai, C.S. and Abel, J.F., Introduction to the Finite Element
Method, Van Nostrand Reinhold Co., New York, 1972.

Haug, E.J., and Arora, J.S., Applied Optimal Design, John Wiley
& Sons, New York, 1976.

Jacobi, C.G.J., Crelle's Journal Far Die Reine Und Angewand
Mathematik, Vol. 30, Berlin: De Gruyter, 1846, pp. 51-95.

Fox, R.L. and Kapoor, M.P., ”Rates of Change of Eigenvalues
and Eigenvectors", AIAA Journal, Vol. 6, Dec. 1968, pp. 2426-2429.

Plaut, R.H. and Huseyin, K., "Derivatives of Eigenvalues and
Eigenvectors in Non-Self Adjoint Systems", AIAA Journal, Vol. II,
Feb. 1973, pp. 250-251.

Rogers, L.C., "Derivatives of Eigenvalues and Eigenvectors",
AIAA Journal, Vol. 8, May 1970, pp. 943-944.

Rudisill, C.S., "Derivatives of Eigenvalues and Eigenvectors

for a General Matrix", AIAA Journal, Vol. 12, May 1974, pp. 721-722.

Nelson, R.B., "Simplified Calculation of Eigenvector Derivatives",
AIAA Journal, Vol. 14, Sept. 1976, pp. 1201-1205.

Frame, J.S., "Matrix Functions and Applications, Part I", IEEE
Spectrum, March, 1964.

Frame, J.S., "Matrix Functions and Applications, Part II", IEEE
Spectrum, April, 1964.

Wilkinson, J.H., The Algebraic Eigenvalue Problem, Clarendon
Press, London, 1965.

Dunford, N. and Schwartz, Linear Operators, Part I, II and III,
John Wiley & Sons, New York, 1971.

Kato, T., Perturbation Theory for Linear Operators, Springer-
Verlag, New York, 1976, 2nd Ed.

97

EEK—...— L___‘S:T" Kl

[15] Luenberger, D.G., Optimization by Vector Space Methods, John
Wiley & Sons, Inc., New York, 1969.

[16] Rellich, F., Perturbation Theory of Eigenvalue Problems,
Lecture Notes, New York University, 1953.

[17] Baumgarter, H., Englichdimensionale analytische Storungtheorie,
Academie-Verlag, 1973.

[18] Knopp, K., Theory of Functions, Part II, Dover, New York,
1947 and 1975.

[19] Ben-Israel, A. and Greville, T.N.E., Generalized Inverses:
Theory and Applications, John Wiley & Sons, New York, 1974.

[20] Haug, E.J. and Rousselet, B., "Design Sensitivity Analysis in
Structural Mechanics: Eigenvalue Variations", Technical Report
No. 52, NSF Project No. ENG 77-19967, Division of Material
Engineering, College of Engineering, The University of Iowa,
Iowa 52242.

[21] Choi, Kyung, Personal Conversation, May, 1980.

[22] Wilkinson, J.H. and Reinsch, C., Handbook for Automatic Compu-
tation, Vol. II, Linear Algebra, Springer-Verlag, New York,
1971.

[23] Forsythe, G.F. and Moler, C.B., Computer Solutions of Linear
Algebraic Systems, Prentice-Hall, Inc.,

[24] Strang, G. and Fix, G.J., An Analysis of the Finite Element
Method, Prentice Hall, Inc., Englewood Cliffs, N.J., 1973.

[25] Prenter, P.M., Splines and Variational Methods, John Wiley
& Sons, New York, 1975.

[26] Bathe, K.J., Wilson, E.L., Numerical Methods in Finite Element
Analysis, Prentice-Hall, Englewood Cliffs, N.J., 1976.

[27] Haug, E.J. and Rousselet, B., Personal Conversation, May, 1980.

[28] Meirovitch, L., Analytic Methods in Vibrations, MacMillan,
New York, 1967.

[29] Lancaster, P. Lambda-Matrices and Vibrating Systems, Pergamon
Press, New York, 1965.

[30] Caughey, T.K. and O'Kelley, M.E.S., "Classical Normal Modes in
Damped Linear Dynamic Systems", Journal of Applied Mechanics,
Vol. 32, 1965, pp. 583-588.

[31] Inman, D.J. and Andry, A.N., "Some Results on the Nature of

Eigenvalues of Discrete Damped Linear Systems", Journal of
Applied Mechanics, ASME, to appear.

98