Krylov-Newton algorithm

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.

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 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.0000000000000018, 1.0], โ”Œโ”€ Bifurcation problem with uType Vector{Float64}
โ”œโ”€ Inplace: false
โ”œโ”€ Dimension: 10
โ”œโ”€ Symmetric: false
โ””โ”€ Parameter: p, [2.491196913123572, 4002.148882269823, 1185.4751271211933, 350.9059321261782, 103.6279551477223, 30.364595174994815, 8.668087834880135, 2.2724082338835228, 0.4585252864710755, 0.04396555817809969, 0.0006089818094604965, 1.2352515502864533e-7, 5.329070518200751e-15], true, 12, 12)

Example

The (basic) tutorial Temperature model presents all cases (direct and iterative ones).