I'm working on modeling two phase immiscible flow in a porous medium. When I setup the system of equations, I obtain a non-linear system of equations that can be expressed in the form:


where the matrix $A(x)$ is a function of the unknown $x$. I know that I can use Newton's method to solve this system, but it is a bit tedious to explicitly write the the derivatives of the left hand side of this equation. Instead, I'd like to use a fixed point iteration of the form:


However, I'm not sure that the iteration will converge to the correct solution. I'm familiar with the contractive mapping theorem as it applies to fixed point iterations in a single dimension. Does this theory extend into higher dimensions as well? Are there any specific conditions on $A(x)$ that guarantee this fixed point iteration will converge to the solution?


The sequence generated by $A(x_n)x_{n+1}=b$ can be written in the form $x_{n+1}=F(x_n):=A(x_n)^{-1}b$, and contracts locally to a solution $x^*$ iff the spectral radius $r$ of the matrix $F'(x^*)=-A(x^*)^{-1}A'(x^*)x^*$ is smaller than one. In this case, $r$ gives the rate of linear convergence. Thus your method will be fast if the $A(x_n)$ are reasonably well-conditioned and depend only weakly on $x_n$.

In all other cases, it is probably far better if you compute a suitable preconditioner $B\approx A(x_0)$ and then solve by one of the standard methods (e.g., Broyden's good method http://en.wikipedia.org/wiki/Broyden%27s_method) the system of equations $F(x)=0$, where $F(x):=B^{-1}(A(x)x-b)$.


The Banach fixed point theorem is extremely general, needing only a metric space. If the map $x \to A(x)^{-1} b$ is globally contractive, the fixed point iteration converges. For your problem (if the formulation is typical of the problem class you describe), it will likely converge linearly, with a basin of attraction that is somewhat larger than basic Newton (without line search, trust region, or continuation methods).

Note that you can still write the system as a defect correction method:

$$\begin{align} \tilde A(x_n) \delta x &= - (A(x_n) x_n - b) = -F(x) \\ x_{n+1} &= x_n + \delta x \end{align}$$

(where $\tilde A = A$ at this point is still your "Picard" linearization). This has identical convergence properties in exact arithmetic, but has a few attractive features in practice:

  1. An explicit residual $F(x) = A(x) x - b$ is available for nonlinear convergence tests, thus you don't have to use the ad-hoc "my solution didn't change much" test which chronically misdiagnoses stagnation for convergence.
  2. The convergence tolerance of the linear solve can usually be greatly relaxed since you are only solving for the defect.
  3. It is also more natural to apply line search or trust region globalization in this context.
  4. It has the structure of a Newton or modified Newton method, thus you can incrementally enrich your approximation $\tilde A$. If you add all terms from Newton linearization, you arrive at precisely a Newton iteration.

As a practical matter, I recommend starting with a Jacobian-free Newton-Krylov method preconditioned by this Picard linearization. If you are using PETSc, just pass your Picard $A$ in both slots of SNESSetJacobian(), then pass -snes_mf_operator. Compare the convergence rate (which should be quadratic) with the (linear) convergence rate you see without -snes_mf_operator.

You should also be aware that there are many nonlinear solution methods, most notably nonlinear GMRES and quasi-Newton, that can accommodate approximate Jacobians such as your Picard linearization. This latter class does not suffer from the finite differencing error of matrix-free finite difference Jacobian application.

  • $\begingroup$ regarding your last remark: The Jacobian $A(x)+A'(x)x$ may be very far from $A(x)$ unless $A(x)$ is almost constant. $\endgroup$ – Arnold Neumaier Sep 28 '12 at 17:18
  • $\begingroup$ Indeed, however, especially for quasi-linear problems, $A(x)$ contains the extremal part of the spectrum and the preconditioner is fairly effective. Note that FAS multigrid usually converges equally fast with a "Picard" correction as with Newton correction. The ability to incrementally add terms is certainly useful unless you have an automatic system to compute all Jacobian terms. $\endgroup$ – Jed Brown Sep 29 '12 at 16:49

Your problem is reminiscent of the self-consistent field (SCF) procedure in electronic structure methods. although these are eigenvalue problems, there are some similarities. Most SCF implementations use the DIIS [1] method to accelerate or enforce convergence. It turns out that DIIS is applicable to nearly any iterative procedure in which convergence is required.

[1] http://en.wikipedia.org/wiki/DIIS

  • $\begingroup$ Are there any iterative procedures in which convergence is not required? $\endgroup$ – chris Sep 28 '12 at 8:01
  • $\begingroup$ Sure, e.g. random number generators require no convergence but are often iterative. $\endgroup$ – Toon Verstraelen Oct 1 '12 at 16:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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