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 \[\vec{F}(\vec{x})=\vec{0},\] 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})$. E.g., for 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}$.

For each Newton iteration we make the improvement $\vec{x}_{i+1} = \vec{x}_i + \vec{\delta x}$. The linear problem for the step $\vec{\delta x}$ is \[\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}=-\vec{F}(\vec{x}_i) \, . \quad\qquad (*) \] In the Hookstep approach, this needs to be approximately solved 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.

The problem $(*)$ is in the form $A\vec{x}=\vec{b}$, where $A$ is an $n\times n$ matrix, and is solved using the Krylov-subspace method, GMRES(m). See File:GMRESm.f90. Iterations of the GMRES(m) algorithm involve calculating products $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}$ for given $\vec{\delta x}$, which may be approximated, by e.g. $\frac{1}{\epsilon}(\vec{F}(\vec{x}_i+\epsilon\,\vec{\delta x})-\vec{F}(\vec{x}_i))$ for some small value $\epsilon$. (Try $\epsilon$ such that $(\epsilon||\vec{\delta x}||)\,/\,||\vec{x}_i||=10^{-6}$.) The important point is that the method does not need to know the matrix $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}$ itself; only a routine for calculating $\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}$.

Preconditioning. The GMRES implementations can be supplied with a preconditioner routine (see File:GMRESm.f90), but this is not necessary if the method is combined with time integration -- strongly decaying modes decay during the time integration (see comment at File:Arnoldi.f).

For further details of the method, see e.g. Newton-Krylov-Hookstep (channelflow.org).