7
$\begingroup$

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:

$A(x)x=b$

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:

$A(x_{n})x_{n+1}=b$

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?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
5
$\begingroup$

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

Improve this answer
$\endgroup$
8
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\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, 2012 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, 2012 at 16:49
1
$\begingroup$

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

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

Your Answer

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

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