Krylov-Newton algorithm

• Krylov-Newton algorithm
• • Convergence criterion
  • Options of the Newton algorithm
  • Space of solutions
  • Different Jacobians
  • Simple example
  • Other flavours of newton
  • Example

BifurcationKit is built upon the newton algorithm for solving (large-dimensional) nonlinear equations

\[F(x)=0\in\mathbb R^n,\quad x\in\mathbb R^n.\]

Writing $J(x)\in\mathcal L(\mathbb R^n)$ the jacobian, the algorithm reads

\[x_{n+1} = x_n - J(x_n)^{-1}F(x_n)\]

with initial guess $x_0$.

The crux of the algorithm is to solve the linear system in $y$:

\[J(x_n)\cdot y = F(x_n).\]

To this end, we never form $J^{-1}$ like with pinv(J) but solve the linear system directly.

Convergence criterion

The Newton iteration is stopped when the (scaled) residual becomes smaller than a user tolerance. More precisely, at iteration $k$ we compute $r_k:=N\big(F(x_k)\big)$ where $N$ is a norm that can be specified by the keyword normN of BifurcationKit.solve (by default normN = norm, the $\ell^2$ norm). The iteration is declared converged as soon as

\[r_k < \mathtt{tol} \quad\text{and}\quad \text{callback}(x_k, p, F(x_k)) \text{ returns true},\]

where tol is the field NewtonPar.tol (default $10^{-12}$) and callback is a keyword of solve which allows to reject (or stop the iteration on) steps based on extra information, e.g. the norm of the residual or the norm of the Newton step. BifurcationKit provides the helpers BifurcationKit.cbMaxNorm and BifurcationKit.cbMaxNormAndΔp to build such callbacks, and they can be passed to solve/continuation (as callback/callback_newton respectively). We refer to the output of continuation iterations for examples of use.

Options of the Newton algorithm

The Newton iterations are controlled by the composite type NewtonPar. The most useful fields are:

Space of solutions

For the algorithm to be defined, a certain number of operations on x need to be available. If you pass x::AbstractArray, you should not have any problem. Otherwise, your x must comply with the requirements listed in Required methods for custom arrays.

Different Jacobians

There are basically two ways to specify the jacobian:

1. Matrix based
2. Matrix-free.

In case you pass a matrix (in effect an AbstractMatrix like a sparse one,...), you can use the default linear solver from LinearAlgebra termed the backslash operator \. This is a direct method. This is the case 1 above.

Another possibility is to pass a function J(dx) and to use iterative linear solvers. In this case, this is termed a Krylov-Newton method. This is the case 2 above. In comparison to the Matrix-based case, there is no restriction to the number of unknowns $n$.

The available linear solvers are explained in the section Linear solvers (LS).

One can find a full description of the Krylov-Newton method in the docstring of BifurcationKit.solve.

Simple example

Here is a quick example to show how the basics work. In particular, the problem generates a matrix based jacobian using automatic differentiation.

using BifurcationKitF(x, p) = x.^3 .- 1x0 = rand(10)prob = ODEBifProblem(F, x0, nothing)sol = BifurcationKit.solve(prob, Newton(), NewtonPar(verbose = true))

NonLinearSolution{Vector{Float64}, ODEBifProblem{BifFunction{typeof(Main.F), BifurcationKit.var"#169#171"{typeof(Main.F)}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Float64, BifurcationKit.Jet{Float64, Nothing, BifurcationKit.AutoDiff, BifurcationKit.FiniteDifferences, Nothing, Nothing, BifurcationKit.AutoDiff, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}}, Vector{Float64}, Nothing, typeof(identity), typeof(BifurcationKit.plot_default), typeof(BifurcationKit.record_sol_default), typeof(BifurcationKit.save_solution_default), typeof(BifurcationKit.update_default)}, Vector{Float64}, Int64}([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], ┌─ Bifurcation problem with uType Vector{Float64}
├─ Inplace: false
├─ Dimension: 10
├─ Symmetric: false
└─ Parameter: p, [2.4486611753738186, 2727.25227619637, 807.6720295152037, 238.90879391155843, 70.38995838757765, 20.469732538821, 5.711099442422088, 1.4102483208678482, 0.24530852292277322, 0.015493037334322219, 7.86145885635573e-5, 2.0599035810421294e-9, 0.0], true, 12, 12)

The function solve is not exported so we call it through BifurcationKit.solve. The type ODEBifProblem is a shortcut to define a problem of the form $\dot x = F(x,p)$ ; the generic (autonomous) problem type is BifurcationProblem. Both are equivalent for the Newton algorithm above.

The returned object sol is a NonLinearSolution; it contains the solution u, the residuals at each iteration, the number of iterations and whether the algorithm converged (use BifurcationKit.converged(sol)).

Other flavours of newton

The Newton algorithm is also implemented for specific functionals, in which case a simplified entry point newton (exported) is provided which dispatches on the type of the problem:

• deflated problems: pass a DeflationOperator to solve in order to find other solutions than the known ones, see Deflated problems;
• periodic orbits computed by finite differences, shooting or collocation, see the corresponding pages Periodic orbits based on Trapezoidal rule, Periodic orbits based on the shooting method, Periodic orbits based on orthogonal collocation;
• travelling waves, see Freezing problems, symmetries and waves;
• refining branch points of codimension 2, see for example newton_hopf, newton_fold, newton_bt and the pages Fold / Hopf Continuation and Bogdanov-Takens refinement.

Example

The tutorial Temperature model presents the various jacobians (direct and iterative ones) in a concrete setting.