A faster modified newton-raphson iteration

COMPUTER METHODS IN APPLIED MECHANICS @ NORTH-HOLLAND PUBLISHING COMPANY

A FASTER MODIFIED

AND ENGINEERING

NEWTON-RAPHSON

20 (1979) 267-278

ITERATION

M.A. CRISFIELD Transport and Road Research ~borato~,

Depar~ent

of the Environment,

Received 12 October 1978 Revised manuscript received 29 December

Crowt~o~e,

Berks. RGll

6AU,

U.K.

1978

The paper describes an accelerated modified Newton-Raphson iteration in which the iterative deflection change is a scalar times the previous iterative change plus a further scalar times the usual unaccelerated change. These scalars are automatically recalculated at each iteration. They are related to inner products involving the iterative deflections and the present and past out-of-balance force vectors. The extra computation for each iteration is negligble, the only penalty being the storage of two extra vectors, The method is based on a secant approach and leads to a significant reduction in the required number of iterations. Examples are presented in which the method is applied to the n,oniinear analysis of thin plates and shells. The technique is used in conjunction with the finite eiement method.

Introduction The modified Newton-Raphson iterative procedure [l] is frequently used for nonlinear structural analysis with finite elements. It is often used in conjunction with an incremental loading procedure. In comparison with the Newton-Raphson procedure, the modified method has two distinct advantages: (i) the structure stiffness matrix need only be formed and factorised once per increment, and, (ii) provided the stresses are in (stable) equilibrium at the beginning of the increment, the tangent stiffness matrix will always be positive definite. In comparison with the Newton-Raphson method, the modified approach suffers the disadvantage of a slower convergence rate, particularly in the later iterations. A number of procedures have been advocated for the acceleration of the modified NewtonRaphson method [2], [3]. Th e “alpha-constant” approach [2] effectively introduces a diagonal accelerating matrix at every second iteration [43. Jennings [3] proposes a modification to Aitken’s acceleration [S] whereby Aitken’s diagonal accelerating .rnat~x is replaced by a scalar to be applied equally to all terms in the iterative vector. Jennings argues that this leads to greater stability. In the present approach, not only is a scalar applied to the usual iterative deflection, but also a proportion of the previous iterative deflection is added into the present vector. In contrast to Jennings’ method, the present accelerations are automatically introduced at every iteration. The method was derived as a result of an examination of the variable metric method f6]-[15]. A brief description of the variable metric method is therefore required in order to introduce the basis of the present approach. However, the technique is first described without derivation. The description relates to structural analysis although, potentially, the field of application is much wider.

M.A. Crisfkld, Faster modified Newton-Raphson

268

iteration

1. The accelerated modified Newton-Raphson method The deflections pi are up-dated by an iterative change Si such that pi+1

=

pi

+

4.

(1)

If K, is an approximate tangent stiffness matrix (Hessian), the iterative deflection given by the standard modified Newton-Raphson procedure is gi = - K,‘gi,

(2)

where gi is the gradient (of the total potential energy) or out-of-balance force vector (with the opposite sign to that usually adopted)- gi is a function of pi. The accelerated iterative deflection is Si = e&1+ J;ii,

(3)

where e, and fi are scalars given by fi = - ajbi,

(4

ei =fi(l-cJbi)-1, where ai, bi and ci are scalar inner products given by ai = S:-lgi-l, bi = S:-l(gi - gi_1) = S:-lyi, Ci =

The procedure

(5)

;i:yim is started (i = 0) with a tangential incremental

solution in which

fcl= 1, e. = 0,

(6)

go = -4,

where Aq is the incremental load vector; also p. is the deflection vector at the end of the previous increment. In comparison with the usual modified Newton-Raphson method the present accelerated technique requires the storage of the two extra vectors S,_1 and gi_1. It can be shown that eq. (3) is a form of secant relationship by multiplying both sides of eq. (3) by yi (or gi - gi-A i.e. S:ri = - Si_,gi. In one dimension

the standard secant method (fig. 1) approximates

(7) the Hessian (stiffness) by

M.A. Crisfield, Faster modifid

Newton-Raphson

269

iteration /

/ / / /

-Yi

/ /

gi-l

/ / /

_P

-hi_1 pi

Fig. 1. Secant

c

pi-l

relationship.

YJGi-1SO that

Eq. (7) can be considered as an n-dimensional form of eq. (8). Further insight can be gained by considering the behaviour of the technique when applied to the minimisation of some quadratic governing function f (in structural terms the total potential energy). In such a case, if the approximate stiffness matrix K, (eq. (2)) is replaced by the identity matrix, the present method is related [15] to the conjugate gradient method [16], [17]. In particular, if eq. (3) is replaced by Si = CYi(fZ?iSi-1 -figi)

(9)

and (Y~is chosen to minimise f(pi+,), the method becomes identical [15] to the conjugate gradient technique. In such a case, the previous “line search” for ai- ensures that S:-*gi = 0,

(10)

and eq. (7) becomes S:ri = ai_16:K&-l = 0,

(11)

where K, is the Hessian (stiffness matrix) related to the governing quadratic function f. Eq. (11) is the basis of the conjugate gradient method. It ensures that Si and Si-1 are orthogonal with respect to K,. When using the accelerated modified Newton-Raphson method, ai- is replaced by unity so

270

M.A. Crisjield, Faster modified Newton-Raphson

iteration

that in general the exact minimum is not located along the direction &1. Some compensation is made in the next iteration by the inclusion of the -S:-lgi term on the right-hand side of eq. (7). A clearer insight into the proposed method can be gained by studying the variable metric method from which the present method was derived.

2. The variable metric method The variable metric method [5]-[15] has been applied to structural analysis by a number of workers [12]-[15]. It is an iterative procedure in which changes are made not only to the displacement vector but also to the flexibility or stiffness matrices. The technique has had little application in conjunction with the finite element method because in its early form [6]-[7] the procedure involved up-dating approximations to the flexibility matrix (inverse Hessian) rather than to the stiffness matrix (Hessian). The storage savings introduced by the banded form of the stiffness matrix and its Cholesky factors were therefore lost. For reasons of numerical stability recent work on the variable metric method [ll] has concentrated on directly up-dating the Cholesky factors of the Hessian (stiffness matrix). Unfortunately, the procedure is still not directly applicable to the finite element method. Although the final (at equilibrium) stiffness matrix is banded, at earlier stages in the iteration the approximating matrix contains nonzero terms outside the band. Consequently, in the present approach the iterative vector is assumed to be related to that generated by the variable metric method although the Hessian (stiffness matrix) is not up-dated. The derivation of the iterative vector requires a description of the variable metric method, which is given in structural terms. For equilibrium, the principle of stationary potential energy requires that g=+

(12)

,

where 45 is the total potential energy, which is a function of the displacements p and the loads q and g is the gradient or out-of-balance force vector (with opposite sign to that usually adopted), which is generally nonzero for some approximate or guessed displacement vector p0 (subscript o means old). It is assumed that a better approximation is given by Pn = PO + 60

(subscript n means new).

(13)

A Taylor’s expansion of eq. (12) gives

where K, is the matrix dg/dp of second partial derivatives of 0 calculated at p0 (i.e. the tangent stiffness matrix). Setting eq. (15) to zero and considering only linear terms in SO gives

M.A. Crisfaeld, Faster modijied Newton-Raphson

the standard Newton-Raphson

271

iteration:

So = - K,‘g,.

Using this approach,

iteration

(16)

a further iteration is

(17)

6, = + K,'gn,

where K, = ilgldp at pn. In both the variable metric and modified Newton-Raphson methods the K, that is used with eq. (16) is generally not the true Hessian (tangent stiffness matrix) at p. but is some approximation to the true Hessian. The variable metric method aims to provide an approximation KN to K, (for use with eq. (17)) that is an improvement on K, but is not found by directly forming the true Hessian (K,) at pn. Consider an exact quadratic function f = a - q’p + $p’K,p

W-3)

with gradient SL=g@)=$=

- q +

0.

(19)

Applying a small variation 6 to p gives

g”=g(p+@=

go + y = - q + K,p + K,6

(20)

so that

y = K,6. In the variable metric method the new approximating that it should also satisfy eq. (21) i.e. Y =

G&.

(21)

matrix KN is required to resemble K, in

(22)

Consequently, the method is based on a secant relationship. One of the simplest and most successful up-dating formulas satisfying eq. (22) is (see [8], [9]) (23) It is easily verified (using eq. (16)) that eq. (23) does satisfy eq. (22). Replacing K, in eq. (17) by KN gives 6, = - K$g,,

(24)

212

M.A.

Crisfield, Faster modified Newton-Raphson

iteration

or (25)

Eq. (25) (see [ll]) may be verified by multiplying both sides by KN from eq. (23). In the variable metric method eq. (23) is used to up-date the Cholesky factors of K,, and eq. (24) is used to find the new iterative displacement vector S,. In the present procedure eq. (25) is used directly and K, is not up-dated. In this way the advantage of a banded stiffness matrix, which is associated with the modified Newton-Raphson method, is combined with the advantage of faster convergence, associated with the variable metric method. If ; = -K,‘g,

(26)

and subscripts u and n in eqs. (25) are replaced with i - 1 and i, respectively, eq. (25) can be replaced by eqs. (3)-(5). Consequently, the proposed modified Newton-Raphson method coincides with the variable metric approach when the approximate Hessian K, (eq. (2)) used in the modified Newton-Raphson method, coincides with the approximate Hessian K, (eq. (26)) used for the previous iteration of the variable metric method. Furthermore, the two methods are only strictly related when the previous iteration of the modified Newton-Raphson method is unaccelerated. Nonetheless, it is proposed that the accelerated procedure of eqs. (3)-(5) should be applied at every iteration. This proposal is justified by the secant relationship of eq. (7) and by the success of the algorithm when subjected to numerical experiments. 2.1. Damping The proposed iterative procedure relates directly to a variable metric method without line searches [S], [lo]. In the original variable metric formulations, line searches were always included, and eq. [16] is replaced by

60= - KJcgo,

(27)

where (Y, is chosen to minimise @(PO+ aJo). The line searches are dispensed with in most recent work, and (Y, is set to unity. However, occasionally, an unsatisfactory iteration will lead to an insufficient reduction (or possibly even an increase) in the potential energy. From eqs. (18) and (19) the change in energy may be approximated by

Approximating

K, by KN and using eq. (22),

A satisfactory iteration may be defined [8] as one for which

M.A.

Crisfield, Faster modified Newton-Raphson

iteration

273

(30)

A@ -=I -/_&go

(note that SAg, should be negative). In the present work A@ is approximated using eq. (29) and p is set to 0.001. As a further check the norm of the residual should be reducing, or not increasing too much. In particular, it is here required that

Ik”ll< l.lll&ll-

(31)

where ]lg]]is a scaled gradient norm which will be defined. If eqs. (30) and (31) are satisfied, the iteration is accepted. If not, (Y,is reduced below unity and a new gradient is calculated. A crude but reasonably effective procedure for reducing LY,is to set it to a

0

=o.qy.

(32)

”

This damping procedure was applied to both the original modified Newton-Raphson method and the accelerated technique. Despite such damping, divergence or very slow convergence may occur if the degree of nonlinearity is too severe for the adopted increment size. Consequently, variable increment sizes should be chosen which relate to both the degree of nonlinearity and the power of the adopted iterative procedure. This selection of increment size should ideally be automatic [18]. However, in the present work the increment sizes were chosen by judgment in an attempt to achieve a reasonable convergence rate with the original modified Newton-Raphson method. Some damping was required with this technique, but none with the more powerful accelerated iterative method. 2.2. Convergence The adopted

convergence

criteria are

Il4l< 4lFll

(33)

lkll < ~2max(11~11~ ll~ll>~

(34)

and

where ]] ]]represents the Euclidean norm, and the iterative (6) and total (p) deflections are scaled by and

&=*8j

pi

=

(35)

*pj,

where K/j is the diagonal tangent stiffness coefficient at the beginning of the increment. gradients g, applied forces q and reactions r are scaled by gi =

gjl-7

In the following examples

qj =

l1 and

qjl*7

3 =

41%.

e2 are each set to 10P4.

The

(36)

274

M.A.

Crisfield, Faster modified Newton-Raphson

iteration

3. Application 3.1. Finite element formulation The accelerated modified Newton-Raphson procedure has been incorporated in an existing finite element program for the large-deflection elastoplastic analysis of thin plates and shallow shells [19], [20]. Although preliminary results including plasticity are encouraging, the present results all relate to elastic plates and shells. The program is based on a Lagrangean formulation for moderately large deflections. The elements are rectangular with quadratic shape functions (eight nodes/element) being used for the in-plane displacement field, while the restricted quartic nonconforming shape functions [21] (which pass the patch test) are used for the out-of-plane deflections. The element stiffness matrices (Hessians) and out-of-balance force vectors (gradients) are all calculated using reduced two-point Gaussian integration, thereby minimising the amount of self-straining inherent in the adopted shape functions [20]. Smoothed bilinear slope functions [20] are used in order to save computer time. 3.2. Simply supported plate subject to lateral pressure The load is plotted against the central deflection in fig. 2. A quarter of the plate was analysed using a five-by-five element mesh. This idealisation leads to 204 (free) degrees of freedom. 300

q = Uniform

pressure

E = Young’s

modulus

t = Plate thickness wc = Central deflection a = length = breadth

of plate

Y = Poisson’s ratio = 0.316

l (i,j) 200

= Present FEM

i = no of accelerated iterations

i = no

of standard

excluding tangential solution

iterations

0

2.0

1.0 Central deflection

3.0

ratio welt

Fig. 2. Simply supported square plate under uniform pressure loading (edges restrained from in-plane movement).

M.A. Chkjield, Faster modijied Newton-Raphson

275

iteration

The results are in excellent agreement with Levy’s classical solution [22]. The number of iterations (given in parenthesis beside each point on the graph) does not include the tangential ‘incremental’ solution. Clearly, the accelerated iterative procedure introduces a very considerable saving in computer time. The reduction in the number of iterations is a function of the increment size and the degree of nonlinearity. 3.3. Clamped cylindrical shell subject to uniform pressure The load is plotted

against the central deflection in fig. 3. A quarter of the plate was analysed using a five-by-five element mesh. This idealisation leads to 195 (free) degrees of freedom. In the previous example the response always involves a hardening structure. For the shell however the response changes from softening to hardening. Again, the accelerated iterative .0020

Range of other FEM solutions(23-26’

~~‘X”’

.0015

(4.7’~

65)

“E 2 : ; :I

,001 0

E b ‘c 5

.0005 4

(i,jl

present FEM

i = noof j = noof

accelerated standard

f = some damped E = 3105

N/mm*,

R = 2540mm.

0

excluding tangential

v = 0.3, t = 3.175mm

L = B = 508mm

I

I 2.0

cylindrical

solution

t

iterations

1 .o Central

Fig. 3. Fully clamped

iterations

iterations

deflection

shell under

- 3.0

ratio we/t

uniform

pressure

loading.

MA.

276

Crisjield, Faster modified Newton-Raphson

iteration

procedure leads to a considerable saving in computer time. The saving is particularly pronounced near the inflection point. No damping was required when the accelerated iterative procedure was used, while damping was required on two increments with the standard modified Newton-Raphson method. In particular, two damped iterations were required at 4 = 0.0015 N/m m* and three at 4 = 0.00155 N/m’. Without such damping the number of standard iterations required at these load levels would have increased to 35 and 65, respectively. 3.4. Hinged cylindrical shell subject to a central point load By changing the boundary conditions, thickness and loading (fig. 4) the response of the cylindrical shell can be altered [23] to exhibit a “snap-through type of instability” (fig. 4). As with the previous cylindrical shell, symmetry was invoked and a quarter of the shell was analysed using a five-by-five element mesh with 234 (free) degrees of freedom. The response was traced numerically by adopting displacement [27] rather than load control at the central node. As with the previous examples, the introduction of the accelerated iterative procedure led to a significant improvement in the convergence characteristics.

3.2

n Sabir and Lock23 + 2.E

.T

l

wc

R

2.4

(i.j) present FEM

I = no of accelerated iterations j = no of standard iterations

t

= some damped

E = 3103

N/mm

R = 2540mm.

excluding tangential

iterations

, v = 0.3,

t = 12.7mm

L = 8 = 508mm

0.8

0

5

10

Fig. 4. Hinged

20

15 Central

cylindrical

deflection

shell under

wc

25

mm

a central

point load.

SOlutlOn

M.A. Crisfield, Faster modified Newton-Raphson

iteration

277

4. Conclusions The accelerated modified Newton-Raphson procedure is fully automatic. It is derived from the variable metric method, but unlike that method, it can be easily applied to nonlinear finite element analysis. The extra cost of each iteration is negligible, and the only penalty is the storage of one vector containing the previous iterative (deflection) change and another containing the previous out-of-balance force vector (gradient). For the examples considered, the procedure led to a very considerable saving in overall computer time due to the reduced number of iterations. At every increment level the convergence characteristics were significantly improved.

References O.C. Zienkiewicz, The finite element method, 3d ed. (McGraw-Hill, 1977) 454. G.C. Nayak and O.C. Zienkiewicz, Note on the “alpha’‘-constant stiffness method for the analysis of non-linear problems, Int. J. Numer. Meths. Eng. 4 (1972) 579-582. the convergence of matrix iterative processes, J. Inst. Math. Appls. 8 (1971) 99-110. I31 A. Jennings, Accelerating [41 A.K. Basu, Letter to the editor, Int. J. Numer. Meths. Eng. (1974) 152. PI A.C. Aitken, On the iterative solution of a system of linear equations, Proc. Roy. Sot. Edinburgh 63 (1950) 52-60. (Argonne Nat. Lab. Report ANL-5990 Rev., 1959). PI W.C. Davidon, Variable metric method for minimisation descent method for minimisation, Comp. J. 6 (1963) [71 R. Fletcher and M.J.D. Powell, A rapidly convergent 163-168. PI R. Fletcher, A new approach to variable metric algorithms, Comp. J. 13 (1970) 317-322. of a class of double-rank minimisation 2: the new algorithm, J. Inst. Math. [91 C.G. Broyden, The convergence Appls. 6 (1970) 222-231. minimisation, Math. [lOI S. Oren, Self-scaling variable metric algorithms without line search for unconstrained Comp. 17 (1973) 873-885. Pll K.W. Brodlie, An assessment of two approaches to variable metric methods, Math. Prog. 12 (1977) 344-355. in structural analysis by direct energy minimisation, A.I.A.A. J. 6 [121 R.L. Fox and E.L. Stanton, Developments (1968) 1036-1042. [I31 D.J. Dawe, A finite deflection analysis of shallow arches by the discrete element method, Int. J. Numer. Meths. Eng. 3 (1971) 529-552. Int. J. Mech. Sci. 19 (1977) 725-743. [I41 G.H. Little, Rapid analysis of plate collapse by live-energy minimisation, M.A. Crisfield, Iterative solution procedures for linear and nonlinear structural analysis. Transport and Road [I51 Research Laboratory Report LR 900 (1979). [W M. Hestenes and E. Stiefel, Method of conjugate gradients for solving linear systems, J. Res. Nat. Bur. Stands. 49 (1952) 409-436. by conjugate gradients, Comp. J. 7 (1964) 149154. [I71 R. Fletcher and C.M. Reeves, Function minimisation and instability problems using the current [W P.G. Bergan and T. H. Soreide, Solution of large displacement stiffness parameter, Finite Elements in Non-linear Mechanics 2 (Tapir, 1978) 647-670. in the non-linear analysis of thin plated structures. Finite iI91 M.A. Crisfield and R.S. Puthli, Approximations Elements in Non-linear Mechanics 1 (Tapir, 1978) 373392. for thin steel plates, in: J. Robinson (ed.), PO1 M.A. Crisfield, Combined material and geometric non-linearity World Congress on Finite Element Methods in Structural Mechanics, Bournemouth, Oct. 1975, v. 2,10.110.25. [21] O.C. Zienkiewicz and Y.K. Cheung, The finite element method for analysis of elastic isotropic and orthotropic slabs, Proc. Inst. Civ. Eng. 28 (1964) 471488. [22] S. Levy, Bending of rectangular plates with large deflection (NACA Tech. Note 846, 1942). [23] A.B. Sabir and A.C. Lock, The application of finite elements to the large deflection geometrically non-linear

[ll PI

278

[24] [25] [26] [27]

M.A.

Crisfield, Faster modified Newton-Raphson

iteration

behaviour of cylindrical shells, in: C.A. Brebbia and H. Tottenham (eds.), Variational methods in engineering 2 (Southampton Univ. Press, 1972) 7/66-7/75. A.Q. Khan, A.A. Mufti and P.I. Harris, Post-buckling of thin plates and shells, in: C. A. Brebbia and H. Tottenham, (eds.), Variational methods in engineering (Southampton Univ. Press, 1972) 7/54-7/65. N. Gass and B. Tabarrok, Large deflection analysis of plates and cylindrical shells by mixed finite element method. Int. J. Numer. Meths. Eng. 10 (1976) 731-736. T.H. Soreide, Collapse behaviour of stiffened plates using alternative finite element formulations (Inst. Statikk, Univ. Trondheim, Norway, Report No 77-3, June 1977). J.H. Argyris, Continua and discontinua. Proceedings of 1st Conference on Matrix Methods for Structural Mechanics (Wright Patterson Air Force Base, Ohio, 1965) 11-189.

Any views expressed in this paper are not necessarily those of the Department of the Environment or of the Department of Transport. Extracts from the text may be reproduced, except for commercial purposes, provided the source is acknowledged.