(p.439) Appendix A2 Solution Procedures for Nonlinear Algebraic Equations
(p.439) Appendix A2 Solution Procedures for Nonlinear Algebraic Equations
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,
(p.440) Here we discuss the following commonly used iterative procedures:

1. The Picard Iteration (or Direct Iteration) Method

2. The Newton–Raphson Iteration Method

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.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
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
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
As applied to the assembly of finite element equations (A2.1.1), the two algorithms take the following forms:
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}_{NL}^{e}]$): (p.444) ${K}^{e}(\{{\overline{u}}^{e}\})$, where $\left\{{\overline{u}}^{e}\right\}=\gamma {\left\{{u}^{e}\right\}}^{(r2)}+(1\gamma ){\left\{{u}^{e}\right\}}^{(r1)}$, and 0 ≤ γ ≤ 1 is called the relaxation or acceleration parameter. In this case Eqs. (A2.2.7) and (A2.2.8) take the form,
The computational algorithm of the direct iteration is summarized below. At each load level follow the steps:

1. Compute element matrices [K ^{e}] and {F ^{e}} (for transient problems, [K ^{e}] and {F ^{e}} are to be replaced by $[{\widehat{K}}^{e}]$ and $\left\{{\widehat{F}}^{e}\right\}$) 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. Assemble the element matrices [K ^{e}] and {F ^{e}} (or $[{\widehat{K}}^{e}]$ and $\left\{{\widehat{F}}^{e}\right\}$).

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

4. Solve the assembled equations.

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

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

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,
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
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 nonzero 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. Evaluate element matrices [K ^{e}] and {F ^{e}} (or $[{\widehat{K}}^{e}]$ and $\left\{{\widehat{F}}^{e}\right\}$ for the transient problems), and compute $[{K}_{T}^{e}]$ and {R ^{e}} using Eqs. (A2.3.7a,b) for an element.

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. 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. Solve the assembled equations.

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

6. Check for convergence.

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.

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 [3–6]. In the modified Riks method (see [3–6]) 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].
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
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
The solution at the rth iteration of the current load step is given by
For the very first iteration of the first load step, we assume u = u _{0}, and a value for the incremental load parameter $\delta {\lambda}_{1}^{0}$, and solve the equation
(p.452) Select the arc length Δs to be the length of the vector
The first arc length Δs is computed using Eqs. (A2.4.7) and (A2.4.8):
For all load steps after the first, the initial incremental load parameter $\delta {\lambda}_{n}^{(0)}$ is calculated from
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,
Now suppose that the solution $({\left\{U\right\}}_{n}^{(r1)},{\lambda}_{n}^{(r1)})$ at the (r – 1)st iteration of the nth load step is known. Expanding {R} in Taylor’s series about $({\left\{U\right\}}_{n}^{(r1)},{\lambda}_{n}^{(r1)})$, we have
For the first iteration of any load step, $\delta {\lambda}_{n}^{(0)}$ is computed from [see Eq. (A2.4.13)]:
The solution increment is updated using
The constants in Eq. (A2.4.9b) are computed using
The computational algorithm of the modified Riks method is summarized next.
(p.455) First iteration of first load step

(i) Choose a load increment $\delta {\lambda}_{1}^{0}$ (say, $\delta {\lambda}_{1}^{0}=1$) and solution vector {U}_{0} (say, {U} = {0}).

(ii) Form element matrices $[{K}_{T}^{e}]$ and

(iii) Assemble element matrices.

(iv) Apply the boundary conditions.

(v) Solve for ${\left\{\delta \widehat{U}\right\}}_{1}$ and ${\left\{\delta \overline{U}\right\}}_{1}^{(1)}$ using Eqs. (A2.4.17a,b).

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

(vii) Update the load increment ${\lambda}_{1}^{(1)}=\delta {\lambda}_{1}^{0}$

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

(ix) Go to Step 9.
First iteration of any load step except the first

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

2. Assemble the element matrices.

3. Apply the boundary conditions.

4. Compute the tangential solution

5. Compute the initial incremental load parameter $\delta {\lambda}_{n}^{(0)}$ by Eq. (A2.4.18):
(A2.4.24) 
6. Compute the incremental solution using Eq. (A2.4.17a):

7. Update the total solution vector and load parameter:
(A2.4.25)
(p.456)

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, …)

9. Update the external load vector
(A2.4.26) 
10. Update the system matrices (skip forming of [K _{T}] for modified NewtonRaphson iteration).

11. Solve for ${\left\{\delta \overline{U}\right\}}_{n}^{(r)}$ and ${\left\{\delta \widehat{U}\right\}}_{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.

12. Compute the incremental load parameter $\delta \lambda [=\delta {\lambda}_{n}^{(r)}]$ from the following quadratic equation:

13. Compute the correction to the solution vector
(A2.4.27) 
14. Repeat steps 9–13 until the following convergence criterion is satisfied:
(A2.4.28) 
15. Adjust the arc length for the subsequent load steps by [see Eq. (A2.4.9a)] $\Delta s=\Delta \widehat{s}({I}_{d}/{I}_{0})$.

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 Nonlinear Problems,” International Journal for Numerical Methods in Engineering, 14, 1262–1267 (1979).
4. Crisfield, M. A., “A Fast Incremental/Iterative Solution Procedure that Handles ‘SnapThrough’,” Computers and Structures, 13, 55–62 (1981).
5. Bergan, P. G., Horrigmoe, G., Krakeland, B. and Soreide, T. H., “Solution Techniques for Nonlinear Finite Element Problems,” International Journal for Numerical Methods in Engineering, 12, 1677–1696 (1978)
6. Padovan, J. and Tovichakchaikul, S., “SelfAdaptive 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)