Transient analysis of structures by stiffly stable methods

Compufers & Slrucrures, Vol. 4, pp. 615-626.

TRANSIENT

Pergamon Press 1974. Printed in Great Britain

ANALYSIS OF STRUCTURES STABLE METHODS

BY STIFFLY

PAUL S. JENSEN Structural Mechanics Laboratory, Lockheed Missiles and Space Co., 3251 Hanover Street, Palo Alto, California 94304, U.S.A. Abstract-The equation of motion for a general structure can be cast as a system of second order ordinary differential equations in discretized form. The locus of the roots for the homogeneous form of the-se equations is analyzed and related to the stability region for general, stifily stable, linear multistep methods. The problem of recasting the second order system as a first order system is analyzed. noting that the final system to be solved is no larger than the original system. A discussion of error control is included and some specific formulas are given. Finally, a new family of third order formulas with a number of advantages over previously published formulas is given.

INTRODUCTION THE EQUATION

of motion for a general structure can be written in the form Mii+ Dti+Su=f

(1)

where u represents displacement and M, D and S represent mass, damping and stiffness of the structure. The forcing function f in general varies both with time and with displacement u. In nonlinear problems, the stiffness and damping can also vary with displacement u. For computational purposes, the displacement is evaluated at a finite set of points on the structure and at discrete points in time. Equation (1) is then taken as a vector equation in terms of the discrete values in time and space with differential quantities in space, which appear in the mass, damping and stiffness functions, replaced by discrete difference quantities. This is generally accomplished by either finite difference or finite element techniques. In vector form, equation (1) is commonly a ‘stiff system of ordinary differential equations in time. By that we mean that some of the components of u vary extremely rapidly with time whereas other components vary slowly with time. This is simply a result of the fact that some portions of a structure are often mechanically stiffer (relative to their masses) than other portions and, hence, tend to vibrate with higher frequency. In the literature, the term ‘stiff’ is applied to any system of ordinary differential equations with widely varying time constants and its relationship to mechanical stiffness is only incidental. The numerical integration of a stiff system of equations presents one with a dilemma, viz:. (a) if the time step is made small enough to accurately treat the stiff (rapidly varying) components, then it is excessively small for the other components which results in excessive computational effort and round-off error in the ‘slower’ components; (b) qn the other hand, if the time step is adjusted to the slowly varying components, instability results for most high order methods and, in any case, the stiff equations are not accurately solved. The challenge of such a dilemma cannot go unmet and, consequently, a large number (over 30) of papers on stiff equation systems have been published since about 1952 when the problem was apparently first attacked. 615

616

PAUL S. JENSEN

One obvious approach is to have different time steps for different portions of the system. Of course, this requires the ability to identify which are stiff and which are not and leads to rather complicated programming. Dahlquist [I] has discussed an approach similar to this for the case in which the high and low frequency components are well separated. In many cases, a great deal of accuracy in the solution of the stiff (and usua!Iy damped) components is not required. For such problems, the objective is to find an integration technique which does not go unstable for large time step (adjusted to the low frequency components). For studying these problems, the classical definition of stability is not strong enough, and so different stability criteria have been devised. In the ciassical sense (Henrici f2], p* ll), a method is stabie if, roughly speaking, it yields a uniformly bounded solution over a fixed time period for a!! step sizes in the range It,,, >h>O, for some positive (problem dependent) limit h,,,. For stiff equation systems, Dahlquist [3] introduced the concept of A-stability. A method is A-stable if for any fixed step size h and )., ReA~0, the method produces a solution {J*")Of the eqUatiOn “C-/Q

(2)

such that y,+O as n+o~. In a general sense, if a method is A-stable, then the error introduced into a solution by the method remains uniformly bounded for any step size. Dahlquist also showed t!rat if a linear multistep method is A-stable, then (i) it is implicit, and (ii) it is of order less than 3, i.e. the error is on the order of h, h2 or h3. A good treatment of three genera! A-stable methods can be found in [4]. The last result is rather discouraging in view of the rather large truncation error. It suggests that perhaps requiring A-stability is too restrictive for many problems. Pursuing this idea, Gear [S, 6f has introduced the concept of stifly stabIe methods which can be of much higher order. These methods have the properties of A-stability for equation (2) for all values of 1-h outside the kidney shaped region of the complex plane as illustrated in Fig. 1. A-stable methods correspond to stiffly stable methods (Fig. 1) with v=O. The second order stiffly stable method is A-stable. Stiffly stable, linear, multistep methods of order up to 11 have been formulated [7].

FIG. I. General form of the unstable region of the complex i. 4 plane for the stiffty stable methods

of Gear.

Transient Analysis of Structures by Stiffly Stable Methods

617

1. STABLLITY REGION OF THE STRUCTURAL EQUATION OF MOTION For anaiysis of the homogeneous ($=O) form of the discretized form of equation (l), we consider first the companion eigenvalue probiem

wfMxi,

Sx, =

i=l,...,II

(3)

where the xi are normalized so that XTMXj = S, =

1 i=j 0

i#j.

Let n x n matrix

have the eigenvectors of (3) as its columns and let A=diag (w:, . . . , a”,) be the diagonal matrix of the corresponding eigenvalues. Then we have sx== MXA and so XTSX= A.

(5)

v(t)= x-‘u(t),

(6)

Substituting a new displacement vector

into the homogeneous vector form of (1) and premultiplying by Xr yields

or, from (4) and (5) P+(XQX)t+Av=O.

(7)

A commonly accepted form of the damping matrix which shall be used for the remainder of this paper is given by D=uM+PS

(8)

where cy:and /I are problem dependent control parameters. Substitution of (8) into (7) yields the uncoupled system of ordinary differential equations i; + (af+ &I)* + Av = 0.

Thus the solution to the&h equation of (9) is given by

(9)

618

PAUL

S. JENSEN

where a, and bj are arbitrary constants and yr and 7, are the roots of

For brevity, we drop the subscript j and consider a typical exponent

(11) Letting z = (a + /302)/2 we have W2= (2z- u)/P. Substitution into (11) then yields

y= -z+iJr2-(Z-1/p)?

(12)

where r=Ji7&. Thus, we see that for all z such that the radicand of (12) is positive, the locus of y is simply an upper semi-circle of radius r centred at x= - l/p. 7 then lies on the reflected semicircle. This locus is illustrated in Fig. 2 for the complex I241plane corresponding to Fig. 1. iy

Flo. 2. Root locus for the structural equation of motion in the complex 14~ plane. Observations. As a-0, the radius r+h//3 which simply enlarges the circle until it reaches the origin. Similarly as b-0, the radius r--to3 and so the locus becomes a vertical line at x= -ha/2. Also, as @*co, the maximum root y-+ -h/j3 as seen from (12). For the stiffly stable methods of Gear [4, one canfairly easily determine approximate limits G, for thejth order method such that the numerical solution will be stable whenever the step size h
j=3,...,6

for any a~0. Similar limits J, can be found for the extended stiffly stable methods of Jain and Srivastava [I 5] which trade improved stability characteristics for slight loss in

Transient Analysis of Structures by Stiffly Stable Methods

619

accuracy. These limits are given in Table 1. Notice that these limits can be increased by specification of a minimum &in > 0. We now discuss the application of stiffly stable methods to structural problems given by equation (1). The methods to be considered generally appear as first order, linear multistep methods which must be suitably modified for the second order structural problem. TABLE 1. Approximate

limits for determining maximum stable step-size At
Jain and Strivastava

3

7.0

6.6

4

2.7

2.6

5

1.3

1.4

6

0.5

0.6

Order j of method

JJ

-

2. LINEAR MULTISTEP

METHODS

Given a system of ordinary differential equations in the form

Y(b>‘Yo,

(13)

a linear m-step method for solving (13) is of the form

i

i=l

%Y”+i=h

i i=l

n= -l,O,

Bifn+i,

1, . . .

(14)

where h is the step size, yn & y(t, +n*h), and a, #O. The initial value y0 is given and the starting values yl, . . . , y,,,_ 2 are determined by a starting procedure such as a Runge-Kutta type method [B]. The values y,,,- r, y,,, . . . are then evaluated using (14). Detailed discussions of general multistep methods can be found in [2 and 81. If pms,=O,equation (14) is an explicit equation for y”+,,, and is frequently called a predictor. When &,,#O, as in all stiffly stable methods, equation (14) is an implicit expression for y.+m and is frequently called a corrector. For practical application, equation (14) can be rewritten in simpler form by Y” -

69”= w,,

n=nz-1,

m,

...

(15)

where &&! a,

(15a)

and the ‘historical’ portion w,= - +mg: (aiy,-,+i-hpiJi.-,+l). mi-

Wb)

620

PAUL

S.

JENSEN

Substitution of (13) into (15) yields the general implicit formula Y” - WY,,

48)= w#n

n=m-1,

m,

...

(16)

which may be solved at each step by any convenient technique for such a (generally nonlinear) system. Householder [9] has accumulated quite an extensive bibliography on these techniques, Gear [lo] discusses solution by variations of simple corrector iteration y~(~*l)-~f(y~*~ f,) = w,,

i=O, 1,. . .

starting with a predicted initial vector y.“I. A discussion on the use of Newton iteration for this purpose may be found in [ll]. A general class of methods iteratively solves an equation of the form Q*y~+%fi,

i=o,

1,. . .

(17)

where y$‘) is an initial approximation of the solution of (16) and Qi and r, are, in general, functions of yn(I). For example, Newton iteration yields

ri=w ” +S(f(y n,n wf

)-J.y rn (‘j)

(184

where Ji is the Jacobian matrix (;~j(X, r,,)/%x,) evaluated at x=yi*). The initial approximation yic) is typically obtained from a predictor formula in the form of equation (14). The form (17) is also very appropriate for linear systems of equations (13) which may be expressed by j, =Jy “I-f(t)

(19)

where matrix J is constant. In this case iteration is not needed, of course, and from (I6 and 17) we note that Q=I-6.f

cm

and r = w, + Sf(f,),

(W

similar to (18). 3. TRANSFORMATION

TO A FIRST ORDER SYSTEM

Since the given problem (1) is a second order system, we consider its recasting as a first order system (13). Let us introduce an auxiliary system of equations of the form v=AMi+Bu

(21)

where the matrices A and B are arbitrary except that A is not singular. Combining (1 and 2 1) we obtain the first order system

Transient

Analysis of Structures by Stiffly Stable Methods

(E-BY)(:)+(:s -;) (:)=(AoI)*

621

(22)

Now combining (22 and 15) and dropping the time subscript n we obtain

or collecting terms, AM +SB A(D +&S)-B

-6;)

(:>=(yy

(23)

where w(i) and w(‘) are the u and v historical vectors (15b) and

It is simple to algebraically solve (23) to obtain E@)u= Mw(‘)+ 6q

(24)

and v=Aq-(d(D+&!T)-B)u

(24a)

E(~)=M+~D+c~~S

(24b)

where

is independent of the choice of A and B. Thus, it is convenient to choose A=& B=D

(25)

so that from (24a) and the second row of (22) we find v=q-GSU

(26)

and ;=i-su.

(26a)

li=(u-w(‘))/s

(26b)

Also, from (15) we have the expression

for ii. Summarizing, at each time step it is our objective to find a zero of n dimensional functional of equation (24) F(u)=E(@u-g(u) where g(u)=Mw”‘+6q

(27)

622

PAUL S. JENSEN

or, from (23a and 25) g(u) =

Mw(‘) + 8w@)+a2f.

(27a)

Hence, we see that even though we had to expand the system (1) to 2n equations because it is of second order, we nevertheless have only to solve n simultaneous (generally nonlinear) equations at each time step.

4. ERROR CONTROL Although the error estimation technique we are about to discuss does not produce rigorous bounds, it does provide reasonably good results in practice. Variations of the technique appear in a number of papers, e.g. Gear [6]. A good, recent study in this area appeared in [12]. Let y(t,) denote the true solution at time t, and y” denote the calculated solution. Assume that Yn_i=Y(tn-i), i=l, . . . , m- 1. Then the calculated solution satisfies Yn =Yw

(28)

+ %~mY(m)(rl)

for known truncation coefficient a,,, and some unknown time t,_l<(Y1< t,, where here ytk) signifies the kth derivative of y with respect to time.1 Using an m step prediction (14) one can obtain a predicted value Y,: =Y(cJ

-

(29

4nh”Y’“‘(t2)

again with known truncation coefficient b, (depending upon the formula used) and some time fn_i<
Y,I-

YL= (arn+ 43~”

(

+ *

-+Y'"'G)

“1

m

=(am+ 4n)hmy(m)(t)

0,

YW2) m > (30)

by the mean value theorem, where e lies between t1 and c2, assuming ytm) continuous. Thus, combining equations (28 and 30) we obtain the error estimate

zn= Y”

-Y(b)

= *(Y.-Y.). m

(31) m

In practice, one periodically calculates this error estimate and increases the step size if it appears to be very small, or decreases the step size if it appears to be getting too large.

Transient Analysis of Structures by Stiffly Stable Methods

5. MULTISTRP

623

FORMULAS

The multistep, stiffly stable formulas devised by Gear [5] can be expressed by m

c aIYn+i=hB3.+,-amamhmy(m)(r), i=l

m-2,.

. .,

(32)

similar to equation (14) where we take yn+ [= y(r,+,), the true value. Several sets of these coefficients are given in Table 2 for the reader’s convenience. 2. Coefficients for stiffly stable methods

TABLE

in

3

4

5

6

7

anI

219

3122

12/125

10/137

6011029

6

12

Pnr

2

al

1

a2

-2 9

-4

3

a3

3

-18

36

11

a4

-16

as

-48 25

a6

60

60

-12 75 -200 300 -300 137

a7

10 -72 225 -400 450 -360 147

For the purposes of error control, some ‘compatible’ predictor formulas in the form

i~lak+r=h/l,-1~;,,-l+ab6,hmy(m~(~), m=2,. . . (33) have been devised. By compatible, we mean that truncation coefficient b, of (33) is positive as required by (29). These are given in Table 3. TABLE3. Compatible

predictor coefficients

m

3

4

5

6

7

bm

l/3

l/4

l/5

l/6

117

Pm-1

2

6

12

60

60

a’1

-1

1

a'2

0

a'3

1

a'4

a’s a'6 a'7

-1

-6 3 2

3 6

-20

10

-120

-18

60

3

65 12

-2 15 -50 100 -150 77 10

624

PAUL

6. A NEW FAMILY

S. JENSEN

OF THZRD ORDER METHODS

Jain and Srivastava f15] showed how the use of one extra ‘historical’ value of the solution could be used to improve either the stability or the accuracy of the Gear type, stiffly stable formulas. As can be seen in Table 1, Section 1, this approach results in a negligible improvement in the stability characteristics for structural applications. The use of an extra historical value of the solution derivative for this purpose; however, yields a substantial improvement in the stability characteristics of the third order formulas. A comprehensive investigation of higher order formulas has not yet been made. Consider the general multistep formula (14) with m=4 where we set PI =f13 =O and _

-

- _

1

0

-2 9

2 -4

-18 11

_

2

+c

0

-1

6 -

1 _ _

(34)

It turns out that this formula is stable for all finite c>O; however, as c--too the magnitude of the largest subdominant and dominant roots of the ‘characteristic’ polynomials f?(X)=

fl

aixi-’

and

respectively approach 1. Thus, the method is unstable in the limit. The values of these magnitudes rP(c) and r,,(c) for the polynomials p(x) and a(x) are given by and

TV = j2c + 7 f ,/‘(2c + 3)’ - 48 f/(4c + 22)

(35)

_I,(C) = ,/c/(6 + c).

(36)

Several values for these characteristic roots are given in Table 4. As these magnitudes approach 1, the solution error due to the spurious roots becomes more and more noticeable for certain problems. For ~220, results appear to be quite satisfactory for a wide range of problems. The stability region for the new formula has the appearance shown in Fig. 1. A comparison of the key factors (see Fig. 1 and Table 1) of the new formula (c = 10) with those of the Gear and the Jain and Srivastava stable formulas are given in Table 5. Of major interest to structural appli~tion is the factors max (A~~~)and 6 appearing in Table 5. For damping coefficients or=0 and /? are given, the new method allow a maximum stable step size of about 3.8 (sy26/7) times that of the other two formulas. On the other hand, if fi=O and a is given, absolute stability for the new formula is guaranteed for step sizes down to about l/5 (R OX@95/0*052)of that for the other two formulas.

Transient Analysis of Structures by Stiffly Stable Methods TABLE4. Several values of the characteristic c

rd4

625

roots

k(c)

0

0*426*

5

0667

0.674

10

0.789

0.791

15

0.845

0.845

20

0.877

0.877

50

0.945

0.945

100

0.971

o-971

0

*complex root

TABLE5. Complex stability region comparison of third order stiWy stable methods (see Fig. 1 and Table 1)

0.084

l-94

7-o

0.052

144

6-6

c=lO

0.0095

1.35

26.4

c=20

0.0044

1.27

46.2

Gear formula Jain and Srivastava Stable Formula New formula

The remaining parameter of considerable importance for multistep formulas is the error coefficient a,,, of (32). For our new family of third order formulas, this coefficient is given by a3(c)=(c+9)/(12c+66).

(37)

It is pleasantly surprising that not only does the stability improve with increasing c, but also the truncation error coefficient a,(c) decreases. The largest value, a,(O), is that of the Gear formula. 7. RESULTS The analysis technique described in this paper has been implemented in a computer program called STUDY (STrUctural Dynamics). It has been applied to non-structural as well as structural problems including a problem describing the behavior of a laser gas. For structural applications it has been used together with a linear analysis finite element program developed by Wilson [13] and a nonlinear analysis finite difference program STAGS [14]. It is also being used in conjunction with a new analysis program for studying ‘n-body’ structures consisting of relatively rigid bodies connected together relatively loosely.

626

PAULS. JENSEN

To date, extensive comparison of the stiffly stable approach with other commonly used techniques in structural analysis has not been made. The linear and nonlinear test problems that have been run have produced very encouraging results from the standpoint of accuracy and stability. Acknowledgements-The author wishes to thank Messrs. B. 0. Almroth and F. A. Brogan and Dr. T. L.

Geers for many lively discussions on the character of transient structural problems and for the many numerical experiments carried out by them. The author also thanks Mr. Theo Winarske for his capable programming of much of the STUDY computer program.

REFERENCES Ill G. DAHLQUIST.A numerical method for some ordinary differential equations with large Lipschitz constants. ZFZP Congress, Edinburgh (1968). [21 PETERHENRICI,Error Propagation for Difference Methods. Wiley, New York (1963). A special stability criterion for linear multistep methods. BIT 3, 22-43 (1963). 131 G. DAHLQUIST, 141W. LINIGERand R. A. WILLOUGHBY,Efficient integration methods for stiff systems of ordinary differential equations. Siam .I. Num. Anal 7 (l), 47-66 (1970). PI C. W. GEAR,Numerical integration of stiff ordinary differential equations. Report No. 221, Dept. of Computer Sciences, University of Illinoi!, Urbana, January (1967). [6] C. W. GEAR, The automatic integration of stiff ordinary differential equations. ZFZP Congress, Edinburgh (1968). [fl M. K. JAINand V. K. SRIVASTAVA, High order stiffly stable methods for ordinary differential equations. Report No. 394; Department of Computer Sciences, University of Illinois, Urbana, April (1970). [8] PETERHENRICI,Discrete Variable Methods in Ordinary Difirential Equations. Wiley, New York (1962). [9] A. S. HOUSEHOLDER, KWIC index for the numerical treatment of nonlinear equations. 0RNL-4595, UC-32-Mathematics and Computers, Oak Ridge National Laboratory, November, (1970). [lo] C. W. GEAR, The automatic integration of ordinary differential equations. Comm. ACM 14 (3), 176-179 (1971). [ll] I. W. SANDBERCI and H. SCHICKMAN, Numerical integration of systems of stiff nonlinear differential equations. Bell Systems Tech. J. 47, 51l-527 (1968). [12] F. T. KRoc+H,Changing stepsize in the integration of differential equations using modified divided differences. Section 914, Tech. Memo. 312, California Institute of Technology Jet Propulsion Laboratory, Pasadena, California, October (1972). [13] E. L. WILSON,A computer program for the dynamic analysis of underground structures. Contact DACA 39-67-C-0020 Report, Waterways Experimental Station, U.S. Army Corps. of Engineers, Berkeley, California, January (1968). [14] B. 0. ALMR~TH,F. A. BROGAN,E. MELLERand F. ZELE,User’s manual for the STAGS computer code. LMSGD266611, Lockheed Missiles and Space Co., Sunnyvale, California, 7 April (1972). [15] M. K. JAINand V. K. SIUVASTAVA, Optimal stiffly stable methods for ordinary differential equations. Report No. 402, Department of Computer Sciences, University of Illinois, Urbana, June (1970). (Received 9 March 1973)