π‘ Temperature model with ApproxFun, no AbstractArray
We reconsider the example Temperature model by relying on the package ApproxFun.jl which allows very precise function approximation. This is an interesting example because we have to change the scalar product of PALC for the method to work well.
This is one example where the state space, the space of solutions to the nonlinear equation, is not a subtype of
AbstractArray. See Required methods for custom arrays for more information.
Code for custom state
We start with some imports:
using ApproxFun, LinearAlgebrausing BifurcationKit, Plotsconst BK = BifurcationKitBK.set_plot_backend!(BK.BK_Plots()) # hideWe then need to add some methods not available in ApproxFun because the state space is not a subtype of AbstractArray:
# specific methods for ApproxFunimport Base: eltype, similar, copyto!, lengthimport LinearAlgebra: mul!, rmul!, axpy!, axpby!, dot, normsimilar(x::ApproxFun.Fun, T) = (copy(x))similar(x::ApproxFun.Fun) = copy(x)mul!(w::ApproxFun.Fun, v::ApproxFun.Fun, Ξ±) = (w .= Ξ± * v)eltype(x::ApproxFun.Fun) = eltype(x.coefficients)length(x::ApproxFun.Fun) = length(x.coefficients)dot(x::ApproxFun.Fun, y::ApproxFun.Fun) = sum(x * y)axpy!(a, x::ApproxFun.Fun, y::ApproxFun.Fun) = (y .= a * x + y)axpby!(a::Float64, x::ApproxFun.Fun, b::Float64, y::ApproxFun.Fun) = (y .= a * x + b * y)rmul!(y::ApproxFun.Fun, b::Float64) = (y.coefficients .*= b; y)rmul!(y::ApproxFun.Fun, b::Bool) = b == true ? y : (y.coefficients .*= 0; y)copyto!(x::ApproxFun.Fun, y::ApproxFun.Fun) = ( (x.coefficients = copy(y.coefficients);x))Problem formulation
We can easily write our functional with boundary conditions in a convenient manner using ApproxFun:
N(x; a = 0.5, b = 0.01) = 1 + (x + a*x^2)/(1 + b*x^2)dN(x; a = 0.5, b = 0.01) = (1-b*x^2+2*a*x)/(1+b*x^2)^2function F_chan(u, p) (;Ξ±, Ξ², Ξ) = p return [Fun(u(0.), domain(u)) - Ξ², Fun(u(1.), domain(u)) - Ξ², Ξ * u + Ξ± * N(u, b = Ξ²)]endfunction Jac_chan(u, p) (;Ξ±, Ξ², Ξ) = p return [Evaluation(u.space, 0.), Evaluation(u.space, 1.), Ξ + Ξ± * dN(u, b = Ξ²)]endWe want to call a Newton solver. We first need an initial guess and the Laplacian operator:
sol = Fun(x -> x * (1-x), Interval(0.0, 1.0))Ξ = Derivative(sol.space, 2)# set of parameterspar_af = (Ξ± = 3., Ξ² = 0.01, Ξ = Ξ)prob = BifurcationProblem(F_chan, sol, par_af, (@optic _.Ξ±); J = Jac_chan, plot_solution = (x, p; kwargs...) -> plot!(x; label = "l = $(length(x))", kwargs...))Finally, we need to provide some parameters for the Newton iterations. This is done by calling
optnewton = NewtonPar(tol = 1e-12, verbose = true)We call the Newton solver:
out = @time BK.solve(prob, Newton(), optnewton, normN = x -> norm(x, Inf64))and you should see
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββ Newton step residual linear iterations ββββββββββββββββ¬βββββββββββββββββββββββ¬βββββββββββββββββ€β 0 β 1.5707e+00 β 0 ββ 1 β 1.1546e-01 β 1 ββ 2 β 8.0149e-04 β 1 ββ 3 β 3.9038e-08 β 1 ββ 4 β 7.9049e-13 β 1 ββββββββββββββββ΄ββββββ-ββββββββββββββββ΄βββββββββββββββββ 0.103869 seconds (362.15 k allocations: 14.606 MiB)Continuation
We can also perform numerical continuation with respect to the parameter $\alpha$. Again, we need to provide some parameters for the continuation:
optcont = ContinuationPar(dsmin = 0.0001, dsmax = 0.05, ds= 0.005, p_max = 4.1, plot_every_step = 10, newton_options = NewtonPar(tol = 1e-8, max_iterations = 20, verbose = true), detect_bifurcation = 0, max_steps = 200)Then, we can call the continuation routine.
# we need a specific bordered linear solver# we use the BorderingBLS one to rely on ApproxFun.\br = continuation(prob, PALC(bls = BorderingBLS(solver = optnewton.linsolver, check_precision = false)), optcont, plot = true, plot_solution = (x, p; kwargs...) -> plot!(x; label = "l = $(length(x))", kwargs...), verbosity = 2, normC = x -> norm(x, Inf64))and you should see
However, if we do that, we'll see that it does not converge very well. The reason is that the default arc-length constraint (see Pseudo arclength continuation) is
\[N(x, p)=\frac{\theta}{\text { length}(x)}\left\langle x-x_{0}, d x_{0}\right\rangle+(1-\theta) \cdot\left(p-p_{0}\right) \cdot d p_{0}-d s=0\]
is tailored for vectors of fixed length. The $\frac{1}{length(x)}$ is added to properly balance the terms in the constraint. Thus, in BifurcationKit, the dot product is in fact (x,y) -> dot(x,y) / length(y).
But here, the vector space is provided with a custom dot product (see above) which depends on the domain, here Interval(0.0, 1.0). Hence, we want to change this constraint $N$ for the following:
\[N(x, p)={\theta}\left\langle x-x_{0}, d x_{0}\right\rangle+(1-\theta) \cdot\left(p-p_{0}\right) \cdot d p_{0}-d s=0.\]
This can be done as follows:
optcont = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, p_max = 4.1, plot_every_step = 10, newton_options = NewtonPar(tol = 1e-8, maxIter = 20, verbose = true), max_steps = 300, ΞΈ = 0.2, detect_bifurcation = 0)br = continuation(prob, PALC(bls=BorderingBLS(solver = optnewton.linsolver, check_precision = false)), optcont, plot = true, # specify the dot product used in PALC dotPALC = BK.DotTheta(dot), # we need a specific bordered linear solver # we use the BorderingBLS one to rely on ApproxFun.\ linear_algo = BorderingBLS(solver = DefaultLS(), check_precision = false), plot_solution = (x, p; kwargs...) -> plot!(x; label = "l = $(length(x))", kwargs...), verbosity = 2, normC = x -> norm(x, Inf64))