Convergence issues for a non-linear system

Question

I have a nasty system of coupled integral equations, which I managed to discretize and recast a non-linear system, i.e. something like:

$$ \vec{w} = F \left( \vec{w} \right) \hspace{32pt} w \in \mathbb{R}^n$$

I am employing a naive iterative solution method. For some reasonable initial choice of $\vec{w}$ I repeatedly calculate $ \vec{w}^{(n+1)} = F \left( \vec{w}^{(n)} \right) $. I define at each step the quantity $x= || \vec{w}^{(n+1)} - \vec{w}^{n} ||$.

I observe that $x$ is decreasing (as expected!) iteration after iteration until I am close the the solution of the system (I know the solution for some particular values of the parameters...) Then $x$ starts rapidly increasing, and it seems like my method diverges progressively from the real solution of the system.

How can I solve my problem? My guess would be taking the solution which produced the lowest $x$, and maybe trying some small perturbation to achieve better convergence starting from there, but that's just a guess. Just taking the solution with the lowest $x$ won't work, because I would like to have better convergence.

Thank you in advance!

Answer 1

A common strategy is to employ a damping strategy, i.e., to compute $\vec{w}^{\ast} = F \left( \vec{w}^{(n)} \right)$ and then set $\vec{w}^{(n+1)} = \alpha \vec{w}^{(\ast)} + (1-\alpha)\vec{w}^{(n)}$ where $\alpha\in[0,1]$. You typically choose $\alpha$ in such a way that it minimizes $\|[\alpha \vec{w}^{(\ast)} + (1-\alpha)\vec{w}^{(n)}] - F(\alpha \vec{w}^{(\ast)} + (1-\alpha)\vec{w}^{(n)})\|$ in some sense. Damping strategies like this are often more robust than just taking a full step $\alpha=1$, which is what you are currently doing.

Answer 2

Look into Anderson mixing. It is used for self consistent quantum chemistry applications which are nonlinear. It is faster and more robust than simple damping.

The headline is that you keep a history of n previous iterations and calculate the residual for each, then perform a linear solve on the set of residuals using the approximation that the solution space is linear near the solution. This works very well, and far faster than the similar Broyden update method, which Eyert proved was equivalent to Anderson.

It fell down slightly when I used it, because my transform was so sensitive to small deltas even near the solution, i.e. the linear approximation did hold hold too well.