Newton-Krylov method
$ \renewcommand{\vec}[1]{ {\bf #1} } \newcommand{\bnabla}{ \vec{\nabla} } \newcommand{\Rey}{Re} \def\vechat#1{ \hat{ \vec{#1} } } \def\mat#1{#1} $
CODE (FORTRAN):
Template/Example: File:Newton Lorenz.f90.
Subroutine only: File:NewtonHook.f90.
CODE (MATLAB / GNU Octave):
Template/Example: File:Newton Lorenz m.zip / File:Newton Lorenz m.tgz
Extended overview pdf (arxiv): [1]
__
The codes above implement the Jacobian-free Newton-Krylov (JFNK) method for solving \begin{equation}\label{eq:ROOTF}\vec{F}(\vec{x})=\vec{0},\end{equation} where $\vec{x}$ and $\vec{F}$ are $n$-vectors, supplemented with a Hookstep--Trust-region approach.
This is a powerful method that can solve for $\vec{x}$ for a complicated nonlinear $\vec{F}(\vec{x})$. For example, to find an equilibrium solution or a periodic orbit, let $\vec{F}(\vec{x}) = \vec{X}(\vec{x})-\vec{x}$, where $\vec{X}(\vec{x})$ is the result of time integration of an initial condition $\vec{x}$.
Newton-Raphson method
To find the roots $x$ of a function $f(x)$ in one dimension, given an initial guess $x_0$, the Newton-Raphson method generates improvements using the iteration $x_{i+1}=x_i-f(x_i)/f'(x_i)$, where the dash denotes the derivative. The iteration can be re-expressed as \[
x_{i+1} = x_i + \delta x_i \,
\quad\mbox{where}\quad
f'(x_i)\,\delta x_i = -f(x_i) \, .
\] The extension of Newton's method to an $n$-dimensional system is then \begin{equation} \label{eq:N-R}
(a)\quad \vec{x}_{i+1} = \vec{x}_i + \vec{\delta x}_i \,
\quad\mbox{where}\quad
(b)\quad\left. \frac{\vec{\partial F}}{\vec{\partial x}}\right|_{\vec{x}_i}
\vec{\delta x}_i = -\vec{F}(\vec{x}_i) \, . \end{equation} In order to apply the update \ref{eq:N-R}(a), the linear system \ref{eq:N-R}(b) needs to be solved for the unknown $\vec{\delta x}_i$.
Jacobian-Free Newton method
The $n\times n$ Jacobian matrix $\left.\frac{\vec{\partial F}}{\vec{\partial x}}\right|_{\vec{x}_i}$ is usually difficult to evaluate. However, the problem \ref{eq:N-R}(b) is in the form $A\,\vec{\delta x}=\vec{b}$, which can be solved using the Krylov-subspace method GMRES(m). (See File:GMRESm.f90.) The GMRES algorithm does not need to know the matrix $A$ itself, only the result of multiplying a vector by $A$. For a given starting vector $\vec{\delta x}_0$, e.g. $\vec{\delta x}_0=\vec{b}$, the method seeks solutions for $\vec{\delta x}$ in $\mathrm{span}\{\vec{x}_0,\,A\vec{x}_0,\,A^2\vec{x}_0,...\}$, but uses Gram-Schmidt orthogonalisation to improve the numerical suitability of this basis set. The set of orthogonalised vectors is called the Krylov-subspace, and m is the maximum number of vectors stored.
Iterations of the GMRES algorithm for the problem \ref{eq:N-R}(b) involve calculating products of the Jacobian with given $\vec{\delta x}$, which may be approximated, e.g. \begin{equation}\label{eq:Japprox}
\left.\frac{\vec{\partial F}}{\vec{\partial x}}\right|_{\vec{x}_i}\,\vec{\delta x} \approx
\frac{1}{\epsilon}(\vec{F}(\vec{x}_i+\epsilon\,\vec{\delta x})-\vec{F}(\vec{x}_i))
\end{equation} for some small value $\epsilon$. (Try $\epsilon$ such that $(\epsilon||\vec{\delta x}||)\,/\,||\vec{x}_i||=10^{-6}$.) The important point is that we do not need to know the Jacobian; only a routine for evaluating $\vec{F}(\vec{x})$ is required.
Note that provided that each step of the Newton method, $\vec{\delta x}$, is essentially in the correct direction, the method is expected to converge. Therefore the tolerance specified in the accuracy of the solution for $\vec{\delta x}$ in each Newton step (calculated via the GMRES method) typically need not be so stringent as the tolerance placed on the Newton method itself for the solution $\vec{x}$. For example, we might seek a relative error for the Newton solution $||\vec{F}(\vec{x})||/||\vec{x}||=O(10^{-8})$, but a relative error for the GMRES steps $||A\,\vec{\delta x}-\vec{b}||/||\vec{\delta x}||=O(10^{-3})$ is likely to be sufficient for calculation of the steps $\vec{\delta x}$.
Hookstep approach
To improve the domain of convergence of the Newton method, it is commonplace to limit the size of the step taken. One approach is simply to take a 'damped' step in the direction of the solution to \ref{eq:N-R}(b), i.e. step by $\alpha\ \vec{\delta x}_i$, where $\alpha \in (0,1]$. In the hookstep approach, we minimise subject to the condition that the magnitude of the Newton step is limited, $||\vec{\delta x}_i||<\delta$, where $\delta$ is the size of the trust region: \begin{equation} \label{eq:dxerr} \min_{\vec{\delta x}_i:\ ||\vec{\delta x}_i||<\delta} \ \left|\left|\left. \frac{\vec{\partial F}}{\vec{\partial x}}\right|_{\vec{x}_i} \vec{\delta x}_i + \vec{F}(\vec{x}_i) \right|\right|\ . \end{equation} Given the minimisation, the hookstep $\vec{\delta x}_i$ is expected to produce a better result than a simple damped step of the same size. It is also expected to perform much better in 'valleys', where it produces a bent/hooked step to a point along the valley, rather than jumping from one side of the valley to the other.
The hookstep can be calculated with little extra work to the GMRES method, provided that the size of Krylov-subspace, m, is chosen sufficiently large to solve to the desired accuracy within m GMRES iterations.
For a given $\vec{\delta x}_i$, the reduction in error predicted by the linearisation \ref{eq:dxerr} can be compared with the actual reduction in $||\vec{F}(\vec{x}_i+\vec{\delta x}_i)||$. According to the accuracy of the prediction, the size of the trust region $\delta$ can be adjusted automatically.
Preconditioning
The GMRES implementations can be supplied with a preconditioner routine (see File:GMRESm.f90), but this can be avoided if the method is combined with time integration (see comment at File:Arnoldi.f).
Further details
Extended overview (arxiv) [2].
On the FORTRAN / MATLAB implementations, see links above.
On GMRES and preconditioning, see File:GMRESm.f90.
On the method, see e.g. Newton-Krylov-Hookstep (channelflow.org).