12.3 Nonlinear Iteration Requests
NONLINEAR_ITERATIONS
NONLINEAR_ITERATIONS iteration_type = list , etc...
Specify the parameters for nonlinear iterations to use for all equations in all sub-domains/element groups belonging to any given solution stagger.
|
|
|
|
|
|
|
(1) |
Iteration_type |
list |
[*] |
Nonlinear Iteration procedure: |
|
|
modified_Newton_Raphson |
|
Modified Newton-Raphson iterations |
|
|
|
Newton_Raphson |
|
|
Newton-Raphson iterations |
|
|
quasi_Newton_BFGS |
|
|
Quasi-Newton iterations (BFGS update) |
|
|
quasi_Newton_Strang_BFGS |
|
Quasi-Newton iterations (BFGS product) |
|
|
|
quasi_Newton_Broyden |
|
|
Quasi-Newton iterations (Broyden update) |
|
|
|
|
Linear iterations |
|
|
|
|
|
|
|
|
(2) |
Max_number_of_iterations |
integer |
[0] |
|
|
|
Min_number_of_iterations |
integer |
[0] |
Minimum number of iterations. |
|
|
|
|
|
|
|
|
Convergence_check |
list |
[on] |
Convergence check |
|
|
on / off |
|
|
|
|
|
|
|
|
|
|
|
Convergence_tol_rhs |
real |
[1.E-3] |
Convergence tolerance for residual |
|
|
Convergence_tol_sol |
real |
[1.E-3] |
Convergence tolerance for solution |
|
(7) |
Convergence_norm L2_norm infinite_norm other |
list |
[*] |
Convergence norm |
|
|
|
|
|
|
|
|
Relative_tolerance |
real |
[1.E-3] |
Relative tolerance
Rtol |
|
|
|
|
|
|
|
|
Absolute_tolerance |
real |
[0.0] |
Absolute tolerance
Atol |
|
|
|
|
|
|
|
(3) |
Jacobian_matrix on / off |
list |
[off] |
Jacobian matrix |
|
|
|
|
|
(cont’d) |
(cont’d)
|
|||||||
|
|
|
|
|
|
||
|
(3) |
Jacobian_derivative |
list |
(*) |
Jacobian derivative option |
||
|
|
backward |
|
|
|
||
|
|
forward |
|
|
|
||
|
|
central |
|
|
|
||
|
(4) |
Reform_lhs_step |
integer |
[0] |
Step number at which left-hand side |
||
|
|
|
|
|
(Stiffness) reforms are to be initiated. |
||
|
|
|
|
|
|
||
|
(4) |
Reform_lhs_freq |
integer |
[0] |
Left-hand side (Stiffness) reform frequency |
||
|
(4) |
Reform_lhs_iter |
integer |
[0] |
Left-hand side (Stiffness) reform frequency during iterations |
|||
|
|
|
|
|
|
|||
|
|
Max_update_vectors |
integer |
[0] |
Maximum number of updating |
|||
|
|
|
|
|
vectors (only applicable to Quasi- |
|||
|
|
|
|
|
|
|||
|
(5) |
Line_search_type |
list |
[none] |
Line search requests |
|||
|
|
Strang |
|
|
|
|||
|
|
line_minimization |
|
|
|
|||
|
|
backtracking |
|
|
|
|||
|
|
backtrack_line_minimization |
|
|
||||
|
|
backtrack_Strang |
|
|
|
|||
|
|
|
|
|
|
|||
|
|
Line_search_freq |
integer |
[1] |
Line search frequency |
|||
|
|
|
|
|
||||
|
|
|
|
|
||||
|
Scaling on / off |
list |
[off] |
|||||
EXAMPLE
Iteration_type = Quasi_Newton_BFGS / # Quasi-Newton BFGS iterations request
Max_number_of_iterations = 10 / # up to 10 iterations per time step
Convergence_tol_rhs = 1.E-6 / # specified convergence tolerance for residual
Convergence_tol_sol = 1.E-6 / # specified convergence tolerance for solution
Reform_lhs_freq = 2 / # reform lhs every 2 time step
Max_update_vectors = 5 / # maximum number of BFGS vectors
Line_search_type = Strang / # perform line search with Strang's algorithm
Line_search_freq = 1 # line search every time step
Notes /
(1) Given the solution un at time tn and a time increment 
, the object is to find the solution vector u
at time 
, which satisfies equilibrium. We write the resulting system of equations to
be solved at each time/load step as:
(1)
where
the unknowns u are the nodal values
at time tn+1, and the
equations express a balance between external (fn+1) and internal
(n(u)) forces. Both r and u are vectors, and the residual r is a system of nonlinear functionals of the solution u and of the parameters un
and . This system is
solved for u by performing a
linearization via a truncated
(2)
where
and
(3)
with
. The matrix of the
first derivatives 
in Eq. 2 is the
consistent (tangent) Jacobian matrix (e.g., 
for second-order
systems). The linear system of equations to be solved (Eq. 2) can be written
as:
(4)
· If the system of equations (Eq. 1) is linear, the expansion in Eq. 2 is exact, and the solution at step n+1 is obtained as:
(5)
·
If the system of equations (Eq. 1) is nonlinear,
it is solved for u by performing a
succession of linearizations, leading to
(6)
where the solution increment p(i) is given by:
(7)
(Notes / (cont'd)
u(i+1) and u(i) are approximations of u at iterations i and i+1, respectively. Denote
(8)
(9)
The linear system of equations (Eq.6) can be written as:
(10)
and the new solution is obtained as:
(11)
Use
of the consistent Jacobian J(i) as the left-hand side yields in
general the fastest convergence provided the initial guess is within the radius
of convergence of the
(12)
An effective representation of the
tangent (Jacobian) matrix is required in order to accurately capture the
solution of the nonlinear problem. This
can of course be accomplished through a strict
Various Newton-like iterative procedures and quasi-Newton updates are available as listed.
- Newton-Raphson Iterations: the tangent matrix is recomputed (and assembled if applicable) at every time step and at every iteration. If the consistent Jacobian derivative is returned by all the subdomains/element groups belonging to the given
solution
stagger, then the method reduces to a true
(Notes / (cont'd)
- Modified Newton-Raphson Iterations: the iterations are performed with the original matrix available at the beginning of the time/load step. The matrix may be updated by requesting a reform at the beginning of the time/load step (Reform_lhs_freq command).
- Quasi-Newton Iterations: both BFGS and Broyden rank two updates are available. The number of updating vectors to be stored is controlled by the command: Max_update_vectors. The number of updating vectors must be less or equal to the number of requested maximum number of iterations.
Two implementations of the BFGS rank two update are available. The Strang_BFGS implementation [12] is restricted to symmetric positive definite matrices, whereas the other BFGS implementation follows its original (and equivalent) form but avoids computing square roots which could be troublesome in the indefinite case. The Broyden update [3] is the most general and is not restricted to symmetric matrices. However, if applied to symmetric matrices, it may not preserve the original symmetry of the system of equations.
Quasi-Newton Methods
is the solution of
(13)
where 
(Eq. 11). The quasi-Newton or secant condition is that J(i+1) satisfy the following equation:
(14)
where 
This is the
generalization to multidimension of the one-dimensional secant approximation to
the derivative 
. However, Eq.
14 does not determine J(i+1)
uniquely in more than one dimension.
Many different auxiliary conditions to nail down J(i+1) have been explored, but one of the
best-performing algorithms in practice results from Broyden's formula [3]. Another update formula was devised by
Broyden, Fletcher, Goldfarb and Shanno and is referred to as the BFGS update
(see e.g. [5]).
(Notes / (cont'd)
Broyden Update Formula: The formula is based on the idea of getting J(i+1) by making the least change to J(i) consistent with the secant equation (Eq. 14). Broyden showed that the resulting formula is :
(15)
The formula can be expressed in terms of the inverse by using the Sherman-Morrison formula [14, 15] to invert Eq. 15 analytically as:
(16)
BFGS Formula: It is most conveniently written directly in terms of [J(i)]-1 rather than J(i), and has the form [12]:
(17)
where
v(i)
and w(i) are update vectors.
Note that the product structure of the BFGS update formula implies that
[J(i+1)]-1
inherits the symmetry (if applicable) of [J(i)]-1 since the two
new factors are transposes of one another.
Also, if [J(i)]-1 is positive-definite and 
, then [J(i+1)]-1 is also positive-definite. Multiplying out the right-hand side of Eq. 17
shows that [J(i+1)]-1 is the sum of [J(i)]-1 and two
rank-one matrices. For this reason, the
BFGS update is called a "rank-two update." Formula for the explicit definitions for the
update vectors v(i) and w(i) can be found
in [12] as:
(18)
(19)
where
(20)
(Notes / (cont'd)
Note that the line search parameter
s(i) enters the definition of 
. This factor arises
when a line search procedure is incorporated within the nonlinear solution
strategy [see Note 4]. Otherwise, s(i) = 1.0 by default.
The above definition of the BFGS update formula is referred to as the Strang_BFGS and is restricted to symmetric positive definite matrices. For problems in which the matrix is not positive definite, another implementation of the BFGS update is available which follows its original (and equivalent) form but avoids computing the square root in Eq. 20 which can be troublesome in the indefinite case. The update formula in that case is defined as:
(21)
where
(22)
Implementation Issues: The use of inverse matrices in the update formula (Eqs. 16, 17, and 21) is purely formal. The actual computation of the inverse is never carried out. Only means of replicating its action upon any vector is required. With J(i) unaltered, the only additional storage required consists of two sets of update vectors (e.g. {v(i) , w(i)}i=1, k for the Strang_BFGS). Typically, a strategy is adopted which limits the maximum number of
updates k to some practical value, e.g. kmax = 10, and is a user-specified parameter value Max_update_vectors, restricted to be less than 20. The default is min (max_number_iter, 20).
Attaining this limit is viewed as a signal to the nonlinear solution driver that the matrix should be reformed. Depending upon the situation, the update procedure is either reinitialized and the set of Quasi-Newton update vectors overwritten as new pairs are generated, or the new update vectors are accumulated over the old ones by sliding down the set of vectors loosing the oldest ones.
(2) Up to the prescribed maximum number of iterations will be performed, unless convergence is detected earlier, i.e. if both:
(Notes / (cont'd)
where r(u(i)) = f - n(u(i)) = residual, u(i) = solution at iteration i, || . || = denotes the Euclidean norm and q is a contraction factor:
where tol_rhs and tol_sol are user-specified convergence tolerances for the residual and the solution, respectively (see Command).
(3) Numerical differentiation is used to compute the consistent Jacobian effective stiffness matrix.
(4) In order to speed up convergence, reform of the effective stiffness matrix may be requested to be performed with a certain time step frequency. Note that by default, for Newton-Raphson iterations a reform is performed at every time/load step for every iteration.
(5) Line searches may be used to speed up convergence in nonlinear analyses. Line search iterations involve a single free parameter s which is determined by means of a scalar-valued objective function or metric. The strategy is then to move to a new point u(i) along the direction of the Newton-like step p(i) (not necessarily all the way) as:
such that the scalar-valued function has "decreased" sufficiently. For simplicity in the notation, we drop subscript (step number) and superscript (iteration number) in the following.
Line Searches with Strang Option: The function employed in this case [12] is the inner product of the solution increment with the recomputed residual, viz.,
The goal is then to find the root of:
g(s) = 0
and a method of successive approximations is used to obtain the solution to this equation. Trial values of s are selected based upon an interpolation/extrapolation procedure that Dahlquist and Bjork [4] refer to as the "Illinois Algorithm." Additional logic has been added to increase the robustness of the nonlinear iterations. Note that each evaluation of g(s)
(Notes / (cont'd)
requires the computation of a new residual corresponding to a new trial update. Although the line search iteration is formally motivated by finding the root of g(s) = 0, in practice seeking this condition is too computationally intensive. More typically, a relaxed convergence criterion is used:
where the convergence tolerance
may be of the order of 0.5. In
our implementation the value 
has been adopted. Once the line search iteration has converged,
the last trial update of the solution variables is accepted as the update for
the ongoing Newton-like iteration.
In order to increase the robustness and/or efficiency of the Newton-like algorithm, a control over the step length p may be used during the iterations. Several techniques for unconstrained optimization have been studied in recent years; see Ref.[6,11] for an overview of these techniques. The idea is to determine the step length sp by performing the minimization of a univariate function representing a measure of the residual. The function is defined as follows. Let
and
where u is the current solution, p is the step given by the Newton-like algorithm, and s is the line search parameter. We have trivially:
Moreover the first derivative of f at s = 0 is:
where J(u) is the Jacobian of r(u). The matrix-vector product J(u)p is computed via the central finite difference approximation:
with 
= machine
precision. The approximation yields
higher accuracy for the evaluation of the matrix-vector product than a simpler (and cheaper) forward
difference stencil but is twice as
expensive. The objective of the line
search is to find s such
(Notes / (cont'd)
that f (s) = f(u + sp) has decreased sufficiently. It is
always possible to find such a s when p is
a descent direction, i.e., 
fT(u)p
= 
(0) < 0. When
(0) = fT(u)p = rT(u)J(u)p
= rT(u)J(u)(-J-1r)
= -rT r < 0
However, use of Quasi-Newton
updates or Modified Newton Raphson iterations do not always yield a descent
direction. If p is found not to be a descent direction, the solution is stopped
and a restart is created to allow repeat of the load step with a smaller time
increment 
.
Line
Searches with Backtracking: The
objective of the line search in this case is to find s 
] 0, 1] such
that the Goldstein-Armijo condition
f (s) 
f (0) + 
s(0)
is satisfied. Here 
] 0, 1/2] and we
use = 10-4.
Our implementation of the Backtracking algorithm follows essentially the
one presented in [13] (see also [10]).
Line Searches with Line Minimization: The objective of the line search in this case is to find s such that u = u + sp minimizes f(s) in the direction p. The minimization is performed with Brent's method [1]. Our implementation of Brent's method follows essentially the one presented in [13]. Line minimization searches are initiated if the relaxed convergence condition:
is satisfied. The convergence tolerance may be of
the order of 0.5. In our implementation
the value 
= 0.5 has been
adopted.
(6) Only applicable to linear iterations.
(7) If convergence_norm = other, then convergence is checked according to
References / Bibliography
1. Brent, R.P., Algorithms
for Minimization without Derivatives, Chapter 5, Prentice Hall,
2. Brodlie, K.W., Gourlay, A.R., and Greenstadt, J., "Rank-one and Rank-two Corrections to Positive Definite Matrices Expressed in Product Form," J. Inst. Maths Applics, 11, (1973), 73-82.
3. Broyden, C.G. Mathematics of Computation, 19, (1965), 577-593.
4. Dahlquist, G. and A. Bjork, Numerical
Methods, Prentice Hall,
5. Dennis, J.E. and J.J. More, "Quasi-Newton Methods,
Motivation and Theory,"
6. Dennis, J.E. and R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear
Equations, Prentice Hall,
7 Fletcher, R., Practical Methods of Optimization, John Wiley, (2nd edition), (1987).
8. Gill, P.E. and W. Murray, Numerical
Methods for Constrained Optimization, Academic Press,
9. Gill, P.E., W.
10. Johan, Z., T.J.R. Hughes and F. Shakib, "A Globally Convergent Matrix-Free Algorithm for Implicit Time-Marching Schemes Arising in Finite Element Analysis in Fluids," Comp. Meth. Appl. Mech. Eng., 87, (1991), 281-304.
11. Luenberger, D.G., Linear and Nonlinear Programming, Addison-Wesley, (2nd edition), (1984).
12. Matthies, H. and G. Strang, "The Solution of Nonlinear
Finite Element Equations," Int. J.
Num. Meth.
13. Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, (2nd edition), (1992).
14.
15. Woodbury,
M., Inverting Modified Matrices,
Memorandum 42, Statistics Research Group,
Notes . .