Update docs

docs/src/assets/radiipolynomial.css

@@ -17,70 +17,6 @@ html.theme--documenter-dark .content pre {
   border: 1px solid #282f2f;
   border-radius: 15px; }
 
-pre .copy-button {
-  opacity: 0;
-  transition: opacity 0.2s;
-  position: absolute;
-  right: 0em;
-  top: 0em;
-  padding: 0.5em;
-  width: 2.5em;
-  height: 2.5em;
-  background: transparent;
-  border: none;
-  font-family: "Font Awesome 5 Free";
-  color: #d5c4a1;
-  cursor: pointer;
-  text-align: center; }
-
-pre .copy-button:focus, pre .copy-button:hover {
-  opacity: 1;
-  background: transparent;
-  color: #d5c4a1; }
-
-pre .copy-button.success {
-  color: #259a12;
-  opacity: 1; }
-
-pre .copy-button.error {
-  color: #cb3c33;
-  opacity: 1; }
-
-pre:hover .copy-button {
-  opacity: 1; }
-
-html.theme--documenter-dark pre .copy-button {
-  opacity: 0;
-  transition: opacity 0.2s;
-  position: absolute;
-  right: 0em;
-  top: 0em;
-  padding: 0.5em;
-  width: 2.5em;
-  height: 2.5em;
-  background: transparent;
-  border: none;
-  font-family: "Font Awesome 5 Free";
-  color: #d5c4a1;
-  cursor: pointer;
-  text-align: center; }
-
-html.theme--documenter-dark pre .copy-button:focus, html.theme--documenter-dark pre .copy-button:hover {
-  opacity: 1;
-  background: transparent;
-  color: #d5c4a1; }
-
-html.theme--documenter-dark pre .copy-button.success {
-  color: #259a12;
-  opacity: 1; }
-
-html.theme--documenter-dark pre .copy-button.error {
-  color: #cb3c33;
-  opacity: 1; }
-
-html.theme--documenter-dark pre:hover .copy-button {
-  opacity: 1; }
-
 /* Start page layout */
 
 .single {

docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.md

@@ -4,7 +4,7 @@ In this example, we will prove the existence of a branch of equilibria of the Fi
 
 ```math
 \begin{cases}
-\displaystyle \frac{d}{dt} u(t) = f(\gamma, u(t)) \bydef \begin{pmatrix} u_1(t)(u_1(t) - a)(1 - u_1(t)) - u_2(t) \\ \varepsilon(u_1(t) - \gamma u_2(t)) \end{pmatrix},\\
+\displaystyle \frac{d}{dt} u(t) = f(u(t), \gamma) \bydef \begin{pmatrix} u_1(t)(u_1(t) - a)(1 - u_1(t)) - u_2(t) \\ \varepsilon(u_1(t) - \gamma u_2(t)) \end{pmatrix},\\
 u(0) = u_0 \in \mathbb{R}^2,
 \end{cases}
 ```
@@ -14,21 +14,28 @@ where ``a = 5`` and ``\varepsilon = 1``.
 The vector-field ``f`` and its Jacobian, denoted ``Df``, may be implemented as follows:
 
 ```@example fhn_pseudo_arclength
-function f(x)
+function f(u, γ)
     a, ϵ = 5, 1
-    γ, u₁, u₂ = x
+    u₁, u₂ = u
     return [u₁*(u₁ - a)*(1 - u₁) - u₂, ϵ*(u₁ - γ*u₂)]
 end
 
-function Df(x)
+function Dᵤf(u, γ)
     a, ϵ = 5, 1
-    γ, u₁, u₂ = x
-    return [0 a*(2u₁-1)+(2-3u₁)*u₁ -1 ; -ϵ*u₂ ϵ -ϵ*γ]
+    u₁, u₂ = u
+    return [a*(2u₁-1)+(2-3u₁)*u₁ -1
+            ϵ                    -ϵ*γ]
+end
+
+function ∂γf(u, γ)
+    a, ϵ = 5, 1
+    u₁, u₂ = u
+    return [0, -ϵ*u₂]
 end
 nothing # hide
 ```
 
-To this end, we use the [pseudo-arclength continuation](https://en.wikipedia.org/wiki/Numerical_continuation#Pseudo-arclength_continuation) and prove, at each step, that there exists a box, surrounding the linear numerical approximation, which contains the desired curve.
+We use the [pseudo-arclength continuation](https://en.wikipedia.org/wiki/Numerical_continuation#Pseudo-arclength_continuation) and retrieve a numerical approximation of the curve. By a contraction argument, we then prove that there exists a surrounding region that contains the desired curve.
 
 In a nutshell, the pseudo-arclength continuation consists in computing a sequence of numerical zeros of ``f``. Starting with an initial approximate zero ``x_\text{init} \in \mathbb{R}^3``, we retrieve an approximate tangent vector ``v`` to the curve at ``x_\text{init}`` by looking at ``\ker Df(x_\text{init})``. Then, our predictor for the next zero is set to ``w \bydef x_\text{init} + \delta v`` where ``\delta > 0`` represents the step size. The Newton's method is applied on the mapping ``F_\text{Newton} : \mathbb{R}^3 \to \mathbb{R}^3`` given by
 
@@ -45,9 +52,9 @@ The mapping ``F_\text{Newton}`` and its Jacobian may be implemented as follows:
 ```@example fhn_pseudo_arclength
 import LinearAlgebra: ⋅
 
-F(x, v, w) = [(x - w) ⋅ v ; f(x)]
+F(x, v, w) = [f(x[1:2], x[3]) ; (x - w) ⋅ v]
 
-DF(x, v) = [transpose(v) ; Df(x)]
+DF(x, v) = [Dᵤf(x[1:2], x[3]) ∂γf(x[1:2], x[3]) ; transpose(v)]
 nothing # hide
 ```
 
@@ -57,9 +64,18 @@ Next, we perform Newton's method:
 using RadiiPolynomial
 import LinearAlgebra: nullspace
 
-x_init = [2, 1.129171306613029, 0.564585653306514] # initial point on the branch of equilibria
+# initial point on the branch of equilibria
+
+γ_init = 2.0
 
-v = vec(nullspace(Df(x_init))) # initial tangent vector
+u_init = [1.1, 0.5]
+u_init, success = newton(u -> (f(u, γ_init), Dᵤf(u, γ_init)), u_init)
+
+# next point on the branch of equilibria
+
+x_init = [u_init ; γ_init]
+
+v = vec(nullspace([Dᵤf(x_init[1:2], x_init[3]) ∂γf(x_init[1:2], x_init[3])])) # initial tangent vector
 
 δ = 5e-2 # step size
 
@@ -69,70 +85,137 @@ x_final, success = newton(x -> (F(x, v, w), DF(x, v)), w)
 nothing # hide
 ```
 
-Once the Newton's method converged to some ``x_\text{final} \in \mathbb{R}^3``, we make a linear approximation of the curve of zeros
+Whenever Newton's method is successful, we proceed to the next iteration of the pseudo-arclength continuation by repeating the above strategy. Performing this sufficiently many times, we can construct an order ``N`` polynomial approximation of the curve of zeros:
 
 ```math
-x_0(s) \bydef x_\text{init} + s (x_\text{final} - x_\text{init}), \qquad \text{for all } s \in [0,1].
+\bar{x}(s) \bydef \bar{x}_0 + 2 \sum_{n = 1}^N \bar{x}_n \phi_n (s), \qquad \text{for all } s \in [-1,1],
 ```
 
-Define the mapping ``F : \mathbb{R}^3 \times [0,1] \to \mathbb{R}^3`` by
+where ``\phi_n`` are the [Chebyshev polynomials of the first kind](https://en.wikipedia.org/wiki/Chebyshev_polynomials).
+
+Define the mapping ``F : \mathbb{R}^3 \times [-1,1] \to \mathbb{R}^3`` by
 
 ```math
 F(x, s) \bydef
 \begin{pmatrix}
-(x - x_0(s)) \cdot v\\
-f(x)
-\end{pmatrix}
+f(x) \\
+(x - \bar{x}(s)) \cdot \bar{v}(s)
+\end{pmatrix},
 ```
 
-and the fixed-point operator ``T : \mathbb{R}^3 \times [0,1] \to \mathbb{R}^3`` by
+and the fixed-point operator ``T : \mathbb{R}^3 \times [-1,1] \to \mathbb{R}^3`` by
 
 ```math
-T(x, s) \bydef x - A F(x, s),
+T(x, s) \bydef x - A(s) F(x, s),
 ```
 
-where ``A : \mathbb{R}^3 \to \mathbb{R}^3`` is the injective operator corresponding to a numerical approximation of ``D_x F(x_0(s), s)^{-1}`` for all ``s \in [0, 1]``.
+where ``A(s) : \mathbb{R}^3 \to \mathbb{R}^3`` is the injective operator corresponding to a numerical approximation of ``D_x F(\bar{x}(s), s)^{-1}`` for all ``s \in [-1, 1]``.
 
-Let ``R > 0``. We use a uniform version of the second-order Radii Polynomial Theorem (cf. Section [Radii polynomial approach](@ref radii_polynomial_approach)) such that we need to estimate ``|T(x_0(s), s) - x_0(s)|_\infty``, ``|D_x T(x_0(s), s)|_\infty`` and ``\sup_{x \in \text{cl}( B_R(x_0(s)) )} |D_x^2 T(x, s)|_\infty`` for all ``s \in [0,1]``. In particular, we have
+Let ``R > 0``. We use a uniform version of the second-order Radii Polynomial Theorem (cf. Section [Radii polynomial approach](@ref radii_polynomial_approach)) such that we need to estimate ``\|T(\bar{x}(s), s) - \bar{x}(s)\|_1``, ``\|D_x T(\bar{x}(s), s)\|_1`` and ``\sup_{x \in \text{cl}( B_R(\bar{x}(s)) )} \|D_x^2 T(x, s)\|_1`` for all ``s \in [-1,1]``. In particular, we have
 
 ```math
-|T(x_0(s), s) - x_0(s)|_\infty = \left|A \begin{pmatrix} 0 \\ f(x_0(s)) \end{pmatrix} \right|_\infty, \qquad \text{for all } s \in [0,1].
+\|T(\bar{x}(s), s) - \bar{x}(s)\|_1 = \left\|A \begin{pmatrix} f(\bar{x}(s)) \\ 0 \end{pmatrix} \right\|_1, \qquad \text{for all } s \in [-1,1].
 ```
 
 The computer-assisted proof may be implemented as follows:
 
 ```@example fhn_pseudo_arclength
-R = 1e-1
+N = 500
+N_fft = nextpow(2, 2N + 1)
+npts = N_fft ÷ 2 + 1
+
+arclength = 15.0
+arclength_grid = [0.5 * arclength - 0.5 * cospi(2j/N_fft) * arclength for j ∈ 0:npts-1]
+x_grid = Vector{Vector{Float64}}(undef, npts)
+v_grid = Vector{Vector{Float64}}(undef, npts)
+
+# initialize
+
+direction = [0, 0, -1] # starts by decreasing the parameter
+x_grid[1] = x_init
+v_grid[1] = vec(nullspace([Dᵤf(x_grid[1][1:2], x_grid[1][3]) ∂γf(x_grid[1][1:2], x_grid[1][3])]))
+if direction ⋅ v_grid[1] < 0 # enforce direction
+    v_grid[1] .*= -1
+end
+
+# run continuation scheme
+
+for i ∈ 2:npts
+    δᵢ = arclength_grid[i] - arclength_grid[i-1]
 
-x₀_interval = interval.(x_init) .+ interval(0.0, 1.0) .* (interval.(x_final) .- interval.(x_init))
-x₀R_interval = interval.(x₀_interval, R; format = :midpoint)
+    wᵢ = x_grid[i-1] .+ δᵢ .* v_grid[i-1]
 
-F_interval = F(x₀_interval, v, x₀_interval)
-F_interval[1] = 0 # the first component is trivial by definition
-DF_interval = DF(x₀_interval, v)
+    x, success = newton(x -> (F(x, v_grid[i-1], wᵢ), DF(x, v_grid[i-1])), wᵢ; verbose = true)
+    success || error()
 
-A = inv(mid.(DF_interval))
+    x_grid[i] = x
+    v_grid[i] = vec(nullspace([Dᵤf(x_grid[i][1:2], x_grid[i][3]) ∂γf(x_grid[i][1:2], x_grid[i][3])]))
+    if v_grid[i-1] ⋅ v_grid[i] < 0 # keep the same direction
+        v_grid[i] .*= -1
+    end
+end
+
+# construct the approximations
+
+function cheb2grid(x::VecOrMat{<:Sequence}, N_fft)
+    vals = fft.(x, N_fft)
+    return [real.(getindex.(vals, i)) for i ∈ eachindex(vals[1])]
+end
 
-Y = norm(Sequence(A * F_interval), Inf)
+grid2cheb(x_fft::Vector{<:Vector}, N) =
+    [rifft!(complex.(getindex.(x_fft, i)), Chebyshev(N)) for i ∈ eachindex(x_fft[1])]
 
-Z₁ = opnorm(LinearOperator(A * DF_interval - I), Inf)
+grid2cheb(x_fft::Vector{<:Matrix}, N) =
+    [rifft!(complex.(getindex.(x_fft, i, j)), Chebyshev(N)) for i ∈ axes(x_fft[1], 1), j ∈ axes(x_fft[1], 2)]
+
+x_fft = [reverse(x_grid) ; x_grid[begin+1:end-1]]
+x̄ = map(x -> interval.(x), grid2cheb(x_fft, N))
+
+v_fft = [reverse(v_grid) ; v_grid[begin+1:end-1]]
+v̄ = map(v -> interval.(v), grid2cheb(v_fft, N))
+
+A = map(A -> interval.(A), grid2cheb(inv.(DF.(x_fft, v_fft)), N))
+
+# AF is a polynomial with respect to s of order 3N
+
+N3 = 3N
+N3_fft = nextpow(2, 2N3 + 1)
+
+AF_fft = cheb2grid(A, N3_fft) .* F.(cheb2grid(x̄, N3_fft), cheb2grid(v̄, N3_fft), cheb2grid(x̄, N3_fft))
+AF = grid2cheb(AF_fft, N3)
+
+Y = norm(norm.(AF, 1), 1)
+
+# ADF is a polynomial with respect to s of order 2N
+
+N2 = 2N
+N2_fft = nextpow(2, 2N2 + 1)
+
+I_ADF_fft = Ref(I) .- cheb2grid(A, N2_fft) .* DF.(cheb2grid(x̄, N2_fft), cheb2grid(v̄, N2_fft))
+I_ADF = grid2cheb(I_ADF_fft, N2)
+
+Z₁ = opnorm(norm.(I_ADF, 1), 1)
+
+#
+
+R = 2sup(Y)
 
 a, ϵ = 5, 1
-Z₂ = opnorm(LinearOperator(interval.(A)), Inf) * max(2abs(a) + 2 + 6(abs(x₀_interval[2]) + R), 2abs(ϵ))
+Z₂ = opnorm(norm.(A, 1), 1) * max(abs(2a + 2) + 6(norm(x̄[1], 1) + R) + abs(ϵ), abs(ϵ))
+
+#
 
 setdisplay(:full)
 
 interval_of_existence(Y, Z₁, Z₂, R)
 ```
 
-Whenever the proof is successful, we proceed to the next iteration of the pseudo-arclength continuation and repeat the above strategy.
-
-The following animation[^1] shows the successive steps of a rigorous pseudo-arclength continuation of equilibria of the FitzHugh-Nagumo model.
+The following figure[^1] shows the numerical approximation of the proven branch of equilibria of the FitzHugh-Nagumo model.
 
 [^1]: S. Danisch and J. Krumbiegel, [Makie.jl: Flexible high-performance data visualization for Julia](https://doi.org/10.21105/joss.03349), *Journal of Open Source Software*, **6** (2021), 3349.
 
 ```@raw html
 <video width="800" height="400" controls autoplay>
-  <source src="../fhn_pseudo_arclength.mp4" type="video/mp4">
+  <source src="../fhn_pseudo_arclength.svg" type="video/mp4">
 </video>
 ```