# File:NewtonHook.f90

$ \renewcommand{\vec}[1]{ {\bf #1} } \newcommand{\bnabla}{\vec{\nabla}} \newcommand{\Rey}{Re} \def\vechat#1{\hat{\vec{#1}}} \def\mat#1{#1} $

**Download via link above. See also** File:Newton template.f90.

## The Newton-Krylov method

The code above implements the Jacobian-free Newton-Krylov (JFNK) method for solving \[\vec{F}(\vec{x})=\vec{0},\] where $\vec{x}$ and $\vec{F}$ are $n$-vectors, using 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}$. (See the comment at File:Arnoldi.f on how timestepping or exponentiation can help suitability for Krylov methods.)

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 $\vec{\delta x}$ is smaller than $\delta$, the size of the trust region. The size of $\delta$ is automatically adjusted 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$. Thus, 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}$.

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

## The code

To download, click the links at top of page. A template for calling the routine: File:Newton template.f90. The code also requires File:GMRESm.f90 and Lapack.

- Prior to calling the subroutine
`newtonhook(...)`, the arrays`new_x(:)`$\equiv\vec{x}_i$ and`new_fx(:)`$\equiv\vec{F}(\vec{x}_i)$ need to be allocated the dimension $n$. - Also prior to calling
`newtonhook(...)`,`new_x(:)`must be assigned the initial guess $\vec{x}_0$. (As the calculation proceeds, the vector`new_x(:)`is updated with improved vectors, having smaller corresponding`new_fx(:)`.)

- In addition to scalar and array variables,
`newtonhook(...)`needs to be passed- a function that calculates the dot product of two vectors,
- a function that calculates $\vec{F}(\vec{x})$ for a given $\vec{\delta x}$,
- a function that calculates $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}$ for a given $\vec{\delta x}$. See possible approximation above, where $\vec{F}(\vec{x}_i)$ is already stored in
`new_fx(:)`. - a subroutine that replaces a vector $\vec{x}$ with the solution of $M\vec{x}'=\vec{x}$ for $\vec{x}'$, where $M$ is a preconditioner matrix. This may simply be an empty subroutine if no preconditioner is required, i.e. $M=I$.
- a subroutine that is called at the end of each iteration. This may be used to save the current state after each iteration, if desired.

- The functions above may require auxiliary data, in addition to the given $\vec{x}$ or $\vec{\delta x}$. A few ways to access the data:
- place the required data in a module and access via '
`use mymodule`'. - place the data in a
`common`block. - write the vector to disk, call a script that runs an external programme, then load the result.

- place the required data in a module and access via '

- Extra constraints may be necessary, e.g. that determine an update to the period of an orbit. The function that evaluates $\left.\frac{\vec{dF}}{\vec{dx}}\right|_{\vec{x}_i}\,\vec{\delta x}$ for a given $\vec{\delta x}$ should append to the result the evaluations of the constraints. Correspondingly, for each constraint an extra
`0`should be appended to $\vec{F}(\vec{x}_i)$ (the evaluated constraint should equal zero when converged).

## Parallel use

It is NOT necessary to edit this code for parallel (MPI) use:

- let each thread pass its subsection of the vector $\vec{x}$,
- make the dot product function
`mpi_allreduce`the result of the dot product. - to avoid multiple outputs to the terminal, set
`info=1`on rank 0 and`info=0`for the other ranks.

## File history

Click on a date/time to view the file as it appeared at that time.

Date/Time | Dimensions | User | Comment | |
---|---|---|---|---|

current | 04:32, 13 December 2016 | (6 KB) | Apwillis (Talk | contribs) | == The method == Implements the Newton-Raphson method for solving \[\vec{F}(\vec{x})=\vec{0}\] with a Hookstep-Trust-region approach. At each iteration, the linear problem for the step in $\vec{x}$ is solved subject to the condition that the step be ... |

- You cannot overwrite this file.

## File usage

The following 4 pages link to this file: