6.5 Matrix Methods

Back

In the preceding section, the equations of motion are solved for the system. Often, one is not interested in

the complete solution. Rather, one is interested in the natural frequencies and the normal modes. When

the number of DoF is very high, only the lowest natural frequency needs to be computed. Several

methods to find some, but not all, of the eigenvalues of the system will be presented.

6-14 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

6.5.1 Example: Three-Degree-of-Freedom System

A three-DoF system will be used as an example.

The larger DoF systems work exactly in the same

manner but produce larger calculations.

Let us consider the following spring – mass

system in Figure 6.6.

The equation of motion for this system is

M

2 0 0

0 1 0

0 0 1

2

6664

3

7775

x€1

x€2

x€3

2

6664

3

7775

þ K

3 21 0

21 2 21

0 21 1

2

6664

3

7775

x1

x2

x3

2

6664

3

7775

¼

0

0

0

2

6664

3

7775

Assuming a harmonic solution of the form x ¼ A sin ðvtÞ gives the following eigenvalue problem:

2l

2 0 0

0 1 0

0 0 1

2

664

3

775

þ k

3 21 0

21 2 21

0 21 1

2

664

3

775

0

BB@

1

CCA

x1

x2

x3

2

664

3

775

¼

0

0

0

2

664

3

775

where l ¼ v2m=k: The eigenvalues are found by setting the characteristic equation of the determinant

equal to zero. This gives the following polynomial:

l3 2 4:5l2 þ 5l 2 1 ¼ 0

It is a simple matter to graph this polynomial and see that there are three real roots to this equation.

6.5.2 Bisection Method

One method for calculating the roots of a nonlinear function is called the bisection method. This method

finds the zeros of nonlinear functions by bracketing the zero in an interval ½a; b􀀉: The interval is chosen so

that f ðaÞ and f ðbÞ are of opposite sign. If f is a continuous function and it is positive at one endpoint, say

f ðaÞ . 0; and it is negative at the other endpoint, f ðbÞ , 0; then it has had to go through zero at some

point in the interval. It is possible that it has gone through zero more then once in the interval. From

Figure 6.7, one can see that one root is between [0, 0.5]; another is between [1, 1.5]; and a third root is in

[2.5, 3]. In order to find all three roots, the bisection method has to be used three times; one for each

interval. As an example, the MATLAB code for the first root is given below.

6.5.2.1 MATLAB Code for the Bisection Method

clear

%set an acceptable tolerance for the root

tol ¼ 10e 2 4

%endpoints of the interval

a ¼ 0

b ¼ 0:5

for i ¼ 1:20

c ¼ 0:5pða þ bÞ

if abs( f(c)) , tol, break, end

FIGURE 6.6 Three spring – mass system. (Source:

Thomson and Dahleh 1998. Theory of Vibration

Applications, 5th ed. With permission.)

Numerical Techniques 6-15

© 2005 by Taylor & Francis Group, LLC

if f ðaÞpf ðcÞ , 0

b ¼ c

else

a ¼ c

end

end

The root 0.2554 is found after ten iterations. As long as the initial interval is chosen so that the

function has different signs at the endpoints, this method will always converge. There are nonlinear root

finding methods, such as Newton’s method, which when they converge do so more quickly. However,

they may or may not converge.

6.5.2.2 MATLAB Function for Finding the Roots of a Polynomial

MATLAB has a built in root-finding method which requires the creation of a vector c that contains the

coefficients of the polynomial in descending order. The MATLAB command roots (c), produces the roots

of the polynomial. For our example, c ¼ ½1; 24:5; 5; 21􀀉: The roots of this equation are given by

l ¼ 2.8892, 1.3554, 0.2554.

6.5.2.3 Mode Shapes

The mode shapes are determined by substituting each eigenvalue into the matrix problem and

computing the corresponding eigenvector. As an example of how to compute the eigenvector for a given

eigenvalue, we will compute the eigenvector associated with l ¼ 2:8892: The first step is to substitute this

value of l into the eigenvalue equation. This gives the following problem:

2:489 21 0

21 1:745 21

0 21 0:744

2

664

3

775

x1

x2

x3

2

664

3

775

¼

0

0

0

2

664

3

775

The solution vector found by Gaussian elimination is given by

x ¼ ðx1; x2; x3Þt ¼ ð0:2992; 0:7446; 1:0Þt

6

4

2

0

−2

−4

−6

−8

−10

−12

−1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5

FIGURE 6.7 Graph of the cubic polynomial l3 2 4:5l2 þ 5l 2 1 ¼ 0: (Source: Thomson and Dahleh 1998. Theory

of Vibration Applications, 5th ed. With permission.)

6-16 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

6.5.3 Directly Calculating the Eigenvalues and Eigenvectors from the Matrix

Equation

Ideally one would like to use a computer program such as MATLAB to compute both the eigenvalues and

the eigenvectors directly from the matrix equation. In order to do this, the matrix equation needs to be

rewritten slightly. In general, the matrix equation for the normal mode vibration is given by

½2lM þ K􀀉x ¼ 0

where M and K are the mass and stiffness matrices, respectively, both are square symmetric matrices, and

l is the eigenvalue related to the natural frequency by l ¼ nv2: Premultiplying the preceding equation by

M21, we have another form of the equation:

½2lI þ A􀀉x ¼ 0

where A ¼ M21K and A is called the dynamic matrix. In general, A is not symmetric. In order to use

MATLAB to compute the eigenvectors and eigenvalues, the matrix must be symmetric. There is a

standard transformation of coordinates that converts the matrix problem into a standard eigenvalue

problem which can be solved by a computer. Assume that we have a transformation of coordinates which

has the following form:

x ¼ U 21

When this transformation is substituted into the equation, we get

½2lMU 21 þ KU 21􀀉y ¼ 0

Premultiplying the equation by U2t gives

½2lU 2t MU 21 þ U2t KU 21􀀉y ¼ 0

From this equation, one can see that if either M or K equals Ut U; then the preceding equation is reduced

to standard form, at which time it is possible to compute both the eigenvalues and eigenvectors directly.

To see how this works, let us assume that M ¼ Ut U: The equation becomes

½2lI þ U2t KU21􀀉y ¼ 0

where l ¼ v2: This equation is in standard form since U2t KU21 is symmetric.

6.5.3.1 Example

To illustrate the use of the dynamic matrix and the standard computer form, we can use MATLAB to

calculate the eigenvalues and eigenvectors for Example 6.5.1. The first step is to convert the problem into

standard form. For this example, the dynamic matrix is given by

A ¼

1:5 20:5 0

21:0 2:0 21:0

0 21:0 1:0

2

664

3

775

We can compute both the eigenvalues and the eigenvectors in MATLAB using the command

½U; D􀀉 ¼ eigðAÞ: The result of this command is two matrices, U and D: Matrix U contains the

eigenvectors as column vectors and D is a diagonal matrix which has the eigenvalues on the diagonal.

Continuing with our example, we get the following two matrices:

U ¼

0:7569 0:3031 0:2333

0:2189 20:8422 0:5808

20:6158 0:4458 0:7799

2

664

3

775

Numerical Techniques 6-17

© 2005 by Taylor & Francis Group, LLC

and

D ¼

1:3554 0 0

0 2:8892 0

0 0 0:2544

2

664

3

775

In order to use this transformation then, one has to be able to write M or K as Ut U: This is known

as the Cholesky decomposition. In MATLAB, if one has a positive definite matrix M; the matrix can

be Cholesky-decomposed with the built in function chol. The command U ¼ cholðMÞ produces an

upper triangular matrix U such that Ut U ¼ M:

Cholesky decomposition is a special case of LU factorization where L is a lower triangular matrix and U

is an upper triangular matrix. The computational procedure begins by writing the algebraic equations

that result from the following calculation:

u11 0 0

u12 u22 0

u13 u23 u33

2

664

3

775

u11 u12 u13

0 u22 u23

0 0 u33

2

664

3

775

¼

m11 m12 m13

m21 m22 m23

m31 m32 m33

2

664

3

775

If the matrix M is positive definite, then the above matrix multiplication results in six linearly

independent equations.

6.5.4 Summary of Matrix Methods

1. Bisection method for the roots of a polynomial.

2. MATLAB command roots for finding the roots of a polynomial.

3. Cholesky decomposition.

Навигация