6.2 Single-Degree-of-Freedom System

Back

The simple undamped spring – mass system oscillates

with a frequency known as the natural

frequency. The natural frequency is determined

by the spring constant and the mass in the

following way, v ¼

ffiffiffiffiffi

k=m p : This relationship is

derived by applying Newton’s second law to the

basic spring – mass system. The resulting equation

is given by mx€ ¼ 2kx: The solution to this

equation is harmonic with the frequency specified

by the expression v ¼

ffiffiffiffiffi

k=m p :

A more realistic vibration model of a simple oscillatory system includes a mass, a spring and a damper

(see Figure 6.1). For simplicity, the mass is concentrated at the center of mass of the object, the spring is

assumed to be of negligible mass, and, for the purposes of this discussion, the damping will be modeled

by viscous damping. This is described by a force proportional to the velocity, denoted by x_: In the case of

no external forcing, the system can be described by the following equation:

mx€ þ cx_ þ kx ¼ 0

where m; c; and k are constants and m is the mass, c is the damping coefficient, and k is the spring

constant. This equation can be solved analytically by assuming a solution of the form

x ¼ est

where s is a constant. Substitution into the differential equation yields the following quadratic equation:

ðms2 þ cs þ kÞest ¼ 0

This equation is satisfied for all values of t when the following quadratic equation, known as the

characteristic equation, is satisfied:

s2 þ

c

m

􀀏 􀀐

s þ

k

m ¼ 0

The characteristic equation has two roots

s1 ¼ 2

c

2m þ

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ðc=2mÞ2 2 ðk=mÞ

q

and

s1 ¼ 2

c

2m

2

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ðc=2mÞ2 2 ðk=mÞ

q

The general solution is

x ¼ A es1 t þ B es2 t

The constants A and B are determined by the initial conditions xð0Þ and x_ð0Þ: The single spring–mass

system exhibits three types of behavior, overdamped, underdamped and critically damped.

FIGURE 6.1 Single spring – mass system.

6-2 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

The overdamped state exists when k=m is larger than ðc=2mÞ2: No oscillations exist in this state. The

underdamped case is oscillatory and results when ðc=2mÞ2 is larger than k=m: The limiting case between

oscillatory and nonoscillatory motion occurs when ðc=2mÞ2 ¼ k=m: When this condition, is met, the

system is said to be critically damped.

The next level of complexity results when the system is forced harmonically. The following equation

serves as a model for this system:

mx€ þ cx_ þ kx ¼ F sin vt

The solution of this equation is found by first computing the complementary function, which is the

solution of the homogeneous equation discussed above, and then the particular solution. The details for

computing a particular solution can be found in books on differential equations or mechanical vibrations

(Thomson and Dahleh, 1998).

6.2.1 Forced Single-Degree-of-Freedom System

In general, when the oscillatory system is forced in a nonharmonic way, the resulting differential equation

cannot be solved in closed form, and numerical methods must be employed to predict the behavior of the

system. In this section, we consider two finite difference methods chosen for their simplicity. For a more

detailed discussion of numerical methods for ordinary differential equations see Atkinson (1978) and

Isaacson (1966).

The spring–mass systemthat is subjected to general forcing, Fðx; x_; tÞ; can be modeled by the following

differential equation:

mx€ þ cx_ þ kx ¼ f ðx; x_; tÞ

x0 ¼ xð0Þ

x_0 ¼ x_ð0Þ

where m; c; and k are constants and F is an arbitrary function. The two initial conditions, x0 and x_0, are

known.

In the first method, the second-order equation is integrated without change in form; in the second

method, the second-order equation is rewritten as a system of two first-order equations and then the

system of equations is integrated. Both methods approximate the first and second derivatives with the

centered difference approximation for the derivatives. Finite difference methods are based on the Taylor

expansion. The centered difference method for the first and second derivative results from a combination

of the forward and backward Taylor expansion about the point xi: To get the forward expansion, one

writes the Taylor expansion for xiþ1: Similarly, the backward expression is obtained from the Taylor

expansion about xi21: These are given by

xiþ1 ¼ xi þ hx_i þ

h2

2

x€i þ

h3

6

fflx þ · · ·

xi21 ¼ xi 2 hx_i þ

h2

2

x€i 2

h3

6

fflx þ · · ·

where h ¼ Dt: A second-order approximation is one which matches the Taylor expansion exactly up to

and including terms of order h2; that is, to determine a second-order expansion one neglects terms of

order h3 and higher. A second-order approximation for the first derivative is obtained by subtracting the

backward difference from the forward difference. The resulting centered difference approximation is

given by

x_i ¼

1

2h ðxiþ1 2 xi21Þ 2

h2

6

fflx þ · · ·

Errors result when this expression is truncated after the first term. These errors depend on h2 and the

third derivative of xi: If the error in the computed derivative is larger than order h2 it may well arise from

the neglected third derivative term.

Numerical Techniques 6-3

© 2005 by Taylor & Francis Group, LLC

The centered difference approximation for the second derivative is found by adding the forward

and backward difference expansions and ignoring terms of order h4 and higher. The resulting

approximation is

x€i ¼

1

h2 ðxi21 2 2xi þ xiþ1Þ þ

h2

12

xiv þ · · ·

The truncation error that occurs for this expression depends on h2 and the fourth derivative of xi: Both

the first and second derivative centered difference approximations are order h2: They can be used

together to create a second-order approximation to a second-order ordinary differential equation such as

that describing the spring – mass system.

6.2.1.1 Centered Difference Approximation

After substituting these two centered difference approximations into the differential equation and

rearranging terms, one gets the following recurrence relation for the single spring – mass system:

xiþ1 ¼ h2f ðxi; tiÞ 2 ð fh2 2 2mÞxi 2 ðm 2 ch=2Þxi21

with the initial conditions

x0 ¼ xð0Þ

x_0 ¼ x_ð0Þ

The recurrence relation should be used to compute all values of x from the initial condition. By letting

i ¼ 1 in the recurrence relation, we get

x2 ¼ h2f ðx1; t1Þ 2 ðkh2 2 2mÞx1 2 ðm 2 ch=2Þx0

In order to compute x2; both x0 and x1 are needed. The initial conditions provide x0: However, to start

the calculation, we need an additional equation for x1: This equation is derived by substituting i ¼ 0 into

the Taylor expansion for xiþ1; using the initial conditions and ignoring terms of order h2 and higher.

Since the centered difference approximation is of order h2; it is consistent to compute x1 with an error of

order h2: For the point i ¼ 0; the forward difference is given by

x1 ¼ x0 þ hx_0

This equation allows one to find x1 in terms of the two initial conditions, after which the recurrence

relation can be used to find all subsequent discrete values of x:

6.2.1.2 Pseudocode for Centered Difference Approximation

The following is an example of the MATLABw routine for the solution to the centered difference

approximation of the general single-DoF spring – mass equations. The function f needs to be specified in

a separate function file:

%initial conditions (index has to start at 1)

xð1Þ ¼ x0

xd ¼ x00

%specify a time step

h ¼ H

%specify the constants m; k; c:

m ¼ M

k ¼ K

c ¼ C

%compute xð2Þ from the Taylor expansion

xð2Þ ¼ xð1Þ þ h p xd

%specify the total number of steps

TEND ¼ tfinal

6-4 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

%compute the solution for all remaining times using the recurrence relation

tð1Þ ¼ h

tð2Þ ¼ 2 p h

for I ¼ 3:TEND

Xi ¼ h2 p f ðxi21; ti21Þ 2 ðk p h2 2 2 p mÞ p xi21 2 ðm 2 c p h=2Þ p xi22

TðiÞ ¼ tðI 2 1Þ þ h

end

The method just presented has ignored terms of order h2 and higher. This is known as the truncation

error. The calculation will contain other errors such as round-off error due to the loss of significant

figures. This loss is related to both the machine and the language used for the calculations. The round-off

errors are also related to the time increment h ¼ Dt in a complicated way that is beyond the scope of this

chapter. Better accuracy can be obtained by choosing a smaller Dt: However, the smaller the Dt; the larger

the number of computations needed to reach the solution in a fixed time T: The increased number of

computations affects both the total time of the calculation and the overall accuracy.

6.2.1.3 Runge – Kutta Methods

The centered difference approximation is not a self-starting method. In other words, the calculation that

determines x1 from the initial conditions does not come from the discretized equation, but rather from

the Taylor expansion directly. An alternate way to solve this equation is to use what is known as a Runge –

Kutta method. All of these methods approximate the differential equation with Taylor series expansions.

The advantage of this class of methods is that they are self-starting. The disadvantage is that they only

work on first-order equations (or systems of first-order equations). Before the Runge – Kutta method is

discussed in detail, the second-order spring – mass equation has to be written as a system of two firstorder

equations. The equation for the single-DoF system, subjected to arbitrary forcing f is

mx€ þ cx_ þ kx ¼ f ðx; x_; tÞ

It can be rewritten as the following system of first-order equations

x_ ¼ u

u_ ¼ x€ ¼

1

m ½f ðx; u; tÞ 2 cu þ kx􀀉

These equations are coupled and need to be solved together. Any Runge – Kutta method for the system

of equations requires a Taylor series expansion of x and u about xi and ui: The Taylor expansion is given

by

xiþ1 ¼ xi þ x_ih þ x€

h2

2 þ · · ·

uiþ1 ¼ ui þ u_ ih þ u€ i

h2

2 þ · · ·

As above, the time increment will be denoted by h ¼ Dt: The first-order Runge – Kutta method (also

known as Euler’s method) is obtained by retaining terms of first-order and lower, i.e., the Euler

approximation is given by

xiþ1 ¼ xi þ x_ih 2 oðhÞ

uiþ1 ¼ ui þ u_ ih 2 oðhÞ

The truncation error for Euler’s method is order h: The error depends linearly on h:

Matching more terms in the Taylor expansion generates higher order Runge – Kutta methods. The

most commonly used Runge – Kutta method is the fourth-order method that matches the Taylor

expansion up to terms of order h4: This is a significant reduction in the error. For a derivation of this

method see Cheney and Kincaid (1999).

Numerical Techniques 6-5

© 2005 by Taylor & Francis Group, LLC

For the system of first-order equations given above, the fourth-order Runge – Kutta method requires

four values of t; x; u, and, G where G ¼ 1=m½f ðx; u; tÞ 2 cu þ kx􀀉: It can be computed for each point i as

follows in Table 6.1.

Combining these quantities in the following method gives the fourth-order Runge – Kutta method:

xiþ1 ¼ xi þ h=6ðU1 þ 2U2 þ 2U3 þ U4Þ

uiþ1 ¼ ui þ h=6ðG1 þ 2G2 þ 2G3 þ G4Þ

where it is recognized that the four values of U divided by six represent the average slope dx=dt and the

four values of G divided by six result in an average of du=dt: A way to check the accuracy of a Runge –

Kutta method is to Taylor expand G1; G2; G3; and G4 and collect like terms. One will find that the above

combination produces an expansion which is exact up to order h4:

6.2.1.4 Pseudocode for the Fourth-Order Runge – Kutta Method

For simplicity, the code is given for the single first-order equation

x_ ¼ Gðt; xÞ with the initial data xð0Þ ¼ x0:

The function G should be specified in a function file.

%initial conditions (index has to start at 1)

xð1Þ ¼ x0;

%time step needs to be specified

h ¼ H;

h2 ¼ 2 p H;

%TEND is the total number of time steps

TEND ¼ Tfinal;

T ¼ h;

for I ¼ 1:TEND

G1 ¼ H p Gðt; xÞ;

G2 ¼ H p Gðt þ h2; x ¼ 0:5 p G1Þ;

G3 ¼ H p Gðt þ h2; x ¼ 0:5 p G2Þ;

G4 ¼ h p Gðt þ h; x þ G3Þ;

X ¼ x þ ðG1 þ 2 p G2 þ 2 p G3 þ G4Þ=6:0;

t ¼ ðI þ 1Þ p h;

end

6.2.1.5 Example

Solve numerically the differential equation

4x€ þ 2000x ¼ FðtÞ

with the initial conditions

x0 ¼ x_0 ¼ 0

The forcing is as shown in Figure 6.2.

TABLE 6.1 Fourth-Order Runge – Kutta Method for Spring – Mass Equation

t x X_ ¼ u X€ ¼ G

T1 ¼ ti X1 ¼ xi U1 ¼ ui G1 ¼ GðT1; X1 ; U1 Þ

T2 ¼ ti þ h=2 X2 ¼ xi þ U1 h=2 U2 ¼ ui þ G1 h=2 G2 ¼ GðT2; X2 ; U2 Þ

T3 ¼ ti þ h=2 X3 ¼ xi þ U2 h=2 U3 ¼ ui þ G2 h=2 G3 ¼ GðT3; X3 ; U3 Þ

T4 ¼ ti þh X4 ¼ xi þ U3h U4 ¼ ui þ G3h G4 ¼ GðT4; ; X4 ; U4Þ

Source: Thomson and Dahleh 1998. Theory of Vibration Applications, 5th ed. With permission.

FIGURE 6.2 The forcing function for Example 6.2.1.5.

(Source: Thomson and Dahleh 1998. Theory of Vibration

Applications, 5th ed. With permission.)

6-6 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

6.2.1.5.1 Centered Difference

Using the centered difference approximation for the second derivative, one gets the following discrete

equation

xiþ1 ¼ h2=4FðtÞ 2 500h2xi 2 xi21 þ 2xi

This equation is valid for i $ 1: From the initial conditions x0 ¼ x_0 ¼ 0; x1 is computed by the

following x1 ¼ x0 þ x_0 ¼ 0 þ 0: A time step of h ¼ 0.02 sec has been used in the calculation. The

numerical solution as compared with the exact solution is given in Table 6.2. Table 6.3 contains the

absolute errors for the two discrete solutions.

6.2.1.5.2 Runge – Kutta

In order to use the Runge – Kutta method, the second-order equation needs to be written as a system of

first-order equations.

Let u ¼ x_; then

u_ ¼ Gðx; tÞ ¼ :25 p FðtÞ 2 500x

where xð0Þ ¼ 0 and uð0Þ ¼ 0:

This is now in a form which can be directly input into a fourth-order Runge – Kutta solver.

The Runge – Kutta method is more accurate then the centered difference approximation. This can be

seen in Table 6.3, which gives the absolute value of the difference between the exact solution and the

computed solutions. This is known as the absolute error. Moreover, the Runge – Kutta method is selfstarting

and can be used for a single variable or a system of variables as in the above example. The price

that is paid is that the fourth-order Runge – Kutta method requires four function evaluations for the first

derivative for each time step. This is offset by the fact that this method has higher accuracy and has a large

stability region so one can take larger time steps.

TABLE 6.2 Solution to Example 6.2.1.5

Time t Exact Solution Central Difference Runge – Kutta

0 0 0 0

0.02 0.00492 0.00500 0.00492

0.04 0.01870 0.01900 0.01869

0.06 0.03864 0.03920 0.03862

0.08 0.06082 0.06159 0.06076

0.10 0.08086 0.08167 0.08083

0.12 0.09451 0.09541 0.09447

0.14 0.09743 0.09807 0.09741

0.16 0.08710 0.08712 0.08709

0.18 0.06356 0.06274 0.06359

0.20 0.02949 0.02782 0.02956

0.22 2 0.01005 2 0.01267 2 0.00955

0.24 2 0.04761 2 0.05063 2 0.04750

0.26 2 0.07581 2 0.07846 2 0.07571

0.28 2 0.08910 2 0.09059 2 0.08903

0.30 2 0.08486 2 0.08461 2 0.08485

0.32 2 0.06393 2 0.06171 2 0.06400

0.34 2 0.03043 2 0.02646 2 0.03056

0.36 0.00906 0.01407 0.00887

0.38 0.04677 0.05180 0.04656

0.40 0.07528 0.07916 0.07509

0.42 0.08898 0.09069 0.08886

0.44 0.08518 0.08409 0.08516

0.46 0.06436 0.06066 0.06423

0.48 0.03136 0.02511 0.03157

Numerical Techniques 6-7

© 2005 by Taylor & Francis Group, LLC

Stability is a measure of how quickly errors in the computed solution grow or decay. There are very few

numerical methods that are stable for all choices of time step. For most methods, there is a range of time

steps which produce a stable method. The fourth-order Runge – Kutta method is stable for larger values

of the time step than are the lower-order Runge – Kutta methods. A numerical method will not converge

to a solution if the time step does not produce a stable method. Often, the stability criterion places a

stricter limitation on the time step than accuracy does. For a more complete discussion of stability, see

Strang (1986). In Table 6.3, the absolute error for both methods grows but it does not grow exponentially.

Controlled error growth is the signature of a stable time step.

6.2.2 Summary of Single-Degree-of-Freedom System

* Unforced equation of motion mx€ þ cx_ þ kx ¼ 0:

* Forced equation of motion mx€ þ cx_ þ kx ¼ f ðx; x_; tÞ:

* Centered difference approximation for the first derivative x_i ¼

1

2h ðxiþ1 2 xi21Þ:

* Centered difference approximation for the second derivative x€i ¼

1

h2 ðxi21 2 2xi þ xiþ1Þ:

* Fourth-order Runge – Kutta method.