Newton-Krylov method: Difference between revisions

From openpipeflow.org
Jump to navigation Jump to search
mNo edit summary
mNo edit summary
 
(42 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})$.  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}$.
 
==== 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(m) algorithm for the problem \ref{eq:N-R}(b) then 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, e.g. by
\begin{equation}\label{eq:Japprox} \left.\frac{\vec{dF}}{\vec{dx}}\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 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’