diff --git a/docs/src/assets/radiipolynomial.css b/docs/src/assets/radiipolynomial.css
index ce66b4f3..101dc168 100644
--- a/docs/src/assets/radiipolynomial.css
+++ b/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 {
diff --git a/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.md b/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.md
index 233be1d0..43a57b0b 100644
--- a/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.md
+++ b/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
```
diff --git a/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.mp4 b/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.mp4
deleted file mode 100644
index 76f1a02e..00000000
Binary files a/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.mp4 and /dev/null differ
diff --git a/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.svg b/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.svg
new file mode 100644
index 00000000..ebeb6de6
--- /dev/null
+++ b/docs/src/examples/finite_dimensional_proofs/fhn_pseudo_arclength.svg
@@ -0,0 +1,148 @@
+
+
diff --git a/docs/src/examples/finite_dimensional_proofs/spiderweb.md b/docs/src/examples/finite_dimensional_proofs/spiderweb.md
index 0a0d911e..bc183805 100644
--- a/docs/src/examples/finite_dimensional_proofs/spiderweb.md
+++ b/docs/src/examples/finite_dimensional_proofs/spiderweb.md
@@ -135,7 +135,7 @@ x₀, success = newton(x -> (F(x, m₀, m, λ, l), DF(x, m₀, m, λ, l)), x₀)
nothing # hide
```
-Let ``R > 0``. According to the first-order Radii Polynomial Theorem (cf. Section [Radii polynomial approach](@ref radii_polynomial_approach)), we need to estimate ``|T(x_0) - x_0|_\infty`` and ``\sup_{x \in \text{cl}( B_R(x_0) )} |DT(x)|_\infty`` which can be readily computed with interval arithmetic.
+Let ``R > 0``. According to the first-order Radii Polynomial Theorem (cf. Section [Radii polynomial approach](@ref radii_polynomial_approach)), we need to estimate ``\|T(x_0) - x_0\|_\infty`` and ``\sup_{x \in \text{cl}( B_R(x_0) )} \|DT(x)\|_\infty`` which can be readily computed with interval arithmetic.
The computer-assisted proof may be implemented as follows:
diff --git a/docs/src/examples/infinite_dimensional_proofs/dde/ikeda_W_u.md b/docs/src/examples/infinite_dimensional_proofs/dde/ikeda_W_u.md
index 187c08d8..d281fa1a 100644
--- a/docs/src/examples/infinite_dimensional_proofs/dde/ikeda_W_u.md
+++ b/docs/src/examples/infinite_dimensional_proofs/dde/ikeda_W_u.md
@@ -3,7 +3,7 @@
In this example, we will rigorously compute the unstable manifolds of the equilibria for the cubic Ikeda equation
```math
-\frac{d}{dt} u(t) = f(u(t), u(t-\tau)) := u(t-\tau) - u(t-\tau)^3.
+\frac{d}{dt} u(t) = f(u(t), u(t-\tau)) \bydef u(t-\tau) - u(t-\tau)^3.
```
The linearization at some equilibrium ``c \in \mathbb{R}`` yields
@@ -15,7 +15,7 @@ The linearization at some equilibrium ``c \in \mathbb{R}`` yields
The right-hand side of the above equation is an infinite dimensional endomorphism acting on ``C([-\tau, 0], \mathbb{R})``. Its compactness guarantees that the spectrum is comprised of eigenvalues accumulating at ``0``; in particular, there are finitely many eigenvalues whose real parts are strictly positive. As a matter of fact, an eigenvector ``\xi \in C([-\tau, 0], \mathbb{C})`` associated with an eigenvalue ``\lambda \in \mathbb{C}`` is given by ``\xi(s) = e^{s \lambda} \xi(0)``, for all ``s \in [-\tau, 0]`` and ``\xi(0) \neq 0``, such that
```math
-\Psi(\lambda) := \lambda - (1 - 3c^2) e^{-\tau \lambda} = 0.
+\Psi(\lambda) \bydef \lambda - (1 - 3c^2) e^{-\tau \lambda} = 0.
```
The characteristic function ``\Psi`` and its derivative with respect to ``\lambda``, denoted ``D\Psi``, may be implemented as follows:
@@ -38,10 +38,10 @@ using RadiiPolynomial
R = 1e-14
-τ = interval(1.59)
+τ = I"1.59"
-Y = abs(Ψ(interval(λ̄₀), 0.0, τ))
-Z₁ = abs(1 - DΨ(λ̄₀, 0.0, mid(τ)) \ DΨ(interval(λ̄₀, R; format = :midpoint), 0.0, τ))
+Y = abs(Ψ(interval(λ̄₀), interval(0), τ))
+Z₁ = abs(1 - interval(DΨ(λ̄₀, 0, mid(τ))) \ DΨ(interval(λ̄₀, R; format = :midpoint), interval(0), τ))
ϵ₀ = inf(interval_of_existence(Y, Z₁, R))
λ₀ = interval(λ̄₀, ϵ₀; format = :midpoint)
@@ -55,8 +55,8 @@ Similarly, for the equilibria ``c = 1`` and ``c = -1``, we may use the same stra
```@example ikeda_W_u
λ̄₁, success = newton(λ -> (Ψ(λ, 1.0, 1.59), DΨ(λ, 1.0, 1.59)), 0.3+1.0im)
-Y = abs(Ψ(interval(λ̄₁), 1.0, τ))
-Z₁ = abs(1 - DΨ(λ̄₁, 1.0, mid(τ)) \ DΨ(interval(λ̄₁, R; format = :midpoint), 1.0, τ))
+Y = abs(Ψ(interval(λ̄₁), interval(1), τ))
+Z₁ = abs(1 - interval(DΨ(λ̄₁, 1, τ)) \ DΨ(interval(λ̄₁, R; format = :midpoint), interval(1), τ))
ϵ₁ = inf(interval_of_existence(Y, Z₁, R))
λ₁ = interval(λ̄₁, ϵ₁; format = :midpoint)
@@ -65,21 +65,21 @@ setdisplay(:full)
λ₁
```
-Let ``\lambda_1, \dots, \lambda_d`` be the unstable eigenvalues and ``\xi_1, \dots, \xi_d`` the respective eigenvectors. Denote by ``\Lambda : \mathbb{C}^d \to \mathbb{C}^d`` the diagonal matrix such that ``\Lambda_{i,i} := \lambda_i``; also, denote by ``\Xi : \mathbb{C}^d \to C([-\tau, 0], \mathbb{C})`` the matrix whose ``i``-th column is the eigenvector ``\xi_i``.
+Let ``\lambda_1, \dots, \lambda_d`` be the unstable eigenvalues and ``\xi_1, \dots, \xi_d`` the respective eigenvectors. Denote by ``\Lambda : \mathbb{C}^d \to \mathbb{C}^d`` the diagonal matrix such that ``\Lambda_{i,i} \bydef \lambda_i``; also, denote by ``\Xi : \mathbb{C}^d \to C([-\tau, 0], \mathbb{C})`` the matrix whose ``i``-th column is the eigenvector ``\xi_i``.
Let
```math
-X := \left\{ \{ x_\alpha \}_{\alpha_1 + \ldots + \alpha_d \ge 0} \in \mathbb{C}^{(\mathbb{N} \cup \{0\})^d} \, : \, | x |_X := \sum_{\alpha_1 + \ldots + \alpha_d \ge 0} |x_\alpha| < +\infty \right\}
+X \bydef \left\{ \{ x_\alpha \}_{\alpha_1 + \ldots + \alpha_d \ge 0} \in \mathbb{C}^{(\mathbb{N} \cup \{0\})^d} \, : \, \| x \|_X \bydef \sum_{\alpha_1 + \ldots + \alpha_d \ge 0} |x_\alpha| < +\infty \right\}
```
and ``* : X \times X \to X`` be the Cauchy product given by
```math
-x * y := \left\{ \sum_{\beta_1 + \ldots + \beta_d \ge 0}^\alpha x_{\alpha - \beta} y_\beta \right\}_{\alpha_1 + \ldots + \alpha_d \ge 0}, \qquad \text{for all } x, y \in X.
+x * y \bydef \left\{ \sum_{\beta_1 + \ldots + \beta_d \ge 0}^\alpha x_{\alpha - \beta} y_\beta \right\}_{\alpha_1 + \ldots + \alpha_d \ge 0}, \qquad \text{for all } x, y \in X.
```
-For any sequence ``x \in X``, the Taylor series ``\sum_{\alpha_1 + \ldots + \alpha_d \ge 0} x_\alpha \sigma^\alpha`` defines an analytic function in ``C^\omega(\mathbb{D}^d, \mathbb{C})`` where ``\mathbb{D} := \{ z \in \mathbb{C} \, : \, |z| \le 1 \}``; while the Cauchy product ``*`` corresponds to the product of Taylor series in sequence space.
+For any sequence ``x \in X``, the Taylor series ``\sum_{\alpha_1 + \ldots + \alpha_d \ge 0} x_\alpha \sigma^\alpha`` defines an analytic function in ``C^\omega(\mathbb{D}^d, \mathbb{C})`` where ``\mathbb{D} \bydef \{ z \in \mathbb{C} \, : \, |z| \le 1 \}``; while the Cauchy product ``*`` corresponds to the product of Taylor series in sequence space.
The Banach space ``X`` is a suitable space to represent a parameterization of the unstable manifold. Indeed, it is a standard result from DDE theory that analytic vector fields yield analytic unstable manifolds of equilibria. In the context of this example, it holds that the unstable manifold is parameterized by an analytic function ``P : \mathbb{C}^d \to C([-\tau, 0], \mathbb{C})`` satisfying ``\frac{d}{ds} [P(\sigma)](s) = [DP(\sigma) \Lambda \sigma](s)`` along with ``[DP(\sigma) \Lambda \sigma](0) = f([P (\sigma)](0), [P(\sigma)](-\tau))``.[^1]
@@ -94,7 +94,7 @@ In terms of the Taylor coefficients, the previous equalities yield
where ``\tilde{x} \in X`` is given component-wise by
```math
-\tilde{x}_\alpha :=
+\tilde{x}_\alpha \bydef
\begin{cases}
c, & \alpha_1 = \ldots = \alpha_d = 0,\\
\xi_1, & \alpha_1 = 1, \alpha_2 = \ldots = \alpha_d = 0,\\
@@ -110,12 +110,12 @@ For the equilibrium ``c = 0``, we may implement the ``1``-dimensional recurrence
```@example ikeda_W_u
n₀ = 85
-x̃₀ = Sequence(Taylor(n₀), zeros(Interval{Float64}, n₀+1))
-x̃₀[1] = 5.0
+x̃₀ = zeros(Interval{Float64}, Taylor(n₀))
+x̃₀[1] = interval(5)
ỹ₀ = copy(x̃₀)
ỹ₀[1] *= exp(-τ * λ₀)
for α ∈ 2:n₀
- x̃₀[α] = -Ψ(α*λ₀, 0.0, τ) \ pow_bar(Sequence(Taylor(α), view(ỹ₀, 0:α)), 3)[α]
+ x̃₀[α] = -Ψ(α*λ₀, interval(0), τ) \ pow_bar(Sequence(Taylor(α), view(ỹ₀, 0:α)), 3)[α]
ỹ₀[α] = x̃₀[α] * exp(-τ * α*λ₀)
end
```
@@ -124,15 +124,15 @@ Similarly, for the equilibrium ``c = 1``, we may implement the ``2``-dimensional
```@example ikeda_W_u
n₁ = 25
-x̃₁ = Sequence(Taylor(n₁) ⊗ Taylor(n₁), zeros(Complex{Interval{Float64}}, (n₁+1)^2))
-x̃₁[(0,0)] = 1.0
-x̃₁[(1,0)] = x̃₁[(0,1)] = 0.35
+x̃₁ = zeros(Complex{Interval{Float64}}, Taylor(n₁) ⊗ Taylor(n₁))
+x̃₁[(0,0)] = interval(1)
+x̃₁[(1,0)] = x̃₁[(0,1)] = interval(0.35)
ỹ₁ = copy(x̃₁)
ỹ₁[(1,0)] *= exp(-τ * λ₁)
ỹ₁[(0,1)] *= exp(-τ * conj(λ₁))
for α₂ ∈ 0:n₁, α₁ ∈ 0:n₁-α₂
if α₁ + α₂ ≥ 2
- x̃₁[(α₁,α₂)] = -Ψ(α₁*λ₁ + α₂*conj(λ₁), 1.0, τ) \ pow_bar(Sequence(Taylor(α₁) ⊗ Taylor(α₂), view(ỹ₁, (0:α₁, 0:α₂))), 3)[(α₁,α₂)]
+ x̃₁[(α₁,α₂)] = -Ψ(α₁*λ₁ + α₂*conj(λ₁), interval(1), τ) \ pow_bar(Sequence(Taylor(α₁) ⊗ Taylor(α₂), view(ỹ₁, (0:α₁, 0:α₂))), 3)[(α₁,α₂)]
ỹ₁[(α₁,α₂)] = x̃₁[(α₁,α₂)] * exp(-τ * (α₁*λ₁ + α₂*conj(λ₁)))
end
end
@@ -141,35 +141,35 @@ end
Consider the truncation operator
```math
-(\pi^n x)_\alpha :=
+(\Pi_n x)_\alpha \bydef
\begin{cases} x_\alpha, & \alpha_1 + \ldots + \alpha_d \le n,\\
0, & \alpha_1 + \ldots + \alpha_d > n,
\end{cases}
\qquad \text{for all } x \in X,
```
-as well as the complementary operator ``\pi^{\infty(n)} := I - \pi^n``.
+as well as the complementary operator ``\Pi_{\infty(n)} \bydef I - \Pi_n``.
-Given that ``\pi^n \tilde{x}`` is a finite sequence of known Taylor coefficients, it follows that the remaining coefficients are a fixed-point of the mapping ``T : \pi^{\infty(n)} X \to \pi^{\infty(n)} X`` given component-wise by
+Given that ``\Pi_n \tilde{x}`` is a finite sequence of known Taylor coefficients, it follows that the remaining coefficients are a fixed-point of the mapping ``T : \Pi_{\infty(n)} X \to \Pi_{\infty(n)} X`` given component-wise by
```math
-( T(h) )_\alpha :=
+( T(h) )_\alpha \bydef
\begin{cases}
0, & \alpha_1 + \ldots + \alpha_d \le n,\\
-\Psi(\alpha_1 \lambda_1 + \ldots + \alpha_d \lambda_d)^{-1} \left( -e^{-\tau (\alpha_1 \lambda_1 + \ldots + \alpha_d \lambda_d)} [(\pi^n \tilde{x} +h)*(\pi^n \tilde{x} +h)*(\pi^n \tilde{x} +h)]_{h_\alpha = 0} \right)_\alpha, & \alpha_1 + \ldots + \alpha_d > n.
+\Psi(\alpha_1 \lambda_1 + \ldots + \alpha_d \lambda_d)^{-1} \left( -e^{-\tau (\alpha_1 \lambda_1 + \ldots + \alpha_d \lambda_d)} [(\Pi_n \tilde{x} +h)*(\Pi_n \tilde{x} +h)*(\Pi_n \tilde{x} +h)]_{h_\alpha = 0} \right)_\alpha, & \alpha_1 + \ldots + \alpha_d > n.
\end{cases}
```
-Let ``R > 0``. Since ``T \in C^1(\pi^{\infty(n)} X, \pi^{\infty(n)} X)`` we may use the [first-order Radii Polynomial Theorem](@ref first_order_RPT) for which we use the estimates
+Let ``R > 0``. Since ``T \in C^1(\Pi_{\infty(n)} X, \Pi_{\infty(n)} X)`` we may use the [first-order Radii Polynomial Theorem](@ref first_order_RPT) for which we use the estimates
```math
\begin{aligned}
-|T(0)|_X &\le \max_{\mu \in \{\Re(\lambda_1), \dots, \Re(\lambda_d)\}} \frac{1}{(n+1)\mu - |1 - 3c^2| e^{-τ (n+1)\mu}} |\pi^{\infty(n)} (\pi^n \tilde{y}*\pi^n \tilde{y}*\pi^n \tilde{y})|_X,\\
-\sup_{h \in \text{cl}( B_R(0) )} |DT(h)|_{\mathscr{B}(X, X)} &\le \max_{\mu \in \{\Re(\lambda_1), \dots, \Re(\lambda_d)\}} \frac{3}{(n+1)\mu - |1 - 3c^2| e^{-τ (n+1)\mu}} (|\pi^n \tilde{y}|_X + R)^2,
+\|T(0)\|_X &\le \max_{\mu \in \{\Re(\lambda_1), \dots, \Re(\lambda_d)\}} \frac{1}{(n+1)\mu - |1 - 3c^2| e^{-τ (n+1)\mu}} \|\Pi_{\infty(n)} (\Pi_n \tilde{y}*\Pi_n \tilde{y}*\Pi_n \tilde{y})\|_X,\\
+\sup_{h \in \text{cl}( B_R(0) )} \|DT(h)\|_{\mathscr{B}(X, X)} &\le \max_{\mu \in \{\Re(\lambda_1), \dots, \Re(\lambda_d)\}} \frac{3}{(n+1)\mu - |1 - 3c^2| e^{-τ (n+1)\mu}} (\|\Pi_n \tilde{y}\|_X + R)^2,
\end{aligned}
```
-where ``\tilde{y} := \left\{ \tilde{x}_\alpha e^{-\tau (\alpha_1 \lambda_1 + \ldots + \alpha_d \lambda_d)} \right\}_{\alpha_1 + \ldots + \alpha_d \ge 0}``.
+where ``\tilde{y} \bydef \left\{ \tilde{x}_\alpha e^{-\tau (\alpha_1 \lambda_1 + \ldots + \alpha_d \lambda_d)} \right\}_{\alpha_1 + \ldots + \alpha_d \ge 0}``.
The computer-assisted proof for the ``1``-dimensional unstable manifold of ``c = 0`` may be implemented as follows:
@@ -179,12 +179,12 @@ X = ℓ¹()
R = 1e-12
tail_ỹ₀³ = ỹ₀ ^ 3
-tail_ỹ₀³[0:n₀] .= 0
-C₀ = (n₀+1) * λ₀ - exp(-τ * (n₀+1) * λ₀)
+tail_ỹ₀³[0:n₀] .= interval(0)
+C₀ = interval(n₀+1) * λ₀ - exp(-τ * interval(n₀+1) * λ₀)
Y = C₀ \ norm(tail_ỹ₀³, X)
-Z₁ = C₀ \ (3(norm(ỹ₀, X) + R)^2)
+Z₁ = C₀ \ (interval(3) * (norm(ỹ₀, X) + R)^2)
# error bound for the Taylor coefficients of order α > 85 of the parameterization on the domain [-1, 1]
@@ -198,13 +198,13 @@ Similarly, the computer-assisted proof for the ``2``-dimensional unstable manifo
```@example ikeda_W_u
tail_ỹ₁³ = ỹ₁ ^ 3
for α₂ ∈ 0:n₁, α₁ ∈ 0:n₁-α₂
- tail_ỹ₁³[(α₁,α₂)] = 0
+ tail_ỹ₁³[(α₁,α₂)] = interval(0)
end
-C₁ = (n₁+1) * real(λ₁) - 2exp(-τ * (n₁+1) * real(λ₁))
+C₁ = interval(n₁+1) * real(λ₁) - interval(2) * exp(-τ * interval(n₁+1) * real(λ₁))
Y = C₁ \ norm(tail_ỹ₁³, X)
-Z₁ = C₁ \ (3(norm(ỹ₁, X) + R)^2)
+Z₁ = C₁ \ (interval(3) * (norm(ỹ₁, X) + R)^2)
# error bound for the Taylor coefficients of order α₁ + α₂ > 25 of the parameterization on the domain 𝔻²
diff --git a/docs/src/examples/infinite_dimensional_proofs/ode/logistic_ivp.md b/docs/src/examples/infinite_dimensional_proofs/ode/logistic_ivp.md
index d9f28a90..3a657aa0 100644
--- a/docs/src/examples/infinite_dimensional_proofs/ode/logistic_ivp.md
+++ b/docs/src/examples/infinite_dimensional_proofs/ode/logistic_ivp.md
@@ -4,7 +4,7 @@ In this example, we will prove the existence of a solution of the logistic equat
```math
\begin{cases}
-\displaystyle \frac{d}{dt} u(t) = f(u(t)) := u(t)(1 - u(t)),\\
+\displaystyle \frac{d}{dt} u(t) = f(u(t)) \bydef u(t)(1 - u(t)),\\
u(0) = 1/2.
\end{cases}
```
@@ -12,25 +12,25 @@ u(0) = 1/2.
Let ``\nu > 0``,
```math
-X := \left\{ \{ x_\alpha \}_{\alpha \ge 0} \in \mathbb{R}^{\mathbb{N} \cup \{0\}} \, : \, | x |_X := \sum_{\alpha \ge 0} |x_\alpha| \nu^\alpha < +\infty \right\}
+\ell^1_{\nu, \mathbb{N}} \bydef \left\{ \{ x_\alpha \}_{\alpha \ge 0} \in \mathbb{R}^\mathbb{N} \, : \, \| x \|_{\ell^1_{\nu, \mathbb{N}}} \bydef \sum_{\alpha \ge 0} |x_\alpha| \nu^\alpha < +\infty \right\}
```
-and ``* : X \times X \to X`` be the Cauchy product given by
+and ``* : \ell^1_{\nu, \mathbb{N}} \times \ell^1_{\nu, \mathbb{N}} \to \ell^1_{\nu, \mathbb{N}}`` be the Cauchy product given by
```math
-x * y := \left\{ \sum_{\beta = 0}^\alpha x_{\alpha - \beta} y_\beta \right\}_{\alpha \ge 0}, \qquad \text{for all } x, y \in X.
+x * y \bydef \left\{ \sum_{\beta = 0}^\alpha x_{\alpha - \beta} y_\beta \right\}_{\alpha \ge 0}, \qquad \text{for all } x, y \in \ell^1_{\nu, \mathbb{N}}.
```
-For any sequence ``x \in X``, the Taylor series ``\sum_{\alpha \ge 0} x_\alpha t^\alpha`` defines an analytic function in ``C^\omega([-\nu, \nu], \mathbb{R})``; while the Cauchy product ``*`` corresponds to the product of Taylor series in sequence space.
+For any sequence ``x \in \ell^1_{\nu, \mathbb{N}}``, the Taylor series ``\sum_{\alpha \ge 0} x_\alpha t^\alpha`` defines an analytic function in ``C^\omega([-\nu, \nu], \mathbb{R})``; while the Cauchy product ``*`` corresponds to the product of Taylor series in sequence space.
-The Banach space ``X`` is a suitable space to represent a solution of the logistic equation. Indeed, it is a standard result from ODE theory that analytic vector fields yield analytic solutions.[^1]
+The Banach space ``\ell^1_{\nu, \mathbb{N}}`` is a suitable space to represent a solution of the logistic equation. Indeed, it is a standard result from ODE theory that analytic vector fields yield analytic solutions.[^1]
[^1]: A. Hungria, J.-P. Lessard and J. D. Mireles James, [Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach](https://doi.org/10.1090/mcom/3046), *Mathematics of Computation*, **85** (2016), 1427-1459.
-It follows that the sequence of coefficients of a Taylor series solving the initial value problem is a zero of the mapping ``F : X \to X`` given component-wise by
+It follows that the sequence of coefficients of a Taylor series solving the initial value problem is a zero of the mapping ``F : \ell^1_{\nu, \mathbb{N}} \to \ell^1_{\nu, \mathbb{N}}`` given component-wise by
```math
-( F(x) )_\alpha :=
+( F(x) )_\alpha \bydef
\begin{cases}
x_0 - 1/2, & \alpha = 0,\\
\alpha x_\alpha - (x*(1 - x))_{\alpha-1}, & \alpha \ge 1.
@@ -65,13 +65,13 @@ end
nothing # hide
```
-Consider the fixed-point operator ``T : X \to X`` defined by
+Consider the fixed-point operator ``T : \ell^1_{\nu, \mathbb{N}} \to \ell^1_{\nu, \mathbb{N}}`` defined by
```math
-T(x) := x - A F(x),
+T(x) \bydef x - A F(x),
```
-where ``A : X \to X`` is an injective operator corresponding to an approximation of ``DF(\bar{x})^{-1}`` for some numerical zero ``\bar{x} \in X`` of ``F``.
+where ``A : \ell^1_{\nu, \mathbb{N}} \to \ell^1_{\nu, \mathbb{N}}`` is an injective operator corresponding to an approximation of ``DF(\bar{x})^{-1}`` for some numerical zero ``\bar{x} \in \ell^1_{\nu, \mathbb{N}}`` of ``F``.
Given an initial guess, the numerical zero ``\bar{x}`` of ``F`` may be obtained by Newton's method:
@@ -84,26 +84,26 @@ x̄, success = newton!((F, DF, x) -> (F!(F, x), DF!(DF, x)), x̄)
nothing # hide
```
-Let ``R > 0``. Since ``T \in C^2(X, X)`` we may use the [second-order Radii Polynomial Theorem](@ref second_order_RPT) such that we need to estimate ``|T(\bar{x}) - \bar{x}|_X``, ``|DT(\bar{x})|_{\mathscr{B}(X, X)}`` and ``\sup_{x \in \text{cl}( B_R(\bar{x}) )} |D^2T(x)|_{\mathscr{B}(X^2, X)}``.
+Let ``R > 0``. Since ``T \in C^2(\ell^1_{\nu, \mathbb{N}}, \ell^1_{\nu, \mathbb{N}})`` we may use the [second-order Radii Polynomial Theorem](@ref second_order_RPT) such that we need to estimate ``\|T(\bar{x}) - \bar{x}\|_{\ell^1_{\nu, \mathbb{N}}}``, ``\|DT(\bar{x})\|_{\mathscr{B}(\ell^1_{\nu, \mathbb{N}}, \ell^1_{\nu, \mathbb{N}})}`` and ``\sup_{x \in \text{cl}( B_R(\bar{x}) )} \|D^2T(x)\|_{\mathscr{B}(\ell^1_{\nu, \mathbb{N}}, \mathscr{B}(\ell^1_{\nu, \mathbb{N}}, \ell^1_{\nu, \mathbb{N}}))}``.
To this end, consider the truncation operator
```math
-(\pi^n x)_\alpha :=
+(\Pi_n x)_\alpha \bydef
\begin{cases} x_\alpha, & \alpha \le n,\\
0, & \alpha > n,
-\end{cases} \qquad \text{for all } x \in X,
+\end{cases} \qquad \text{for all } x \in \ell^1_{\nu, \mathbb{N}},
```
-as well as the complementary operator ``\pi^{\infty(n)} := I - \pi^n``.
+as well as the complementary operator ``\Pi_{\infty(n)} \bydef I - \Pi_n``.
Thus, we have
```math
\begin{aligned}
-|T(\bar{x}) - \bar{x}|_X &\le |\pi^n A \pi^n F(\bar{x})|_X + \frac{1}{n+1} |\pi^{\infty(n)} F(\bar{x})|_X,\\
-|DT(\bar{x})|_{\mathscr{B}(X, X)} &\le |\pi^n A \pi^n DF(\bar{x}) \pi^n - \pi^n|_{\mathscr{B}(X, X)} + \frac{\nu}{n+1} |2\bar{x} - 1|_X,\\
-\sup_{x \in \text{cl}( B_R(\bar{x}) )} |D^2T(x)|_{\mathscr{B}(X^2, X)} &\le 2 \nu \left( |\pi^n A \pi^n|_{\mathscr{B}(X, X)} + \frac{1}{n+1} \right).
+\|T(\bar{x}) - \bar{x}\|_{\ell^1_{\nu, \mathbb{N}}} &\le \|\Pi_n A \Pi_n F(\bar{x})\|_{\ell^1_{\nu, \mathbb{N}}} + \frac{1}{n+1} \|\Pi_{\infty(n)} F(\bar{x})\|_{\ell^1_{\nu, \mathbb{N}}},\\
+\|DT(\bar{x})\|_{\mathscr{B}(\ell^1_{\nu, \mathbb{N}}, \ell^1_{\nu, \mathbb{N}})} &\le \max\left(\|\Pi_n A \Pi_n DF(\bar{x}) \Pi_n - \Pi_n\|_{\mathscr{B}(\ell^1_{\nu, \mathbb{N}}, \ell^1_{\nu, \mathbb{N}})}, \frac{\nu}{n+1} \|2\bar{x} - 1\|_{\ell^1_{\nu, \mathbb{N}}}\right),\\
+\sup_{x \in \text{cl}( B_R(\bar{x}) )} \|D^2T(x)\|_{\mathscr{B}(\ell^1_{\nu, \mathbb{N}}, \mathscr{B}(\ell^1_{\nu, \mathbb{N}}, \ell^1_{\nu, \mathbb{N}}))} &\le 2 \nu \left( \|\Pi_n A \Pi_n\|_{\mathscr{B}(\ell^1_{\nu, \mathbb{N}}, \ell^1_{\nu, \mathbb{N}})} + \frac{1}{n+1} \right).
\end{aligned}
```
@@ -112,31 +112,31 @@ In particular, from the latter estimate, we may freely choose ``R = \infty``.
The computer-assisted proof may be implemented as follows:
```@example logistic_ivp
-ν = interval(2.0)
+ν = interval(2)
X = ℓ¹(GeometricWeight(ν))
R = Inf
x̄_interval = interval.(x̄)
-F_interval = Sequence(Taylor(2n+1), similar(coefficients(x̄_interval), 2n+2))
+F_interval = zeros(eltype(x̄_interval), Taylor(2n+1))
F!(F_interval, x̄_interval)
tail_F_interval = copy(F_interval)
-tail_F_interval[0:n] .= interval(0.0)
+tail_F_interval[0:n] .= interval(0)
-DF_interval = LinearOperator(Taylor(n), Taylor(n), similar(coefficients(x̄_interval), n+1, n+1))
+DF_interval = zeros(eltype(x̄_interval), Taylor(n), Taylor(n))
DF!(DF_interval, x̄_interval)
-A = inv(mid.(DF_interval))
+A = interval.(inv(mid.(DF_interval)))
bound_tail_A = inv(interval(n+1))
# computation of the bounds
Y = norm(A * F_interval, X) + bound_tail_A * norm(tail_F_interval, X)
-Z₁ = opnorm(A * DF_interval - I, X) + bound_tail_A * ν * norm(2x̄_interval - 1, X)
+Z₁ = max(opnorm(A * DF_interval - UniformScaling(interval(1)), X), bound_tail_A * ν * norm(interval(2) * x̄_interval - interval(1), X))
-Z₂ = 2ν * (opnorm(interval.(A), X) + bound_tail_A)
+Z₂ = (opnorm(A, X) + bound_tail_A) * ν * interval(2)
setdisplay(:full)
diff --git a/docs/src/examples/infinite_dimensional_proofs/ode/lorenz_po.md b/docs/src/examples/infinite_dimensional_proofs/ode/lorenz_po.md
index 196cfd3a..cb3a7fc7 100644
--- a/docs/src/examples/infinite_dimensional_proofs/ode/lorenz_po.md
+++ b/docs/src/examples/infinite_dimensional_proofs/ode/lorenz_po.md
@@ -3,7 +3,7 @@
In this example, we will prove the existence of a periodic orbit of the Lorenz system
```math
-\frac{d}{dt} u(t) = f(u(t), \sigma, \rho, \beta) :=
+\frac{d}{dt} u(t) = f(u(t), \sigma, \rho, \beta) \bydef
\begin{pmatrix}
\sigma(u_2(t) - u_1(t))\\
u_1(t)(\rho - u_3(t)) - u_2(t)\\
@@ -11,7 +11,7 @@ u_1(t) u_2(t) - \beta u_3(t)
\end{pmatrix}, \qquad \sigma, \rho, \beta \in \mathbb{R}.
```
-The vector field ``f`` and its derivative with respect to ``u``, denoted ``Df``, may be implemented as follows:
+The vector field ``f`` and its derivative with respect to ``u``, denoted ``D_u f``, may be implemented as follows:
```@example lorenz_po
using RadiiPolynomial
@@ -43,27 +43,27 @@ nothing # hide
Let ``\nu > 1``,
```math
-X_\textnormal{F} := \left\{ u \in \mathbb{C}^\mathbb{Z} \, : \, |u|_{X_\textnormal{F}} := \sum_{k \in \mathbb{Z}} |u_k| \nu^{|k|} \right\}
+\ell^1_{\nu, \mathbb{Z}} \bydef \left\{ u \in \mathbb{C}^\mathbb{Z} \, : \, \|u\|_{\ell^1_{\nu, \mathbb{Z}}} \bydef \sum_{k \in \mathbb{Z}} |u_k| \nu^{|k|} \right\}
```
-and ``* : X_\textnormal{F} \times X_\textnormal{F} \to X_\textnormal{F}`` be the discrete convolution given by
+and ``* : \ell^1_{\nu, \mathbb{Z}} \times \ell^1_{\nu, \mathbb{Z}} \to \ell^1_{\nu, \mathbb{Z}}`` be the discrete convolution given by
```math
-u * v := \left\{ \sum_{l \in \mathbb{Z}} u_{k - l} v_l \right\}_{k \in \mathbb{Z}}, \qquad \text{for all } u, v \in X_\textnormal{F}.
+u * v \bydef \left\{ \sum_{l \in \mathbb{Z}} u_{k - l} v_l \right\}_{k \in \mathbb{Z}}, \qquad \text{for all } u, v \in \ell^1_{\nu, \mathbb{Z}}.
```
-For any sequence ``u \in X_\textnormal{F}``, the Fourier series ``\sum_{k \in \mathbb{Z}} u_k e^{i \omega k t}``, for some frequency ``\omega > 0``, defines an analytic ``2\pi\omega^{-1}``-periodic function in ``C^\omega(\mathbb{R}, \mathbb{C})``; while the discrete convolution ``*`` corresponds to the product of Fourier series in sequence space.
+For any sequence ``u \in \ell^1_{\nu, \mathbb{Z}}``, the Fourier series ``\sum_{k \in \mathbb{Z}} u_k e^{i \omega k t}``, for some frequency ``\omega > 0``, defines an analytic ``2\pi\omega^{-1}``-periodic function in ``C^\omega(\mathbb{R}, \mathbb{C})``; while the discrete convolution ``*`` corresponds to the product of Fourier series in sequence space.
-The Banach space ``X_\textnormal{F}`` is a suitable space to represent each component of a periodic solution of the Lorenz system. Indeed, it is a standard result from ODE theory that analytic vector fields yield analytic solutions.[^1]
+The Banach space ``\ell^1_{\nu, \mathbb{Z}}`` is a suitable space to represent each component of a periodic solution of the Lorenz system. Indeed, it is a standard result from ODE theory that analytic vector fields yield analytic solutions.[^1]
[^1]: A. Hungria, J.-P. Lessard and J. D. Mireles James, [Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach](https://doi.org/10.1090/mcom/3046), *Mathematics of Computation*, **85** (2016), 1427-1459.
-Define the Banach space ``X := \mathbb{C} \times X_\textnormal{F}^3`` endowed with the norm ``|x|_X := \max(|\gamma|, |u_1|_{X_\textnormal{F}}, |u_2|_{X_\textnormal{F}}, |u_3|_{X_\textnormal{F}})`` for all ``x = (\gamma, u_1, u_2, u_3) \in X``. It follows that the sequence of coefficients of a ``2\pi\gamma``-periodic Fourier series solving the Lorenz equations is a zero of the mapping ``F : X \to X`` given by
+Define the Banach space ``X \bydef \mathbb{C} \times (\ell^1_{\nu, \mathbb{Z}})^3`` endowed with the norm ``\|x\|_X \bydef |\gamma| + \|u_1\|_{\ell^1_{\nu, \mathbb{Z}}} + \|u_2\|_{\ell^1_{\nu, \mathbb{Z}}} + \|u_3\|_{\ell^1_{\nu, \mathbb{Z}}}`` for all ``x = (\gamma, u_1, u_2, u_3) \in X``. It follows that the sequence of coefficients of a ``2\pi\gamma``-periodic Fourier series solving the Lorenz equations is a zero of the mapping ``F : X \to X`` given by
```math
-F(x) :=
+F(x) \bydef
\begin{pmatrix}
-\sum_{j = 1}^3 (\sum_{k = -n}^n (u_j)_k - \xi_j)\eta_j\\
+\sum_{j = 1}^3 (\sum_{k = -K}^K (u_j)_k - \xi_j)\eta_j\\
\left\{ \gamma ( f(u, \sigma, \rho, \beta) )_k - i k u_k \right\}_{k \in \mathbb{Z}}
\end{pmatrix}, \qquad \text{for all } x = (\gamma, u_1, u_2, u_3) \in \text{domain}(F),
```
@@ -73,27 +73,27 @@ where ``\xi \in \mathbb{R}^3`` is a chosen approximate position of the periodic
The mapping ``F`` and its Fréchet derivative, denoted ``DF``, may be implemented as follows:
```@example lorenz_po
-function F!(F, x, σ, ρ, β)
+function F!(F, x, σ, ρ, β, ξ, η)
γ, u = x[1], component(x, 2)
F[1] =
- (sum(component(u, 1)) - 10.205222700615433) * 24.600655549587863 +
- (sum(component(u, 2)) - 11.899530531689562) * (-2.4927169722923335) +
- (sum(component(u, 3)) - 27.000586375896557) * 71.81142025024573
+ (component(u, 1)(0) - ξ[1]) * η[1] +
+ (component(u, 2)(0) - ξ[2]) * η[2] +
+ (component(u, 3)(0) - ξ[3]) * η[3]
project!(component(F, 2), γ * f!(component(F, 2), u, σ, ρ, β) - differentiate(u))
return F
end
-function DF!(DF, x, σ, ρ, β)
+function DF!(DF, x, σ, ρ, β, η)
γ, u = x[1], component(x, 2)
DF .= 0
- component(component(DF, 1, 2), 1)[1,:] .= 24.600655549587863
- component(component(DF, 1, 2), 2)[1,:] .= -2.4927169722923335
- component(component(DF, 1, 2), 3)[1,:] .= 71.81142025024573
+ component(component(DF, 1, 2), 1)[1,:] .= η[1]
+ component(component(DF, 1, 2), 2)[1,:] .= η[2]
+ component(component(DF, 1, 2), 3)[1,:] .= η[3]
f!(component(DF, 2, 1), u, σ, ρ, β)
@@ -107,7 +107,7 @@ nothing # hide
Consider the fixed-point operator ``T : X \to X`` defined by
```math
-T(x) := x - A F(x),
+T(x) \bydef x - A F(x),
```
where ``A : X \to X`` is an injective operator corresponding to an approximation of ``DF(\bar{x})^{-1}`` for some numerical zero ``\bar{x} = (\bar{\gamma}, \bar{u}_1, \bar{u}_2, \bar{u}_3) \in X`` of ``F``.
@@ -117,109 +117,116 @@ Given an initial guess, the numerical zero ``\bar{x}`` of ``F`` may be obtained
```@example lorenz_po
σ, ρ, β = 10.0, 28.0, 8/3
-n = 400
-
-x̄ = Sequence(ParameterSpace() × Fourier(n, 1.0)^3, zeros(ComplexF64, 1+3*(2n+1)))
-x̄[1] = 9.150971830259179/2π # γ, i.e. approximate inverse of the frequency
-component(component(x̄, 2), 1)[0:14] =
- [6.25, -0.66 - 1.45im, 0.6 - 1.2im, 1.11 - 0.26im, 0.77 + 0.57im,
- 0.08 + 0.76im, -0.35 + 0.45im, -0.39 + 0.13im, -0.37 - 0.0008im, -0.44 - 0.23im,
- -0.18 - 0.68im, 0.65 - 0.61im, 0.80 + 0.50im, -0.53 + 0.43im, 1.25 - 0.07im]
-component(component(x̄, 2), 2)[0:14] =
- [6.25, -0.56 - 1.5im, 0.76 - 1.12im, 1.17 - 0.03im, 0.62 + 0.78im,
- -0.18 + 0.76im,-0.54 + 0.3im, -0.45 - 0.06im, -0.37 - 0.2im, -0.3 - 0.51im,
- 0.29 - 0.8im, 1.11 - 0.13im, 0.4 + 1.16im, -0.91 - 0.05im, 1.31 + 1.13im]
-component(component(x̄, 2), 3)[0:14] =
- [24.45, -0.22 - 1.62im, 1.13 - 0.83im, 1.2 + 0.53im, 0.14 + 1.28im,
- -1.03 + 0.75im, -1.14 - 0.52im, -0.08 - 1.21im, 0.98 - 0.57im, 0.79 + 0.59im,
- -0.27 + 0.69im, -0.34 - 0.23im, 0.57 + 0.22im, -1.23 + 1.02im, 0.75 - 2.69im]
-component(component(x̄, 2), 1)[-14:-1] .= conj.(component(component(x̄, 2), 1)[14:-1:1])
-component(component(x̄, 2), 2)[-14:-1] .= conj.(component(component(x̄, 2), 2)[14:-1:1])
-component(component(x̄, 2), 3)[-14:-1] .= conj.(component(component(x̄, 2), 3)[14:-1:1])
-
-newton!((F, DF, x) -> (F!(F, x, σ, ρ, β), DF!(DF, x, σ, ρ, β)), x̄)
+K = 60
+
+x̄ = zeros(ComplexF64, ParameterSpace() × Fourier(K, 1.0)^3)
+x̄[1] = 1.5/(2π) # γ, i.e. approximate inverse of the frequency
+component(component(x̄, 2), 1)[1:2:5] =
+ [-2.9 - 4.3im,
+ 1.6 - 1.1im,
+ 0.3 + 0.4im]
+component(component(x̄, 2), 2)[1:2:5] =
+ [-1.2 - 5.4im,
+ 3.0 + 0.8im,
+ -0.4 + 1.1im]
+component(component(x̄, 2), 3)[0:2:4] =
+ [ 23,
+ 3.8 + 4.7im,
+ -1.8 + 0.9im]
+component(component(x̄, 2), 1)[-5:2:-1] .= conj.(component(component(x̄, 2), 1)[5:-2:1])
+component(component(x̄, 2), 2)[-5:2:-1] .= conj.(component(component(x̄, 2), 2)[5:-2:1])
+component(component(x̄, 2), 3)[-4:2:0] .= conj.(component(component(x̄, 2), 3)[4:-2:0])
+
+ξ = component(x̄, 2)(0)
+η = differentiate(component(x̄, 2))(0)
+
+newton!((F, DF, x) -> (F!(F, x, σ, ρ, β, ξ, η), DF!(DF, x, σ, ρ, β, η)), x̄)
# impose that x̄[1] is real and component(x̄, 2) are the coefficients of a real Fourier series
x̄[1] = real(x̄[1])
for i ∈ 1:3
component(component(x̄, 2), i)[0] = real(component(component(x̄, 2), i)[0])
- component(component(x̄, 2), i)[-n:-1] .= conj.(component(component(x̄, 2), i)[n:-1:1])
+ component(component(x̄, 2), i)[-K:-1] .= conj.(component(component(x̄, 2), i)[K:-1:1])
end
```
-Let ``R > 0``. Since ``T \in C^2(X, X)`` we may use the [second-order Radii Polynomial Theorem](@ref second_order_RPT) such that we need to estimate ``|T(\bar{x}) - \bar{x}|_X``, ``|DT(\bar{x})|_{\mathscr{B}(X, X)}`` and ``\sup_{x \in \text{cl}( B_R(\bar{x}) )} |D^2T(x)|_{\mathscr{B}(X^2, X)}``.
+Let ``R > 0``. Since ``T \in C^2(X, X)`` we may use the [second-order Radii Polynomial Theorem](@ref second_order_RPT) such that we need to estimate ``\|T(\bar{x}) - \bar{x}\|_X``, ``\|DT(\bar{x})\|_{\mathscr{B}(X, X)}`` and ``\sup_{x \in \text{cl}( B_R(\bar{x}) )} \|D^2T(x)\|_{\mathscr{B}(X, \mathscr{B}(X, X))}``.
To this end, consider the truncation operator
```math
-(\pi^n u)_k :=
+(\Pi_K u)_k \bydef
\begin{cases}
-u_k, & |k| \le n,\\
-0, & |k| > n,
+u_k, & |k| \le K,\\
+0, & |k| > K,
\end{cases}
-\qquad \textnormal{for all } u \in X_\textnormal{F}.
+\qquad \text{for all } u \in \ell^1_{\nu, \mathbb{Z}}.
```
-Using the same symbol, this projection extends naturally to ``X_\textnormal{F}^3`` and ``X`` by acting on each component as follows ``\pi^n u := (\pi^n u_1, \pi^n u_2, \pi^n u_3)``, for all ``u = (u_1, u_2, u_3) \in X_\textnormal{F}^3``, and ``\pi^n x := (\gamma, \pi^n u_1, \pi^n u_2, \pi^n u_3)``, for all ``x = (\gamma, u_1, u_2, u_3) \in X``. For each of the Banach spaces ``X_\textnormal{F}, X_\textnormal{F}^3, X``, we define the complementary operator ``\pi^{\infty(n)} := I - \pi^n``.
+Using the same symbol, this projection extends naturally to ``(\ell^1_{\nu, \mathbb{Z}})^3`` and ``X`` by acting on each component as follows ``\Pi_K u \bydef (\Pi_K u_1, \Pi_K u_2, \Pi_K u_3)``, for all ``u = (u_1, u_2, u_3) \in (\ell^1_{\nu, \mathbb{Z}})^3``, and ``\Pi_K x \bydef (\gamma, \Pi_K u_1, \Pi_K u_2, \Pi_K u_3)``, for all ``x = (\gamma, u_1, u_2, u_3) \in X``. For each of the Banach spaces ``\ell^1_{\nu, \mathbb{Z}}, (\ell^1_{\nu, \mathbb{Z}})^3, X``, we define the complementary operator ``\Pi_{\infty(K)} \bydef I - \Pi_K``.
Thus, denoting ``\bar{u} = (\bar{u}_1, \bar{u}_2, \bar{u}_3)``, we have
```math
\begin{aligned}
-|T(\bar{x}) - \bar{x}|_X &\le
-|\pi^n A \pi^n F(\bar{x})|_X + \frac{\bar{\gamma}}{n+1} | \pi^{\infty(n)} f(\bar{u}) |_{X_\textnormal{F}^3},\\
-|DT(\bar{x})|_{\mathscr{B}(X, X)} &\le
-|\pi^n A \pi^n DF(\bar{x}) \pi^{2n} - \pi^n|_{\mathscr{B}(X, X)} + \frac{1}{n+1} | \pi^{\infty(n)} f(\bar{u}) |_{X_\textnormal{F}^3} +\\
-&\qquad \frac{\bar{\gamma}}{n+1}
-\max \left(2 \sigma, 1 + |\bar{u}_1|_{X_\textnormal{F}} + |\rho - \bar{u}_3|_{X_\textnormal{F}}, \beta + |\bar{u}_1|_{X_\textnormal{F}} + |\bar{u}_2|_{X_\textnormal{F}}\right),\\
-\sup_{x \in \text{cl}( B_R(\bar{x}) )} |D^2T(x)|_{\mathscr{B}(X^2, X)} &\le
-2\left(|\pi^n A \pi^n|_{\mathscr{B}(X, X)} + \frac{1}{n+1}\right) \Big(\bar{\gamma} + R +\\
-&\qquad \max \left(2 \sigma, 1 + \rho + |\bar{u}_1|_{X_\textnormal{F}} + |\bar{u}_3|_{X_\textnormal{F}} + 2R, \beta + |\bar{u}_1|_{X_\textnormal{F}} + |\bar{u}_2|_{X_\textnormal{F}} + 2R\right)\Big).
+\|T(\bar{x}) - \bar{x}\|_X &\le
+\|\Pi_K A \Pi_K F(\bar{x})\|_X + \frac{\bar{\gamma}}{n+1} \|\Pi_{\infty(K)} f(\bar{u}, \sigma, \rho, \beta)\|_{(\ell^1_{\nu, \mathbb{Z}})^3},\\
+\|DT(\bar{x})\|_{\mathscr{B}(X, X)} &\le
+\|\Pi_K A \Pi_K DF(\bar{x}) \Pi_{2K} - \Pi_K\|_{\mathscr{B}(X, X)} + \frac{1}{n+1} \max\Big( \|\Pi_{\infty(K)} f(\bar{u}, \sigma, \rho, \beta)\|_{(\ell^1_{\nu, \mathbb{Z}})^3},\\
+&\qquad \bar{\gamma} \max\left(\sigma + \|\rho-\bar{u}_3\|_{\ell^1_{\nu, \mathbb{Z}}} + \|\bar{u}_2\|_{\ell^1_{\nu, \mathbb{Z}}}, \sigma + 1 + \|\bar{u}_1\|_{\ell^1_{\nu, \mathbb{Z}}}, \|\bar{u}_1\|_{\ell^1_{\nu, \mathbb{Z}}} + \beta\right) \Big),\\
+\sup_{x \in \text{cl}( B_R(\bar{x}) )} \|D^2T(x)\|_{\mathscr{B}(X, \mathscr{B}(X, X))} &\le
+\left(\|\Pi_K A \Pi_K\|_{\mathscr{B}(X, X)} + \frac{1}{n+1}\right) \max\Big( 2 (\bar{\gamma} + R),\\
+&\qquad \max\left(\sigma + \|\rho-\bar{u}_3\|_{\ell^1_{\nu, \mathbb{Z}}} + \|\bar{u}_2\|_{\ell^1_{\nu, \mathbb{Z}}} + 2R, \sigma + 1 + \|\bar{u}_1\|_{\ell^1_{\nu, \mathbb{Z}}} + R, \|\bar{u}_1\|_{\ell^1_{\nu, \mathbb{Z}}} + R + \beta\right) \Big).
\end{aligned}
```
The computer-assisted proof may be implemented as follows:
```@example lorenz_po
-ν = interval(1.01)
+ν = interval(1.05)
X_F = ℓ¹(GeometricWeight(ν))
-X_F³ = NormedCartesianSpace(X_F, ℓ∞())
-X = NormedCartesianSpace((ℓ∞(), X_F³), ℓ∞())
-R = 1e-8
+X_F³ = NormedCartesianSpace(X_F, ℓ¹())
+X = NormedCartesianSpace((ℓ¹(), X_F³), ℓ¹())
+R = 1e-10
-σ_interval, ρ_interval, β_interval = interval(10.0), interval(28.0), 8.0/interval(3.0)
+σ_interval, ρ_interval, β_interval = interval(10), interval(28), interval(8)/interval(3)
-x̄_interval = Sequence(ParameterSpace() × Fourier(n, interval(1.0))^3, interval.(coefficients(x̄)))
+x̄_interval = Sequence(ParameterSpace() × Fourier(K, interval(1))^3, interval.(coefficients(x̄)))
γ̄_interval = real(x̄_interval[1])
ū_interval = component(x̄_interval, 2)
-F_interval = Sequence(ParameterSpace() × Fourier(2n, interval(1.0))^3, similar(coefficients(x̄_interval), 1+3*(4n+1)))
-F!(F_interval, x̄_interval, σ_interval, ρ_interval, β_interval)
+ξ_interval = interval.(ξ)
+η_interval = interval.(η)
+
+F_interval = zeros(eltype(x̄_interval), ParameterSpace() × Fourier(2K, interval(1))^3)
+F!(F_interval, x̄_interval, σ_interval, ρ_interval, β_interval, ξ_interval, η_interval)
tail_γ̄f_interval = copy(component(F_interval, 2))
for i ∈ 1:3
- component(tail_γ̄f_interval, i)[-n:n] .= 0
+ component(tail_γ̄f_interval, i)[-K:K] .= interval(0)
end
-DF_interval = LinearOperator(space(F_interval), space(x̄_interval), similar(coefficients(x̄_interval), length(x̄_interval), length(F_interval)))
-DF!(DF_interval, x̄_interval, σ_interval, ρ_interval, β_interval)
+DF_interval = zeros(eltype(x̄_interval), space(F_interval), space(x̄_interval))
+DF!(DF_interval, x̄_interval, σ_interval, ρ_interval, β_interval, η_interval)
-A = inv(mid.(project(DF_interval, space(x̄_interval), space(x̄_interval))))
-bound_tail_A = inv(interval(n+1))
+A = interval.(inv(mid.(project(DF_interval, space(x̄_interval), space(x̄_interval)))))
+bound_tail_A = inv(interval(K+1))
# computation of the bounds
Y = norm(A * F_interval, X) + bound_tail_A * norm(tail_γ̄f_interval, X_F³)
-Z₁ = opnorm(A * DF_interval - I, X) + bound_tail_A * norm(tail_γ̄f_interval, X_F³) / γ̄_interval +
- bound_tail_A * γ̄_interval * max(2σ_interval,
- 1 + norm(component(ū_interval, 1), X_F) + norm(ρ_interval-component(ū_interval, 3), X_F),
- β_interval + norm(component(ū_interval, 1), X_F) + norm(component(ū_interval, 2), X_F))
+opnorm_Df = max(σ_interval + norm(ρ_interval-component(ū_interval, 3), X_F) + norm(component(ū_interval, 2), X_F),
+ σ_interval + 1 + norm(component(ū_interval, 1), X_F),
+ norm(component(ū_interval, 1), X_F) + β_interval)
+
+Z₁ = opnorm(A * DF_interval - UniformScaling(interval(1)), X) +
+ bound_tail_A * max(norm(tail_γ̄f_interval / γ̄_interval, X_F³), γ̄_interval * opnorm_Df)
-Z₂ = (opnorm(interval.(A), X) + bound_tail_A) * 2 * (γ̄_interval + R +
- max(2σ_interval,
- 1 + ρ_interval + norm(component(ū_interval, 1), X_F) + norm(component(ū_interval, 3), X_F) + 2R,
- β_interval + norm(component(ū_interval, 1), X_F) + norm(component(ū_interval, 2), X_F) + 2R))
+Z₂ = (opnorm(A, X) + bound_tail_A) * max(2 * (γ̄_interval + R),
+ max(σ_interval + norm(ρ_interval - component(ū_interval, 3), X_F) + R + norm(component(ū_interval, 2), X_F) + R,
+ σ_interval + 1 + norm(component(ū_interval, 1), X_F) + R,
+ norm(component(ū_interval, 1), X_F) + R + β_interval))
setdisplay(:full)
diff --git a/docs/src/examples/infinite_dimensional_proofs/ode/lorenz_po.mp4 b/docs/src/examples/infinite_dimensional_proofs/ode/lorenz_po.mp4
index 08d093ff..0291d7cc 100644
Binary files a/docs/src/examples/infinite_dimensional_proofs/ode/lorenz_po.mp4 and b/docs/src/examples/infinite_dimensional_proofs/ode/lorenz_po.mp4 differ
diff --git a/docs/src/radii_polynomial_approach.md b/docs/src/radii_polynomial_approach.md
index 004c8daf..7a123e19 100644
--- a/docs/src/radii_polynomial_approach.md
+++ b/docs/src/radii_polynomial_approach.md
@@ -31,8 +31,8 @@ Let ``X`` be a Banach space, ``U`` an open subset of ``X``, ``T \in C^1(U, X)``
- (First-order) Suppose ``Y, Z_1 \ge 0`` satisfy
```math
\begin{aligned}
-|T(\bar{x}) - \bar{x}|_X &\le Y,\\
-\sup_{x \in \text{cl}( B_R(\bar{x}) )} |DT(x)|_{\mathscr{B}(X, X)} &\le Z_1,
+\|T(\bar{x}) - \bar{x}\|_X &\le Y,\\
+\sup_{x \in \text{cl}( B_R(\bar{x}) )} \|DT(x)\|_{\mathscr{B}(X, X)} &\le Z_1,
\end{aligned}
```
and define the *radii polynomial* by ``p(r) \bydef Y + (Z_1 - 1) r``.
@@ -40,9 +40,9 @@ If there exists a *radius* ``\bar{r} \in [0, R]`` such that ``p(\bar{r}) \le 0``
- (Second-order) Suppose ``Y, Z_1, Z_2 \ge 0`` satisfy
```math
\begin{aligned}
-|T(\bar{x}) - \bar{x}|_X &\le Y,\\
-|DT(\bar{x})|_{\mathscr{B}(X, X)} &\le Z_1,\\
-|DT(x) - DT(\bar{x})|_{\mathscr{B}(X, X)} &\le Z_2 |x - \bar{x}|_X, \qquad \text{for all } x \in \text{cl}( B_R(\bar{x}) ),
+\|T(\bar{x}) - \bar{x}\|_X &\le Y,\\
+\|DT(\bar{x})\|_{\mathscr{B}(X, X)} &\le Z_1,\\
+\|DT(x) - DT(\bar{x})\|_{\mathscr{B}(X, X)} &\le Z_2 \|x - \bar{x}\|_X, \qquad \text{for all } x \in \text{cl}( B_R(\bar{x}) ),
\end{aligned}
```
and define the *radii polynomial* by ``p(r) \bydef Y + (Z_1 - 1) r + \frac{Z_2}{2} r^2``.
diff --git a/docs/src/sequence_spaces/sequences.md b/docs/src/sequence_spaces/sequences.md
index 4b6f3516..a3a9e972 100644
--- a/docs/src/sequence_spaces/sequences.md
+++ b/docs/src/sequence_spaces/sequences.md
@@ -73,7 +73,7 @@ To improve performance, the FFT algorithm may be used to compute discrete convol
```@repl sequences
x = Sequence(Taylor(3), interval.([inv(10_000.0 ^ i) for i ∈ 0:3]))
x³ = x ^ 3
-x³_fft = rifft!(similar(x³), fft(x, fft_size(space(x), 3)) .^ 3)
+x³_fft = rifft!(zero(x³), fft(x, fft_size(space(x), 3)) .^ 3)
```
To circumvent machine precision limitations, the `banach_rounding!` method enclose rigorously each term of the convolution beyond a prescribed order.[^1]