Newton-Krylov method: Difference between revisions

From openpipeflow.org
Jump to navigation Jump to search
mNo edit summary
mNo edit summary
 
(11 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{latexPreamble}}
=== Jacobian-Free Newton-Krylov (JFNK) method ===
Free to download:


'''CODE (FORTRAN): Template-Example: [[File:Newton_Lorenz.f90]].'''
'''This page and the codes have been moved to https://github.com/apwillis1/JFNK-Hookstep '''


'''CODE (MATLAB / GNU Octave): Template-Example: [[File:Newton_Lorenz_m.zip]] / [[File:Newton_Lorenz_m.tgz]]'''
Codes are in MATLAB (Octave compatible) and Fortran90.   
 
They can be integrated with codes in other languages via command-line calls.
'''Info on main subroutine:''' (essentially same for FORTRAN/MATLAB versions)  [[File:NewtonHook.f90]]
 
___ __
 
'''(An Extended overview of the Newton-Krylov method is here (pdf,arxiv): [https://arxiv.org/pdf/1908.06730.pdf]  See section 4 on using the code.)'''
 
___NOTOC___ __
 
 
 
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) [https://arxiv.org/pdf/1908.06730.pdf].
 
Usage of NewtonHook subroutine [[File:NewtonHook.f90]] (essentially same for FORTRAN and MATLAB implementations).
 
On GMRES and preconditioning, see [[File:GMRESm.f90]].
 
More on the Newton-Krylov method, see e.g. [http://channelflow.org/dokuwiki/doku.php?id=docs:math:newton_krylov_hookstep Newton-Krylov-Hookstep] (channelflow.org).

Latest revision as of 08:50, 1 February 2023

Jacobian-Free Newton-Krylov (JFNK) method

This page and the codes have been moved to https://github.com/apwillis1/JFNK-Hookstep

Codes are in MATLAB (Octave compatible) and Fortran90. They can be integrated with codes in other languages via command-line calls.

Retrieved from ‘http:///index.php?title=Newton-Krylov_method&oldid=1059’