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 Requested methods for Custom State.
Different Jacobians
There are basically two ways to specify the jacobian:
- Matrix based
- 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 BifurcationKit
F(x, p) = x.^3 .- 1
x0 = 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"#152#154"{typeof(Main.F)}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Float64, 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.822974993391147, 3.6720618744401145e6, 1.0880180731059401e6, 322375.46530502324, 95518.39633243749, 28301.487022406745, 8385.365784744188, 2484.292935472064, 735.8271824298597, 217.7637624767373, 64.26434372541303, 18.78544139894608, 5.318123666180427, 1.3525810845206885, 0.24265626976332633, 0.015451981769049716, 7.824431957104139e-5, 2.040547064652287e-9, 0.0], true, 18, 18)Example
The (basic) tutorial Temperature model presents all cases (direct and iterative ones).