• Equilibria of the Lorenz system
• • • Step 1: Formulation
    • Step 2: Approximation (floating-point arithmetic)
    • • The approximate zero
      • The approximate inverse
    • Step 3: Bounds (interval arithmetic)
    • • First-order validation
      • The second-order bound
    • Step 4: Conclusion

Equilibria of the Lorenz system

We prove the existence of the equilibria of the Lorenz system

\[\begin{aligned} \dot{u}^{(1)} &= \sigma (u^{(2)} - u^{(1)}), \\ \dot{u}^{(2)} &= u^{(1)}(\rho - u^{(3)}) - u^{(2)}, \\ \dot{u}^{(3)} &= u^{(1)} u^{(2)} - \beta u^{(3)}, \end{aligned} \qquad (\sigma, \rho, \beta) = (10, 28, 8/3).\]

The unknowns are a point of $\mathbb{R}^3$ rather than a number, which is the simplest possible setting in which to meet the library's cartesian spaces.

Step 1: Formulation

The equilibria are the zeros of the vector field itself, $F(u) = 0$ with $u = (u^{(1)}, u^{(2)}, u^{(3)}) \in \mathbb{R}^3$. We take $\mathscr{X} = \mathbb{R}^3$ with the norm $\|u\|_{\mathscr{X}} = \max\left(|u^{(1)}|, |u^{(2)}|, |u^{(3)}| \right)$.

In the library, $\mathbb{R}^3$ is constructed via ScalarSpace()^3: a CartesianPower of the space of plain numbers ScalarSpace. A Sequence in that space just wraps a 3-coefficient vector and a LinearOperator acting on them wraps a 3-by-3 matrix.

using RadiiPolynomial, LinearAlgebrafunction F(u, params)    σ, ρ, β = params    u1, u2, u3 = u[1], u[2], u[3]    return Sequence(space(u), [σ*(u2 - u1), u1*(ρ - u3) - u2, u1*u2 - β*u3])endfunction DF(u, params)    σ, ρ, β = params    u1, u2, u3 = u[1], u[2], u[3]    return LinearOperator(space(u), space(u),        [        -σ            σ       zero(u1)             ρ - u3      -one(u1)            -u1                 u2            u1            -β ])end

The norm is prescribed separately, and the maximum of the components absolute value is embodied by $\ell^\infty$.

X = ScalarSpace()^3X_norm = EllInf()

Step 2: Approximation (floating-point arithmetic)

The approximate zero

params = (10.0, 28.0, 8/3)u_bar, converged = newton(u -> (F(u, params), DF(u, params)), Sequence(X, [8.0, 8.0, 26.0]))

(Sequence(ScalarSpace()³, [8.48528137423857, 8.48528137423857, 27.0]), true)

The components are reachable with component, which is how to retrieve a block of a cartesian unknown:

component(u_bar, 3)

Sequence in 𝕂 with coefficients SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}:
 27.0

The approximate inverse

$DF(\num)$ is a $3 \times 3$ matrix, so A is simply an approximate matrix inverse.

A = inv(DF(u_bar, params))

Step 3: Bounds (interval arithmetic)

$\beta = 8/3$ is not a representable floating-point number, so the parameters must be enclosed rather than converted. interval(8)/interval(3) produces the tightest interval containing $8/3$; writing interval(8/3) instead would enclose the double nearest $8/3$, which is a different real number.

params_i = (interval(10), interval(28), interval(8)/interval(3))u_i = interval(u_bar)A_i = interval(A)Y = norm(A_i * F(u_i, params_i), X_norm)

[0.0, 1.67477e-15]_com

First-order validation

The first-order Radii Polynomial Theorem needs a bound on $\|I - A DF(u)\|$ over an entire ball. As in A first validation, interval arithmetic supplies it by making every component an interval of radius $R$ and evaluate DF there.

R = 10sup(Y)u_R = Sequence(X, interval.(coefficients(u_bar), R; format = :midpoint))Z₁_R = opnorm(I - A_i * DF(u_R, params_i), X_norm)ie_first, proved_first = interval_of_existence(Y, Z₁_R, R)inf(ie_first), sup(ie_first), proved_first

(1.6747652892627036e-15, 1.6747652892626938e-14, true)

That is already a complete proof, and it needed only $F$ and $DF$. But the chosen value of $R$, which bounds the uniqueness radius, must be small for $A$ to accuratly approximate uniformly $DF(u)^{-1}$ for all $u \in B(\num, R)$. The first-order theorem is not ideal to certify a large uniqueness radius.

The second-order bound

To do better, we put the limiting constraint $Z_1 < 1$ at the center of the ball $\|I - A DF(\num)\| < 1$ and we control how fast $DF$ varies nearby.

Z₁ = opnorm(I - A_i * DF(u_i, params_i), X_norm)

[0.0, 6.66134e-16]_com_NG

Writing $h = u - \num$, every entry of $DF$ is affine in $u$, so the difference is exactly

\[DF(u) - DF(\num) = \begin{pmatrix} 0 & 0 & 0\\ -h^{(3)} & 0 & -h^{(1)}\\ h^{(2)} & h^{(1)} & 0 \end{pmatrix}.\]

The $\ell^\infty$-induced norm of a matrix is its largest row sum, here $\max\left(0, |h^{(3)}| + |h^{(1)}|, |h^{(2)}| + |h^{(1)}|\right) \le 2 \|h\|_X$. Hence

\[\|A (DF(u) - DF(\num))\|_{\mathscr{L}(X,X)} \le 2 \|A\|_{\mathscr{L}(X,X)} \|u - \num\|_X, \qquad Z_2 \bydef 2 \|A\|_{\mathscr{L}(X,X)}.\]

This $Z_2$ is independent of $R$. That is to be expected since $F$ is quadratic, so $DF$ is affine and $DF(u) - DF(\num)$ is linear in $u - \num$. Thus, its Lipschitz constant is global. We may therefore take $R = \infty$.

R  = InfZ₂ = exact(2) * opnorm(A_i, X_norm)ie, proved = interval_of_existence(Y, Z₁, Z₂, R)inf(ie), sup(ie), proved

(3.399208571853739e-15, 2.18695349393644, true)

Step 4: Conclusion

There exists an exact equilibrium of the Lorenz system within inf(ie) of u_bar, and it is the only one within sup(ie).

The next table compares the two validations.