# Newton-Raphson method fails!

I am trying to solve an equation like $R(x) = 0$, using Newton-Raphson method. To obtain the $x$ increment in each iteration I solve $dx = -(A)^{-1}\cdot R$ where $A = dR/dx$. But the convergence criteria, here $||R||$ reduces to a certain amount, depending on the input, e.g. 1.e-8, and it remains so. After much trying I did not find any special bug, so i am starting to have second thoughts about my implementation.

I am wondering now that the approach to calculate the inverse of matrix $A$ might cause inaccuracy so that the convergence fails. What do you think?

EDIT: The code is written in fortran; I calculate the inverse using MKL's routines: dgetrf and dgetri. I am using double-precision variables.

The problem is to solve $R(\sigma)$ in n steps using Backward-Euler method. So functions $Z^p$ and $Z^d$ are also dependent on $\sigma$. In each step $\Delta\epsilon_{n+1}$ is updated and we solve an implicit problem to update $\sigma$

where for a plane strain problem the stress tensor reads

$${\boldsymbol\sigma} = \left[ \begin{array}{c} \sigma_{11}\\ \sigma_{22}\\ \sigma_{33}\\ \sigma_{12} \end{array} \right]$$

The function $R$ is: $$\mathbf{R}_{\boldsymbol\sigma} = {\boldsymbol\sigma}_{n+1} - {\boldsymbol\sigma}_n - \mathbf{D}\Delta{\boldsymbol\epsilon_{n}} + \Delta t\,\mathbf{D}\left[\mathbf{Z}^p({\boldsymbol\sigma}^{\rm{dev}}) + \mathbf{Z}^d({\boldsymbol\sigma}^{\rm{dev}})\right]$$ where \begin{align} \left[\mathbf{Z}^d({\boldsymbol\sigma}^{\rm{dev}})\right]_i &= \alpha\,\left[{\boldsymbol\sigma}^{\rm{dev}}\right]_i \\ \left[\mathbf{Z}^p({\boldsymbol\sigma}^{\rm{dev}})\right]_i &= \beta\,||{\boldsymbol\sigma}^{\rm{dev}}||^{\frac{1-m}{m}} \left[{\boldsymbol\sigma}^{\rm{dev}}\right]_i \\ \left[{\boldsymbol\sigma}^{\rm{dev}}\right]_i &= \left[\mathbf{I}^{\rm{dev}}\right]_{ij}\,\left[{\boldsymbol\sigma}\right]_j \end{align}

$I^{dev}$ is defined as:

$$I_{ij}^{dev} = \delta_{ij} - \frac{1}{3}m_im_j \quad\text{where}\quad \mathbf{m} = \left[ \begin{array}{c} 1\\ 1\\ 1\\ 0 \end{array} \right]$$

and $\delta$ is the Kronecker delta.

The tangent modulus is given by $$d\mathbf{R}_{\boldsymbol{\sigma}} = \frac{\partial \mathbf{R}_{\boldsymbol{\sigma}}}{\partial {\boldsymbol{\sigma}}}\,d{\boldsymbol{\sigma}} = \left[\mathbf{I} + \Delta t\,\mathbf{D}\left(\frac{\partial \mathbf{Z}^d({\boldsymbol\sigma}^{\rm{dev}})}{\partial \boldsymbol{\sigma}} + \frac{\partial \mathbf{Z}^p({\boldsymbol\sigma}^{\rm{dev}})}{\partial \boldsymbol{\sigma}}\right) \right]\,d{\boldsymbol{\sigma}} = \mathbf{A}\,d{\boldsymbol{\sigma}}$$ If we expand $\mathbf{A}$ we get $$\mathbf{A} = \mathbf{I} + \Delta t\,\mathbf{D}\left[\left(\alpha+\beta\,||{\boldsymbol\sigma}^{\rm{dev}}||^{\frac{1-m}{m}}\right)\,\mathbf{I}^{\rm{dev}} + \frac{(1-m)\beta}{m}\,||{\boldsymbol\sigma}^{\rm{dev}}||^{\frac{1-3m}{m}}{\boldsymbol\sigma}^{\rm{dev}}\otimes{\boldsymbol\sigma}^{\rm{dev}}\right]$$ Here $x=\sigma_{n+1}$. Also, $\alpha, \beta, \Delta t, I_{ij}^{dev}, D$ and $m$ are constant($m=0.3$ and $i,j=1,...,4$.). $\sigma_n$ and $\Delta{\boldsymbol\epsilon_{n}}$ are constant in step n.

Here you can see the results for $R$ and $\sigma_{n+1}$.

• What is your expression for $A$? Newton iterations often don't converge (or converge slowly) if the Jacobian has been derived from continuum mechanics considerations but is not consistent with the numerical integration method. Apr 11, 2014 at 3:47
• I am not sure I get your point. Can you explain more? (I added $A$ to the post.) Apr 11, 2014 at 7:56
• The primary reason for Newton iterations failing to converge, in my experience, is that the Jacobian is wrong. If the Jacobian is correct, then the cause is an instability in which case a line search or arc-length method can fix the problem. In your case, you haven't explained what your quantities are - matrices, vectors, what size etc. You have also not defined $\sigma^{dev}$ and $I^{dev}$. So it's hard to say what's wrong. For example, if $\sigma$ is the stress tensor you have to consider all nine components when you take derivatives. Have you done that? Apr 13, 2014 at 2:05
• stress and strain are rank-1 tensors with dimension of 4 for the case of plane strain. I don't believe that is an issue. Maybe the problem is the instability. Apr 13, 2014 at 7:51
• I always like to think of stress and strain as rank-2 tensors for physical reasons - independently of whether a 2D assumption is used. You model is relatively simple and should be easy to test. One ways of finding out the source of the problem is to test your model for a 1D state of stress (uniaxial stress) to find out whether there actually is an instability. Also, you can look at what characteristic the acoustic tensor has. The initial rate of convergence seems to be quadratic. Do you have a convergence plot? What are the eigenvalues of A? Apr 13, 2014 at 21:45

As a check we can look at the problem in terms tensors. Let $$\boldsymbol{\sigma}^{\text{dev}} = \boldsymbol{\sigma} - \frac{1}{3}(\boldsymbol{\sigma}:\mathbf{I})\,\mathbf{I}$$ where $\mathbf{I}$ is the second-order identity tensor. Then the fourth-order tensor $\mathsf{I}^{\rm dev}$ is given by $$\mathsf{I}^{\rm dev} = \mathsf{I}^{(4s)} - \frac{1}{3} \mathbf{I}\otimes\mathbf{I}$$ where $\mathsf{I}^{(4s)}$ is the symmetric fourth-order identity tensor with components $1/2(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})$ and $\delta_{ij}$ is the Kronecker delta.
We can also express $\mathbf{Z}^p$ in tensor form as $$\boldsymbol{Z}^p = \beta\,\left\|\boldsymbol{\sigma}^{\rm dev}\right\|^{\frac{1-m}{m}} \boldsymbol{\sigma}^{\rm dev} \,.$$ Then $$\frac{\partial \boldsymbol{Z}^p}{\partial\boldsymbol{\sigma}} = \beta\left[\frac{\partial}{\partial \boldsymbol{\sigma}}\left(\left\|\boldsymbol{\sigma}^{\rm dev}\right\|^{\frac{1-m}{m}}\right)\otimes\boldsymbol{\sigma}^{\rm dev} + \left\|\boldsymbol{\sigma}^{\rm dev}\right\|^{\frac{1-m}{m}} \frac{\partial \boldsymbol{\sigma}^{\rm dev}}{\partial\boldsymbol{\sigma}}\right]$$ Let us look at the first term $$\frac{\partial}{\partial \boldsymbol{\sigma}}\left(\left\|\boldsymbol{\sigma}^{\rm dev}\right\|^{\frac{1-m}{m}}\right) = \frac{1-m}{m}\left\|\boldsymbol{\sigma}^{\rm dev}\right\|^{\frac{1-2m}{m}} \frac{\partial\left\|\boldsymbol{\sigma}^{\rm dev}\right\|}{\partial \boldsymbol{\sigma}}$$ Next look at the part $$\frac{\partial\left\|\boldsymbol{\sigma}^{\rm dev}\right\|}{\partial \boldsymbol{\sigma}} = \frac{1}{2\left\|\boldsymbol{\sigma}^{\rm dev}\right\|}\,\frac{\partial}{\partial\boldsymbol{\sigma}}(\boldsymbol{\sigma}^{\rm dev}:\boldsymbol{\sigma}^{\rm dev})$$ Finally, we look at the last bit above, using $\boldsymbol{\sigma}^{\rm dev}:\mathbf{I} = 0$, $$\frac{\partial}{\partial\boldsymbol{\sigma}}(\boldsymbol{\sigma}^{\rm dev}:\boldsymbol{\sigma}^{\rm dev}) = \mathsf{I}^{\rm dev}:\boldsymbol{\sigma}^{\rm dev} + \boldsymbol{\sigma}^{\rm dev}:\mathsf{I}^{\rm dev} = 2\boldsymbol{\sigma}^{\rm dev}$$ Therefore, $$\frac{\partial \boldsymbol{Z}^p}{\partial\boldsymbol{\sigma}} = \frac{\beta(1-m)}{m}\left\|\boldsymbol{\sigma}^{\rm dev}\right\|^{\frac{1-3m}{m}} \boldsymbol{\sigma}^{\rm dev}\otimes\boldsymbol{\sigma}^{\rm dev} + \beta\,\left\|\boldsymbol{\sigma}^{\rm dev}\right\|^{\frac{1-m}{m}}\mathsf{I}^{\rm dev} \,.$$ This is different from the expression you have used. So you will need to recheck your algebra. So the tangent modulus appears to be correct.