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}$.
Newton-Raphson method fails!
fortran;newton method