## J. N. Reddy

Print publication date: 2004

Print ISBN-13: 9780198525295

Published to Oxford Scholarship Online: January 2010

DOI: 10.1093/acprof:oso/9780198525295.001.0001

Show Summary Details
Page of

PRINTED FROM OXFORD SCHOLARSHIP ONLINE (www.oxfordscholarship.com). (c) Copyright Oxford University Press, 2018. All Rights Reserved. Under the terms of the licence agreement, an individual user may print out a PDF of a single chapter of a monograph in OSO for personal use (for details see www.oxfordscholarship.com/page/privacy-policy). Subscriber: null; date: 23 October 2018

# (p.439) Appendix A2 Solution Procedures for Nonlinear Algebraic Equations

Source:
An Introduction to Nonlinear Finite Element Analysis
Publisher:
Oxford University Press

(p.439) Appendix A2

Solution Procedures for Nonlinear Algebraic Equations

# A2.1 Introduction

Finite element formulations of nonlinear differential equations lead to nonlinear algebraic equations for each element of the finite element mesh. The element equation is of the form,

(A2.1.1)
where
(A2.1.2)
(A2.1.2) When [K e] is independent of {u e}, the matrix coefficients can be evaluated for all elements, and the assembled equations can be solved for the global nodal values {u} after imposing boundary conditions. When [K e] depends on the unknown solution vector {u e}, the matrix coefficients cannot be evaluated. If we can find an approximation to {u e}, say {u e}1, then [K e({u e}1)] can be evaluated and assembled. This amounts to linearizing the nonlinear equations, Eq. (A2.1.1). Then a next approximation to the solution can be obtained by solving the assembled equations,
(A2.1.3)
This procedure can be repeated until the approximate solution comes close to the actual solution in some measure. Such a procedure is called an iterative procedure.

(p.440) Here we discuss the following commonly used iterative procedures:

1. 1. The Picard Iteration (or Direct Iteration) Method

2. 2. The Newton–Raphson Iteration Method

3. 3. The Riks Method

The details of these methods are discussed next, with the aid of a single nonlinear equation. For a more complete presentation of these methods, the reader may consult the references at the end of this appendix.

Consider the nonlinear equation,

(A2.1.4)
where u is the solution to be determined, K(u) is a known function of u, F is the known ‘force’, and R is the residual
(A2.1.5)
A plot of the equilibrium path, R(u, F) = 0, is shown in Figure A2.1.1. For any value u 1,K(u 1) denotes the secant of the curve at u = u 1, and $( ∂ R ∂ u ) | u 1$ denotes the tangent of the curve at u = u 1.

# A2.2 Picard Iteration Method

In the Picard iteration, also known as the direct iteration method, we begin with an initial guess for u, say u (0), (u (0) = 0 in Figure A2.1.1) and determine a first approximation to u by solving the equation

(A2.2.1)

Figure A2.1.1 Typical force–displacement curve.

(p.441) u (1)u, and a second approximation for u is sought by using the last approximation to evaluate K
(A2.2.2)
This procedure is continued until the difference between two consecutive approximations of u differ by a preselected value. Thus, the algorithm and criterion for convergence may be written as
(A2.2.3)
(A2.2.4)
where ʥ denotes the convergence tolerance and r denotes the iteration number.

A geometric interpretation of the procedure described above is illustrated in Figure A2.2.1(a) for an initial guess of u (0) = 0. At the beginning of iteration r, the secant of the curve R(u) = 0 is found at the point u = u (r – 1) and the solution u (r) is computed using Eq. (A2.2.3). Figure A2.2.1(a) shows the convergence to the true solution u c, whereas Figure A2.2.1(b) shows a possible divergence of the algorithm. Thus, the success of the algorithm depends on the nature of the nonlinear curve R(u) = 0, the initial guess, and the load increment.

In the direct iteration method discussed above, the secant is evaluated at each iteration and inverted to obtain the next approximate solution. This can be computationally very expensive when the number of algebraic equations to be solved is large, that is, when K is a matrix and [K]–1 is its inverse. When K has a linear portion, and in most problems of interest to us it does, an alternative direct iteration algorithm can be formulated. Let

(A2.2.5)
where K L and K N(u) are the linear and nonlinear parts of K. Note that K Lis the slope at u = 0 of the curve R(u) = 0. Then we can write
(A2.2.6)
This scheme involves evaluating the nonlinear part K N at each iteration, which is computationally less expensive when compared to evaluating [K(u r – 1))] and inverting it. The inversion of K L is required only once, and it should be saved for subsequent use. The criterion in Eq. (A2.2.4) can be used to check for convergence.

Geometrically, the alternative direct iteration algorithm uses the initial slope, that is, slope at the origin of the curve for all iterations, while updating the effective force

(p.442) at each iteration (see Figure A2.2.2). The rate of convergence of this algorithm, if at all it converges, is slower than that in Eq. (A2.2.3).

Figure A2.2.1 Direct iteration scheme. (a) Case of convergence. (b) Case of divergence.

(p.443)

Figure A2.2.2 Modified direct iteration scheme.

As applied to the assembly of finite element equations (A2.1.1), the two algorithms take the following forms:

(A2.2.7)
(A2.2.8)
(A2.2.9)
The global matrices [K], [K l] and [KNL] are assembled from the corresponding element matrices:
(A2.2.10)
where
(A2.2.11)

The rate of convergence of the iterative procedure can be accelerated, in certain type of nonlinear behavior, by a relaxation procedure in which a weighted average of the last two solutions is used to evaluate [K e] (or $[ K N L e ]$): (p.444) $K e ( { u ¯ e } )$, where ${ u ¯ e } = γ { u e } ( r − 2 ) + ( 1 − γ ) { u e } ( r − 1 )$, and 0 ≤ γ‎ ≤ 1 is called the relaxation or acceleration parameter. In this case Eqs. (A2.2.7) and (A2.2.8) take the form,

(A2.2.12)
(A2.2.13)
(A2.2.14)
The actual value of γ‎ varies from problem to problem.

The computational algorithm of the direct iteration is summarized below. At each load level follow the steps:

1. 1. Compute element matrices [K e] and {F e} (for transient problems, [K e] and {F e} are to be replaced by $[ K ^ e ]$ and ${ F ^ e }$) using the solution {u}(r–1) from the previous iteration (of current load and/or time). For the first iteration of the subsequent load steps, use {u c}, the converged solution of the last load step. For the first load step use {u c} = 0, provided [K e]–1 exists, to compute the linear solution.

2. 2. Assemble the element matrices [K e] and {F e} (or $[ K ^ e ]$ and ${ F ^ e }$).

3. 3. Apply the boundary conditions on the assembled set of equations.

4. 4. Solve the assembled equations.

5. 5. Check for convergence using Eq. (A2.2.9).

6. 6a. If the convergence criterion is satisfied, increase the load to next level, initialize the counter on iterations, and repeat Steps 1–5.

7. 6b. If the convergence criterion is not satisfied, check if the maximum number of allowable iterations is exceeded. If yes, terminate the computation printing a message to that effect. If the number of iterations did not exceed the maximum allowed, update {u}(r – 1) and {u}(r) and repeat Steps 1–5.

# A2.3 Newton–Raphson Iteration Method

Suppose that we know solution of Eq. (A2.1.1) at (r – 1)st iteration and interested in seeking solution at the rth iteration. We expand R(u) about the known solution u (r – 1) in Taylor’s series,

(A2.3.1)
where δ‎u is the increment,
(A2.3.2)
(p.445) Assuming that the second- and higher-order terms in δ‎u are negligible, we can write Eq. (A2.3.1) as
(A2.3.2)
where K T is the slope (tangent) of the curve R(u) at u (r – 1):
(A2.3.3)
The residual or imbalance force, R(u (r – 1)) is gradually reduced to zero if the procedure converges. Equation (A2.3.2) gives the increment of u at the rth iteration so that the total solution is
(A2.3.4)
The iteration is continued until a convergence criterion, say Eq. (A2.2.4), is satisfied. Other convergence criteria include checking the magnitude of the imbalance force.

A geometrical interpretation of the Newton–Raphson procedure is shown in Figure A2.3.1(a). For most problems, the method has faster convergence characteristics. Figure A2.3.1(b) illustrates possible divergence of the iterative procedure for certain problems.

The Newton–Raphson method requires that the tangent K T be computed at each iteration. This can be very expensive when many degrees of freedom are involved. A modified Newton–Raphson technique involves, for a fixed load step, either keeping K T fixed while updating the imbalance force at each iteration (see Figure A2.3.2) or updating K T only at each preselected number of iterations while updating the imbalance force at each iteration. There are several other modifications of the procedure.

The Newton–Raphson and modified Newton–Raphson procedures take the following forms when applied to the assembly of element equations (A2.3.2): Newton–Raphson Procedure

(A2.3.5)
Modified Newton–Raphson Procedure
(A2.3.6)
where [K T] is the tangent matrix
(A2.3.7a)
(p.446) and
(A2.3.7b)
(A2.3.8)
and ${ U ¯ }$ is the solution at the beginning of the current load step, and
(A2.3.9)
The relaxation procedure described in Eq. (A2.3.9) can be used to accelerate convergence for certain nonlinear problems.

Figure A2.3.1 The Newton–Raphson scheme. (a) A case of convergence. (b) A case of divergence.

(p.447)

Figure A2.3.2 Modified Newton-Raphson scheme.

Since only the increment of the solution is computed in each iteration of the Newton–Raphson iteration, the incremental equations in (A2.3.5) and (A2.3.6) are subject to homogeneous form of specified essential boundary conditions of the problem. Thus after the first iteration of the first load step, any specified non-zero values of {U} should be set to zero so that in the subsequent iterations and loads, {δ‎u} is subjected to homogeneous boundary conditions.

For each load step, the following computations are required for the Newton– Raphson or modified Newton–Raphson procedure:

1. 1. Evaluate element matrices [K e] and {F e} (or $[ K ^ e ]$ and ${ F ^ e }$ for the transient problems), and compute $[ K T e ]$ and {R e} using Eqs. (A2.3.7a,b) for an element.

2. 2. Assemble element matrices $[ K T e ]$ and {R e}; for the modified Newton– Raphson iteration procedure, save either assembled [K t] or its inverse for use in subsequent iterations.

3. 3. Apply the boundary conditions on the assembled set of equations. Note: Set the specified boundary conditions on {U} to zero after Step 3 in the first iteration of the first load step.

4. 4. Solve the assembled equations.

5. 5. Update the solution vector using Eq. (A2.3.9).

6. (p.448)
7. 6. Check for convergence.

8. 7a. If the convergence criterion is satisfied, increase the load, initialize the iteration counter, and repeat Steps 1–6. For the modified Newton– Raphson iteration, compute {F e} – [K]{U e}(r – 1) in Step 1 and go to Step 2.

9. 7b. If the convergence criterion is not satisfied, check if the maximum number of iterations allowed is exceeded. If it is, terminate the computation by printing a message. If the maximum allowable number of iterations is not exceeded, go to Step 1.

In order to reduce the number of operations per iteration, in the modified Newton–Raphson method the same system matrices are used for several iterations. These matrices are updated only at the beginning of each load step or only when the convergence rate becomes poor. The modified Newton– Raphson method may require more iterations to reach a new equilibrium point.

# A2.4 Riks and Modified Riks Schemes

The Newton–Raphson method and its modifications are often used to trace nonlinear solution paths. However, the Newton–Raphson method fails to trace the nonlinear equilibrium path through the limit point (see Figure A2.4.1), because in the vicinity of a limit point the tangent matrix [K T] becomes singular, and the iteration procedure diverges. Riks [1] and Wempner [2] suggested a procedure to predict the nonlinear equilibrium path through limit points. The method provides the Newton–Raphson method and its modifications with a technique to control progress along the equilibrium path. The theoretical development of this method and its modification can be found in [36]. In the modified Riks method (see [36]) the load increment for each load step is considered to be an unknown (see [3]) and solved as a part of the solution.

The basic idea of the Riks technique can be described for a single nonlinear equation as follows (see Figure A2.4.2): The length Δ‎s of the tangent to the current equilibrium point is prescribed, and the new point is found as the intersection of the normal to the tangent with the equilibrium path [see Figure A2.4.2(a)]. Then iteration is performed along the normal toward the new equilibrium point, as illustrated in Figure A2.4.2(a). Crisfield [4] suggested using a circular arc in place of the normal [see Figure A2.4.2(b)]. The center of the circle is at the current equilibrium point and Δ‎s is its radius. For multidimensional problems the normal and circular arcs become a plane and sphere, respectively, Crisfield [4] updated the tangent stiffness matrix only at the beginning of each load increment (i.e. modified Newton–Raphson method). In the present study we describe the modified Riks method due to Crisfield [4].

(p.449)

Figure A2.4.1 Load–deflection curve with limit points.

We wish to solve Eq. (A2.1.4) for u as a function of the source term F. If F is independent of the geometry, we can write it as

(A2.4.1)
where λ‎ is a scalar, called load parameter, which is considered as an unknown parameter. Equation (A2.1.5) becomes
(A2.4.2)
Now suppose that the solution $( u n ( r − 1 ) , λ n ( r − 1 ) )$ at (r – 1)st iteration of the nth load step is known and we wish to determine the solution $( u n ( r ) , λ n ( r ) )$ at the rth iteration. Expanding R, which is now a function of λ‎ and u, in Taylor’s series about the known solution, we have,
Omitting the higher-order terms involving the increments $δ λ n ( r )$ and $δ u n ( r )$, we obtain
(A2.4.3)

(p.450)

Figure A2.4.2 The Riks scheme. (a) Normal plane method. (b) Circuar arc method.

where K t =∂R/∂u is the tangent matrix [see Eq. (A2.3.3)]. The incremental solution at the current iteration of the nth load step is given by

(A2.4.4a)
where $δ u n ( r )$ is the usual increment in displacement due to known out-of-balance force vector $R n ( r − 1 )$ with known $λ n ( r − 1 )$ and K T is the tangent at the beginning of the current load increment (i.e. Modified Newton-Raphson method is used)
(A2.4.4b)
(p.451) and $δ u ^ n$ is the tangential solution (see Figure A2.4.3)
(A2.4.4c)
Note that K T is evaluated using the converged solution u n – 1 of the last load step,
(A2.4.5)
and $δ u ^ n$ is computed at the beginning of each load step.

The solution at the rth iteration of the current load step is given by

(A2.4.6a)
(A2.4.6b)

For the very first iteration of the first load step, we assume u = u 0, and a value for the incremental load parameter $δ λ 1 0$, and solve the equation

(A2.4.7)
and compute $δ u 1 ( 0 ) = δ λ 1 0 ⋅ δ u ^ 1$.

A2.4.3 Modified Riks scheme.

(p.452) Select the arc length Δ‎s to be the length of the vector

(A2.4.8)
Substituting for $Δ u n ( r )$ from Eqs. (A2.4.6b) [and $δ u n ( r )$ from Eq. (A2.4.4a)] we obtain the following quadratic equation for the increment in the load parameter, $δ λ n ( r )$:
(A2.4.9a)
where
(A2.4.9b)
Let us denote the roots of this quadratic equation as $δ λ n 1 ( r )$ and $δ λ n 2 ( r )$. To avoid “tracing back” the equilibrium path (i.e. going back on the known equilibrium path), we require the angle between the incremental solution vectors at two consecutive iterations, $Δ u n ( r − 1 )$ and $Δ u n ( r )$ be positive. Corresponding to the two roots, $δ λ n 1 ( r )$ and $δ λ n 2 ( r )$, there correspond to two values of $Δ u n ( r )$, denoted $Δ u n 1 ( r )$ and $Δ u n 2 ( r )$. The root that gives the positive angle is the one we select from $( δ λ n 1 ( r ) , δ λ n 2 ( r ) )$. The “angle” is defined to be product of the vector $Δ u n ( r − 1 )$ and $Δ u n ( r )$. Then we check to see which one of the following two products is positive:
(A2.4.10a)
If both roots are positive, then we select the ones closest to the linear solution,
(A2.4.10b)

The first arc length Δ‎s is computed using Eqs. (A2.4.7) and (A2.4.8):

(A2.4.11)
To control the number of iterations taken to converge in the subsequent load increments, δ‎s can be scaled,
(A2.4.12)
where Δ‎s n – 1 is the arc length used in the last iteration of the (n – 1)st load step, I d is the number of desired iterations (usually < 5) and I 0 is the number of iterations required for convergence in the previous step. Equation (p.453) (A2.4.12) will automatically give small arc lengths in the areas of the most severe nonlinearity and longer lengths when the response is linear or nearly linear. To avoid convergence of the solution at a higher equilibrium path, maximum arc lengths should be specified.

For all load steps after the first, the initial incremental load parameter $δ λ n ( 0 )$ is calculated from

(A2.4.13)
The plus sign is for continuing the load increment in the same direction as the previous load step and the negative sign is to reverse the load step. The sign follows that of the previous increment unless the value of determinant of the tangent matrix has changed in sign.

The modified Riks procedure described above for a single equation can be extended to the finite element equations in (A2.1.1). We introduce a load parameter λ‎ as an additional dependent variable,

(A2.4.14)
In writing Eq. (A2.4.14), loads are assumed to be independent of the deformation. The assembled equations associated with Eq. (A2.1.1) become
(A2.4.15)
The residual vector {R} is now considered as a function of both {U} and λ‎.

Now suppose that the solution $( { U } n ( r − 1 ) , λ n ( r − 1 ) )$ at the (r – 1)st iteration of the nth load step is known. Expanding {R} in Taylor’s series about $( { U } n ( r − 1 ) , λ n ( r − 1 ) )$, we have

where the subscript ‘n’ is omitted for brevity. Omitting second- and higher-order terms involving δλ‎(r) and $δ U j ( r )$, we can write
or
(A2.4.16a)
(A2.4.16b)
where {F} is the load vector, ${ R } n ( r − 1 )$ the unbalanced force vector at iteration (r – 1),
(A2.4.17a)
(p.454) and $δ λ n ( r )$ the load increment [given by Eq. (A2.4.10)], and
(A2.4.17b)

For the first iteration of any load step, $δ λ n ( 0 )$ is computed from [see Eq. (A2.4.13)]:

(A2.4.18)
where [see Eq. (A2.4.12)]
(A2.4.19a)
and $( Δ s ^ ) n$ is the arc length computed using the relation
(A2.4.19b)
{Δ‎U}n – 1 being the converged solution increment of the previous load step. For the first iteration of the first load step, we use:
(A2.4.20a)
(A2.4.20b)
where $δ λ 1 0$ is an assumed load increment and {U 0} is an assumed solution vector (often we assume $δ λ 1 0 = 1$ and {U 0} = {0}).

The solution increment is updated using

(A2.4.21)
and the total solution at the current load step is given by
(A2.4.22)

The constants in Eq. (A2.4.9b) are computed using

(A2.4.23)

The computational algorithm of the modified Riks method is summarized next.

(p.455) First iteration of first load step

1. (i) Choose a load increment $δ λ 1 0$ (say, $δ λ 1 0 = 1$) and solution vector {U}0 (say, {U} = {0}).

2. (ii) Form element matrices $[ K T e ]$ and

3. (iii) Assemble element matrices.

4. (iv) Apply the boundary conditions.

5. (v) Solve for ${ δ U ^ } 1$ and ${ δ U ¯ } 1 ( 1 )$ using Eqs. (A2.4.17a,b).

6. (vi) Compute the solution increment [see Eqs. (A2.4.16) and (A2.4.21)] and update the solution

7. (vii) Update the load increment $λ 1 ( 1 ) = δ λ 1 0$

8. (viii) Compute the arc length [see Eq. (A2.4.20a)]

9. (ix) Go to Step 9.

First iteration of any load step except the first

1. 1. Calculate the system matrices [K e], $[ K T e ]$ and {F e}.

2. 2. Assemble the element matrices.

3. 3. Apply the boundary conditions.

4. 4. Compute the tangential solution

5. 5. Compute the initial incremental load parameter $δ λ n ( 0 )$ by Eq. (A2.4.18):

(A2.4.24)

6. 6. Compute the incremental solution using Eq. (A2.4.17a):

7. 7. Update the total solution vector and load parameter:

(A2.4.25)

8. (p.456)
9. 8. Check for convergence [see Eq. (A2.4.28)]. If convergence is achieved, go to step 15 below. If not, continue with the 2nd iteration by going to step 9.

## The rth iteration of any load step (r = 2, 3, …)

1. 9. Update the external load vector

(A2.4.26)

2. 10. Update the system matrices (skip forming of [K T] for modified Newton-Raphson iteration).

3. 11. Solve for ${ δ U ¯ } n ( r )$ and ${ δ U ^ } n$ from the two sets of equations in (A2.4.17a,b); for the modified Newton–Raphson method Eq. (A2.4.17b) need not be resolved.

4. 12. Compute the incremental load parameter $δ λ [ = δ λ n ( r ) ]$ from the following quadratic equation:

and Δ‎s is the arc length of the current load step. Two solutions δλ‎1 and δλ‎2 of this quadratic equation are used to compute two corresponding vectors ${ Δ U } n 1 ( r )$ and ${ Δ U } n 2 ( r )$ The δλ‎ that gives positive value to the product ${ Δ U } n ( r − 1 ) ⋅ { Δ U } n ( r )$ is selected. If both δλ‎1 and δλ‎2 give positive values of the product, we use the one giving the smallest value of (–a 3/a 2).

5. 13. Compute the correction to the solution vector

and update the incremental solution vector, the total solution vector and the load parameter:
(A2.4.27)

6. 14. Repeat steps 9–13 until the following convergence criterion is satisfied:

(A2.4.28)

7. 15. Adjust the arc length for the subsequent load steps by [see Eq. (A2.4.9a)] $Δ s = Δ s ^ ( I d / I 0 )$.

8. 16. Start a new load step by returning to Step 1.

(p.457) 1. Riks, E., “The Application of Newton’s Method to the Problem of Elastic Stability,” Journal of Applied Mechanics, 39, 1060–1066 (1972).

2. Wempner, G. A., “Discrete Approximations Related to Nonlinear Theories of Solids,” International Journal of Solids and Structures, 7, 1581–1599 (1971).

3. Batoz, J. L. and Dhatt, G.,“Incremental Displacement Algorithms for Non-linear Problems,” International Journal for Numerical Methods in Engineering, 14, 1262–1267 (1979).

4. Crisfield, M. A., “A Fast Incremental/Iterative Solution Procedure that Handles ‘Snap-Through’,” Computers and Structures, 13, 55–62 (1981).

5. Bergan, P. G., Horrigmoe, G., Krakeland, B. and Soreide, T. H., “Solution Techniques for Non-linear Finite Element Problems,” International Journal for Numerical Methods in Engineering, 12, 1677–1696 (1978)

6. Padovan, J. and Tovichakchaikul, S., “Self-Adaptive Predictor–Corrector Algorithms for Static Nonlinear Structural Analysis,” Computers & Structures, 15, 365–377 (1982).

7. Newton, I., “De analysis per aequationes infinitas (1690),” in The Mathematical Papers of Isaac Newton, Vol. 11 (1667–1670), D. T. Whiteside (ed.), Cambridge University Press, Cambridge, UK, 207–247 (1968).

8. Raphson, J., “Analysis Aequationum universalis seu ad aequationes algebraicas resolvendas methodus generalis et expedita, ex nove infinitarum serierum doctrina deducta ac demonstrata,” London (1690) (original in British Library, London).

9. Cajori, F., “Historical Note on the Newton–Raphson Method of Approximation,” American Mathematical Monthly, 18, 29–32 (1911).

10. Bicanic, N. and Johnson, K. H., “Who was ‘-Raphson’?,” International Journal for Numerical Methods in Engineering, 148–152 (1979). (p.458)