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}$.