very slow convergence of Picard method for solving nonlinear system of equations I have a nonlinear system of equations as
$$
\left(\mathbf{K}_{\mathbf{L}}+\mathbf{K}_{\mathbf{N L}}(\mathbf{X})\right) \mathbf{X}=\mathbf{F}
$$
in which $\mathbf{K}_{\mathbf{N L}}(\mathbf{X})$ represents the nonlinear stiffness matrix which is dependent to $\mathbf{X}$.
I'm solving with Picard iteration like this:

*

*first ignore the nonlinear stiffness matrix and solve the linear
matrix for $\mathbf{X}$.

*put the resulting $\mathbf{X}$ in the nonlinear stiffness matrix and solve the full equation for $\mathbf{X}$.

*check convergence and repeat 2 if the convergence is not satisfied.

the problem i have here is when the force vector($\mathbf{F}$) is small the nonlinear equation solves very fast but when i increase the force beyond some threshold it gets ages to converge.
i have tried to solve it using Matlab fsolve function with algorithms like 'trust-region' and 'levenberg-marquardt' but the same thing happens with large force vectors.
is there any way i can improve the convergence speed ?
p.s.
heres a gif of the result vector $\mathbf{X}$ inside the convergence loop with a force vector slighly over the threshold.

edit(more details):
so my problem is bending of a nonlinear timoshenko beam that has three governing equations as below:
$$
-\frac{d}{d x}\left\{A_{x x}\left[\frac{d u}{d x}+\frac{1}{2}\left(\frac{d w}{d x}\right)^{2}\right]+B_{x x} \frac{d \phi_{x}}{d x}\right\}=0
$$
$$
-\frac{d}{d x}\left\{A_{x x} \frac{d w}{d x}\left[\frac{d u}{d x}+\frac{1}{2}\left(\frac{d w}{d x}\right)^{2}\right]+B_{x x} \frac{d w}{d x} \frac{d \phi_{x}}{d x}\right\}-\frac{d}{d x}\left[S_{x x}\left(\frac{d w}{d x}+\phi_{x}\right)\right]=q
$$
$$
-\frac{d}{d x}\left\{D_{x x} \frac{d \phi_{x}}{d x}+B_{x x}\left[\frac{d u}{d x}+\frac{1}{2}\left(\frac{d w}{d x}\right)^{2}\right]\right\}+S_{x x}\left(\frac{d w}{d x}+\phi_{x}\right)=0
$$
along with the proper boundary conditions and using finite difference, when assembled they form:
$$
\left(\mathbf{K}_{\mathbf{L}}+\mathbf{K}_{\mathbf{N L}}(\mathbf{X})\right) \mathbf{X}=\mathbf{F}
$$
 A: I assume that $K_{NL}(0) = 0$.
Currently you are using the iteration
$$
(K_L + K_{NL}(X_n))X_{n+1} = F
$$
with $X_0 = 0$. Instead, first solve
$$
(K_L + K_{NL}(X_\sigma))X_\sigma = \sigma F
$$
for small $\sigma$, using this method. This converges quickly as you noticed. Then  solve
$$
(K_L + K_{NL}(X_{\sigma'}))X_{\sigma'} = \sigma' F
$$
with the same iteration for some $\sigma' > \sigma$, using the same iteration but now starting with the previously found $X_\sigma$ and not with the zero vector. And so on until the right hand side is $F$.
For example, choose $\sigma = N^{-1}, \, \sigma' = 2 N^{-1}$ and so on, for sufficiently large $N$.
Of course Anderson acceleration  is also a good idea here :)
A: You can try some of the convergence acceleration algorithms which can work very well for fixed-point (Picard-type) iterations.  If you are using R, there is a package called SQUAREM which implements a reliable, convergence acceleration scheme.  It is based on the paper (Varadhan and Roland, Simple and Globally Convergent Methods for Accelerating the Convergence of Any EM Algorithm, Scandinavian Journal of Statistics, 2008).  EM algorithms are essentially Picard-like algorithms - they are contraction mappings which converge slowly.
