3.9 Newton-Raphson Method

The Newton-Raphson method (also known as “Newton's method”, “Full Newton method” or “the Newton method”) is the basis for solving the nonlinear finite element equations. This section will describe the Full Newton method and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method ¹. The latter variation is actually a quasi-Newton method. It is important since it provides several advantages over the full Newton method and it is this method that is implemented in FEBio ¹.

Full Newton Method

The Newton-Raphson equation (3.1-3) can be written in terms of the discretized equilibrium equations that were derived in the previous section as follows:

\[ \begin{equation} \delta\mathbf{v}^{T}\cdot\mathbf{K}\cdot\mathbf{u}=-\delta\mathbf{v}^{T}\cdot\mathbf{R}.\label{eq355} \end{equation} \]

Since the virtual velocities \(\delta\mathbf{v}\) are arbitrary, a discretized Newton-Raphson scheme can be formulated as follows:

\[ \begin{equation} \begin{array}{cc} \mathbf{K}\left(\mathbf{x}_{k}\right)\cdot\mathbf{u}=-\mathbf{R}\left(\mathbf{x}_{k}\right); & \mathbf{x}_{k+1}=\mathbf{x}_{k}+\mathbf{u}\end{array}.\label{eq356} \end{equation} \]

This is the basis of the Newton-Raphson method. For each iteration \(k\), both the stiffness matrix and the residual vector are re-evaluated and a displacement increment u is calculated by pre-multiplying both sides of the above equation by \(\mathbf{K}^{-1}\). This procedure is repeated until some convergence criteria are satisfied.

The formation of the stiffness matrix and, especially, calculation of its inverse, are computationally expensive. Quasi-Newton methods do not require the reevaluation of the stiffness matrix for every iteration. Instead, a quick update is calculated. One particular method that has been quite successful in the field of computational solid mechanics is the BFGS method, which is described in the next section.

BFGS Method

The BFGS method updates the stiffness matrix (or rather its inverse) to provide an approximation to the exact matrix. A displacement increment is defined as

\[ \begin{equation} \mathbf{d}_{k}=\mathbf{x}_{k}-\mathbf{x}_{k-1}\,,\label{eq357} \end{equation} \]

and an increment in the residual is defined as

\[ \begin{equation} \mathbf{G}_{k}=\mathbf{R}_{k-1}-\mathbf{R}_{k}\,.\label{eq358} \end{equation} \]

The updated matrix \(\mathbf{K}_{k}\) should satisfy the quasi-Newton equation:

\[ \begin{equation} \mathbf{K}_{k}\mathbf{d}_{k}=\mathbf{G}_{k}\,.\label{eq359} \end{equation} \]

In order to calculate this update, as displacement increment is first calculated:

\[ \begin{equation} \mathbf{u}=\mathbf{K}_{k-1}^{-1}\mathbf{R}_{k-1}\,.\label{eq360} \end{equation} \]

This displacement vector defines a “direction” for the actual displacement increment. A line search (see next section) can now be applied to determine the optimal displacement increment:

\[ \begin{equation} \mathbf{x}_{k}=\mathbf{x}_{k-1}+s\mathbf{u}\,,\label{eq361} \end{equation} \]

where \(s\) is determined from the line search. With the updated position calculated, \(\mathbf{R}_{k}\) can be evaluated. Also, using equations \eqref{eq357} and \eqref{eq358}, \(\mathbf{d}_{k}\) and \(\mathbf{G}_{k}\) can be evaluted. The stiffness update can now be expressed as

\[ \begin{equation} \mathbf{K}_{k}^{-1}=\mathbf{A}_{k}^{T}\mathbf{K}_{k-1}^{-1}\mathbf{A}_{k}\,,\label{eq362} \end{equation} \]

where the matrix \(\mathbf{A}\) is an \(n\times n\) matrix of the simple form:

\[ \begin{equation} \mathbf{A}_{k}=\mathbf{1}+\mathbf{v}_{k}\mathbf{w}_{k}^{T}\,.\label{eq363} \end{equation} \]

The vectors \(\mathbf{v}\) and \(\mathbf{w}\) are given by

\[ \begin{equation} \mathbf{v}_{k}=-\left(\frac{\mathbf{d}_{k}^{T}\mathbf{G}_{k}}{\mathbf{d}_{k}^{T}\mathbf{K}_{k-1}\mathbf{d}_{k}}\right)^{1/2}\mathbf{K}_{k-1}\mathbf{d}_{k}-\mathbf{G}_{k}\,,\label{eq364} \end{equation} \]

\[ \begin{equation} \mathbf{w}_{k}=\frac{\mathbf{d}_{k}}{\mathbf{d}_{k}^{T}\mathbf{G}_{k}}\,.\label{eq365} \end{equation} \]

The vector \(\mathbf{K}_{k-1}\mathbf{d}_{k}\) is equal to \(s\mathbf{R}_{k-1}\) and has already been calculated.

To avoid numerically dangerous updates, the condition number \(c\) of the updating matrix \(\mathbf{A}\) is calculated:

\[ \begin{equation} c=\left(\frac{\mathbf{d}_{k}^{T}\mathbf{G}_{k}}{\mathbf{d}_{k}^{T}\mathbf{K}_{k-1}\mathbf{d}_{k}}\right)^{1/2}\,.\label{eq365b} \end{equation} \]

The update is not performed when this number exceeds a preset tolerance.

Considering the actual computations involved, it should be noted that using the matrix updates defined above, the calculation of the search direction in \eqref{eq360} can be rewritten as,

\[ \begin{equation} \mathbf{u}=\left(\mathbf{1}+\mathbf{w}_{k-1}\mathbf{v}_{k-1}^{T}\right)\cdots\left(\mathbf{1}+\mathbf{w}_{1}\mathbf{v}_{1}^{T}\right)\mathbf{K}_{0}^{-1}\left(\mathbf{1}+\mathbf{v}_{1}\mathbf{w}_{1}^{T}\right)\cdots\left(\mathbf{1}+\mathbf{v}_{k-1}\mathbf{w}_{k-1}^{T}\right)\mathbf{R}_{k-1}\,.\label{eq366} \end{equation} \]

Hence, the search direction can be computed without explicitly calculating the updated matrices or performing any additional costly matrix factorizations as required in the full Newton-Raphson method.

Line Search Method

A powerful technique often used to improve the convergence rate of Newton based methods is the line search method. In this method, the direction of the displacement vector \(\mathbf{u}\) is considered as optimal, but the magnitude is controlled by a parameter \(s\):

\[ \begin{equation} \mathbf{x}_{k+1}=\mathbf{x}_{k}+s\mathbf{u}\,.\label{eq367} \end{equation} \]

The value of \(s\) is usually chosen so that the total potential energy \(W\left(s\right)=W\left(\mathbf{x}_{k}+s\mathbf{u}\right)\) at the end of the iteration is minimized in the direction of \(\mathbf{u}\). This is equivalent to the requirement that the residual force \(\mathbf{R}\left(\mathbf{x}_{k}+s\mathbf{u}\right)\) at the end of the iteration is orthogonal to \(\mathbf{u}\):

\[ \begin{equation} R\left(s\right)=\mathbf{u}^{T}\mathbf{R}\left(\mathbf{x}_{k}+s\mathbf{u}\right)=0\,.\label{eq368} \end{equation} \]

However, in practice it is sufficient to obtain a value of \(s\) such that,

\[ \begin{equation} \left|R\left(s\right)\right|<\rho\left|R\left(0\right)\right|\,,\label{eq369} \end{equation} \]

where typically a value of \(\rho=0.9\) is used. Under normal conditions the value \(s=1\) automatically satisfies equation \eqref{eq369} and therefore few extra operations are involved. However, when this is not the case, a more suitable value for \(s\) needs to be obtained. For this reason it is convenient to approximate \(R\left(s\right)\) as a quadratic in \(s\):

\[ \begin{equation} R\left(s\right)\approx\left(1-s\right)R\left(0\right)+R\left(1\right)s^{2}=0\,,\label{eq370} \end{equation} \]

which yields a value for \(s\) as

\[ \begin{equation} s=\frac{r}{2}\pm\sqrt{\left(\frac{r}{2}\right)^{2}-r}\,,\quad r=\frac{R\left(0\right)}{R\left(1\right)}\,.\label{eq371} \end{equation} \]

If \(r<0\), the square root is positive and a first improved value for \(s\) is obtained:

\[ \begin{equation} s_{1}=\frac{r}{2}+\sqrt{\left(\frac{r}{2}\right)^{2}-r}\,.\label{eq372} \end{equation} \]

If \(r>0\) the \(s\) can be obtained by using the value that minimizes the quadratic function, that is, \(s_{1}=r/2\). This procedure is now repeated with \(R\left(1\right)\) replaced by \(R\left(s_{1}\right)\) until equation \eqref{eq369} is satisfied.

Conjugate gradient solution method

The Newton-Raphson solution algorithms are highly effective for nearly linear problems, but they have two major limitations. Firstly, for very large models, the size of the stiffness matrix K can grow very large as well. The factorization of this matrix, which is needed when using direct linear solvers, can become very time consuming and requires a large amount of memory. This is the main factor that limits the size of model that can be solved on a given computer. Secondly, where the stiffness is small, it can overestimate the displacement, resulting in excessive distortion of the mesh and the familiar negative Jacobian failure. This can occur when a material has a low initial stiffness that increases at higher strains (we work around this by using a large number of small timesteps to reduce the nonlinearity on each increment), when buckling causes a temporary low or negative stiffness, when a part of the model is not sufficiently constrained, or when contacts cause sudden changes of stiffness. For these reasons the Newton-Raphson algorithms often fail to converge for very nonlinear materials, models with contacts or unstable structures that can buckle or crumple under load.

An alternative strategy is to estimate a direction to move each node directly from the residual R:

\[ \boldsymbol{u=-}\boldsymbol{R\left(x_{k}\right)} \]

Then a line search is needed to find the appropriate distance to move in that direction. This is the steepest descent method, which has the disadvantage that it converges slowly in a series of decreasing zig-zag steps. Its convergence is much improved by adding a fraction β of the previous step:

\[ \boldsymbol{u}_{k}=-\boldsymbol{R}\left(\boldsymbol{x}_{k}\right)+\beta\boldsymbol{u}_{k-1} \]

There are many different ways to calculate β; one method which has excellent convergence for nonlinear problems is due to Hager and Zhang:

\[ \boldsymbol{y}_{i}=\boldsymbol{R}_{i-1}-\boldsymbol{R}_{i} \]

\[ \beta_{i}=\frac{-1}{\boldsymbol{u}_{i-1}\cdot\boldsymbol{y}_{i}}\left(\boldsymbol{y}_{i}-2\boldsymbol{u}_{i-1}\frac{\left\Vert \boldsymbol{y}_{i}\right\Vert ^{2}}{\boldsymbol{u}_{i-1}\cdot\boldsymbol{y}_{i}}\right)\cdot\boldsymbol{R}_{i} \]

\[ \eta_{i}=\frac{-1}{\left\Vert \boldsymbol{u}_{i-1}\right\Vert min\left\{ \eta,\left\Vert \boldsymbol{R}_{i-1}\right\Vert \right\} } \]

\[ \beta_{i}=max\left\{ \beta_{i},\eta_{i}\right\} \]

A limitation of this method is that the distance each node moves is proportional to the residual; stiffer nodes produce a larger residual for a given displacement and therefore move further, whereas in fact they need to move less than more compliant nodes. This results in mesh distortion where the element size varies, for midside nodes in quadratic elements or for mixed materials. This problem can be reduced by using a preconditioner, dividing the residual at each node by an estimate of its stiffness. This requires changes also to the calculation of β:

\[ \beta_{i}=\frac{-1}{\boldsymbol{u}_{i-1}\cdot\boldsymbol{y}_{i}}\left(\boldsymbol{y_{i}^{T}P}-2\boldsymbol{u}_{i-1}\frac{\boldsymbol{y_{i}^{T}Py}_{i}}{\boldsymbol{u}_{i-1}\cdot\boldsymbol{y}_{i}}\right)\cdot\boldsymbol{R}_{i} \]

\[ \eta_{i}=0.4\frac{\Delta\boldsymbol{u}_{i-1}^{T}\boldsymbol{R}_{i}}{\boldsymbol{u_{i-1}^{T}P^{-1}u_{i-1}}} \]

\[ \beta_{i}=max\left\{ \beta_{i},\eta_{i}\right\} \]

\[ \boldsymbol{u}_{i+1}=\boldsymbol{P}\boldsymbol{R}_{i}+\beta_{i}\boldsymbol{u}_{i} \]

The preconditioner P can be found in various ways; a convenient approximation is to use the reciprocal of the main diagonal of the stiffness matrix:

\[ \boldsymbol{P}^{-1}=diag\left(\boldsymbol{K}\right) \]

This is calculated once at the beginning of the solution.

---

1. Matthies, H.; Strang, G.. "The solution of nonlinear finite element equations." Intl J Num Meth Eng, vol. 14, pp. 1613-26 (1979). ↩↩