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 API.
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 = BifurcationProblem(F, x0, nothing)
sol = newton(prob, NewtonPar(verbose = true))
NonLinearSolution{Vector{Float64}, BifurcationProblem{BifFunction{typeof(Main.F), BifurcationKit.var"#8#24", Nothing, BifurcationKit.var"#6#22", Nothing, BifurcationKit.var"#11#28"{BifurcationKit.var"#d1Fad#26"}, BifurcationKit.var"#13#30", BifurcationKit.var"#15#32", BifurcationKit.var"#17#34", Bool, Float64}, Vector{Float64}, Nothing, Setfield.IdentityLens, typeof(BifurcationKit.plot_default), typeof(BifurcationKit.record_sol_default), typeof(BifurcationKit.save_solution_default)}, Vector{Float64}, Int64}([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0000000000000002], ┌─ Bifurcation Problem with uType Vector{Float64}
├─ Inplace: false
├─ Symmetric: false
└─ Parameter: p, [2.2319442467943955, 2.099437757410674e9, 6.220556315661258e8, 1.8431277946403712e8, 5.4611193656011015e7, 1.61810941573366e7, 4.794398009581227e6, 1.4205621139500395e6, 420907.0337631312, 124712.93593034455 … 284.43439995053757, 84.01837898037451, 24.637694191379186, 7.0495224197447115, 1.8576666121761036, 0.3734593743839818, 0.03282662928849622, 0.00034655126182991935, 4.0017183566831704e-8, 6.661338147750939e-16], true, 23, 23)
Example
The (basic) tutorial Temperature model presents all cases (direct and iterative ones).