Newton-Krylov method: Difference between revisions

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


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


'''MATLAB / Octave: [[File:NewtonHook.m]].  Template/Example [[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.
___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{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{x}_0$, the method seeks solutions for $\vec{x}$ in $\mathrm{span}\{\vec{x}_0,\,A\vec{x}_0,\,A^2\vec{x}_0,...\}$, but uses Gram-Schmidt orthogonalisation to improve the 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 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 ====
 
To improve the domain of convergence of the Newton method, it is commonplace to limit the size of the step taken.  A simple approach is simply to take a smaller step in the direction of the solution to \ref{eq:N-R}(b).  In the hookstep approach, the error in solving \ref{eq:N-R}(b) is minimised:
\begin{equation} \label{eq:dxerr}
\min_{\vec{\delta x_i}:\ ||\vec{\delta x}_i||<\delta} \  ||\left. \frac{\vec{\partial F}}{\vec{\partial x}}\right|_{\vec{x}_i}
\vec{\delta x}_i + \vec{F}(\vec{x}_i) ||
\end{equation}
i.e. we search over the unknown $\vec{\delta x}_i$ given 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'''.  For a given $\vec{\delta x}_i$, the reduction in error predicted by the linearisation \ref{eq:dxerr} is 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.
 
The hookstep can be calculated with little extra work to the GMRES method, provided that m for the GMRES method is chosen sufficiently large to solve to the desired accuracy within m iterations.
 
==== 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. [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’