I am trying to solve two non-linear equations self-consistently in a Gummel loop. Sometimes (every once in a while), I get to a situation when the loop repeats itself with wrong solutions and a certain error persists. As a simple example, consider the following two equations:

\begin{align} &y = -x + 1 \\ &y = \sqrt{x} \end{align} and suppose the loop reaches to $x=0$ for the input of first equation, which leads to $y=1$ for the input of the second equation. This results in $x=1$ for the input of the first equation, leading to $y=0$, and the situation repeats itself (and of course, does not converge to a correct solution).

I was wondering if there is a good and comprehensive reference on this particular problem and on the properties of equations which lead to such behaviour. Also, what is the best way to avoid such difficulties in general?

  • 1
    $\begingroup$ If you in mind a setting in which two sets of unknowns are not only coupled by a nonlinear system but are unknowns drawn from high-dimensional spaces, it would improve your Question to have a more faithful illustration. Identify those aspects of the problem which are considered advantageous to numerical solutions, so that Readers have a chance to reply with methods to exploit those advantages. $\endgroup$
    – hardmath
    Sep 9, 2015 at 22:12

1 Answer 1


In general, the iterative method is not required to converge. This is commonly discussed in numerical analysis textbooks in relation to nonlinear root-finding (e.g., here it is in Analysis of Numerical Methods by Isaacson and Keller, p.85--).

You are iterating the function $$ f(x) = (x-1)^2 $$ and solving the equation $x = f(x)$ with the roots $x_{1,2} = \frac12(3\pm \sqrt{5}) = 0.38, 2.62$. (This form of the equation is easier to analyse than your version with both variables $x$ and $y$.)

For the iterative method to succeed, according to a standard theorem, it is sufficient for the initial guess to lie in (or to eventually converge to a value in) an interval $(x_-,x_+)$ surrounding a root such that $|f'(x)|\leq \lambda < 1$ for each $x\in(x_-,x_+)$. Since $f'(x) = 2(x-1)$, we can see that near the roots $|f'(x)|>1$ ($f'(x_{1,2}) = -1.24, 3.24$). This means that each root is unstable: for each root the interval surrounding it from which initial guesses are guaranteed to converge is empty. This unstable behaviour is much easier to see if you start from an initial guess that is already near a root, $x_{1,2}+\epsilon$, instead of $0$.

This issue can be solved, for example, by rewriting the equations in an equivalent form that happens to lead to stable roots.

  • $\begingroup$ By constructing this simple system of equations, I just wanted to present a simplified explanation of my problem. The actual system that I am trying to solve consists of non-linear Poisson equation and the Schrodinger equation, which cannot be rewritten in an equivalent more stable form. Therefore, a more general discussion or reference would be of great help to me. $\endgroup$ Sep 9, 2015 at 13:57
  • 2
    $\begingroup$ @MaziarNoei I feel like you gotta put such things explicitly into the question, it's a bit misleading otherwise - there's no way for me to guess what you actually meant. The same theory applies without too many modifications to the general case: $f$ has to be Lipschitz near the root with Lipschitz constant $<1$ ($f$ has to be a contraction, so that a fixed point theorem applies). $\endgroup$
    – Kirill
    Sep 9, 2015 at 15:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.