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} $

FORTRAN: File:NewtonHook.f90. Template/Example File:Newton Lorenz.f90.

MATLAB / Octave: File:NewtonHook.m. Template/Example File:Newton Lorenz m.tgz.

__

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{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$. Iterations of the GMRES algorithm for the problem \ref{eq:N-R}(b) then 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 of $||\vec{F}(\vec{x})||/||\vec{x}||=O(10^{-8})$, but a relative error for the GMRES steps of $||A\vec{x}-\vec{b}||/||\vec{x}||=O(10^{-3})$ is likely to be sufficient for calculation of the steps $\vec{\delta x}$.

Hookstep approach

In this approach, \ref{eq:N-R}(b) is approximately solved for the unknown $\vec{\delta x}_i$ subject to the condition that the magnitude of the Newton step $\vec{\delta x}$ is smaller than $\delta$, the size of the trust region. The size of $\delta$ is adjusted automatically according to the accuracy of the predicted reduction in the magnitude of $\vec{F}(\vec{x})$ at each Newton step.

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

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