Asymptotic numerical method (ANM)

To access this algorithm, you have to use the package AsymptoticNumericalMethod.jl

The method ^[Rubbert],^{[Charpentier]} seeks a Taylor approximation of the solutions of

\[\mathbf F(X,p)=0\in\mathbb R^n\]

where $X\in\mathbb R^n, p\in\mathbb R$. The solution is found by solving

\[F(x(s),p(s))= 0\]

\[\langle x(s)-x(0),\frac{\partial u}{\partial s}(0)\rangle + (p(s)-p(0))\frac{\partial p}{\partial s}(0) = s\]

where

\[x(s)=\sum\limits_{k=0}^K x_ks^k,\quad p(s)=\sum\limits_{k=0}^K p_ks^k\]

for some user passed $K>0$. It gives

\[F(x(s),p(s)) = \sum\limits_{k=0}^K F_ks^k+h.o.t.\]

from which we deduce the equations $F_k=0$. We then find:

\[F_{1, x_1=Id, p_1=1}(x_k,p_k) = -F_{k,x_k=0,p_k=0}.\]

The validity range of the solution can be estimated by

\[r_K = \left(\frac{\epsilon}{\lvert\lvert F_K\lvert\lvert}\right)^{1/K}\]

This allows to iterate and find a sequence of series which spans the parameter range.

Implementation

The method is based on the package TaylorSeries.jl which makes it easy to manipulate Taylor series based on Automatic Differentiation.

Method

See AsymptoticNumericalMethod.ANM for more information.

Example

We provide an example of use. We define a BifurcationProblem as usual and pass the continuation algorithm ANM.

using AsymptoticNumericalMethod, Plots, Parameters, Setfield
using LinearAlgebra: norm
using BifurcationKit
const BK = BifurcationKit

norminf(x) = norm(x, Inf)

function F(x, p)
	@unpack α = p
	f = similar(x)

	f[1] = (-2x[1]+x[2]) + α * exp(x[1])
	f[2] = ( x[1]-2x[2]) + α * exp(x[2])

	return f
end

sol0 = zeros(2)
par = (α = 0.0, )
prob = BifurcationProblem(F, sol0, par, (@lens _.α); record_from_solution = (x,p) -> norminf(x))

┌─ Bifurcation Problem with uType Vector{Float64}
├─ Inplace:  false
├─ Symmetric: false
└─ Parameter: α

optanm = ContinuationPar(dsmin = 0.01, dsmax = 0.15, detect_bifurcation = 3, ds= 0.01, newton_options = NewtonPar(tol = 1e-9, verbose = false), n_inversion = 6, max_bisection_steps = 15, max_steps = 15, )

branm = continuation(prob, ANM(20, 1e-8), optanm, normC = norminf, verbosity = 2)

 ┌─ Curve type: EquilibriumCont
 ├─ 15 series of degree 20, tol = 1.0e-8
 ├─ Number of points: 301
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter α starts at 6.666666666666667e-5, ends at 1.1157858298101163e-9
 ├─ Algo: ANM
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1,       bp at α ≈ +0.36763519 ∈ (+0.36787944, +0.36787944), |δp|=2e-11, [converged], δ = ( 1,  0), step =  27, eigenelements in eig[ 28], ind_ev =   1
- #  2,       bp at α ≈ +0.14935757 ∈ (+0.14936120, +0.14936121), |δp|=8e-09, [converged], δ = ( 1,  0), step =  76, eigenelements in eig[ 77], ind_ev =   2

You can plot the result as usual:

plot(branm)

You can also show the radius of convergence of each series:

plot(branm, plotseries = true)

Finally, for each series, we ca evaluate the residual norm:

plot()
for ii in eachindex(branm.polU)
	s = LinRange(-0*branm.radius[ii], branm.radius[ii], 20)
	plot!([branm.polp[ii].(s)], [norminf(F(branm.polU[ii](_s), BK.setparam(prob,branm.polp[ii](_s)))) for _s in s], legend = false, linewidth=5)#, marker=:d)
end
title!("")

References