• Parameter continuation for the cubic root
• • • Step 1: Problem definition
    • Step 2: Approximate zero (floating-point arithmetic)
    • Step 3: Approximate inverse (floating-point arithmetic)
    • Step 4: Bounds estimation (interval arithmetic)

Parameter continuation for the cubic root

In this example, we prove the existence of a one-parameter family of solutions to

\[0 = f(u, \lambda) \bydef u^3 - \lambda - 2, \qquad \lambda \in [-1, 1].\]

Step 1: Problem definition

The map $f$ and its derivative $D_u f$ are implemented as follows:

f(u, λ) = u^3 - λ - exact(2)

Duf(u, λ) = exact(3) * u^2

Consider the fixed-point operator $T : \mathbb{R} \times [-1,1] \to \mathbb{R}$ given by

\[T(u, \lambda) \bydef u - A(\lambda) f(u, \lambda),\]

where $A(\lambda) : \mathbb{R} \to \mathbb{R}$ is an injective operator corresponding to a numerical approximation of $D_u f(\bar{u}(\lambda), \lambda)^{-1}$, for some approximate zero $\bar{u}(\lambda) \in \mathbb{R}$ of $f$, for all $\lambda \in [-1, 1]$.

Step 2: Approximate zero (floating-point arithmetic)

We use the numerical continuation method to retrieve a numerical approximation of the curve.

We construct a grid of parameters and iterate Newton's method for each step, using the previous approximate zero as the predictor of the solution at the next step.

using RadiiPolynomial

K = 10
K_fft = fft_size(Chebyshev(K))
npts = K_fft ÷ 2 + 1

λ_grid = [-cospi(2j/K_fft) for j = 0:npts-1]
u_grid = Vector{Float64}(undef, npts)

# initialize

u_init = 1.0
u_init, success_newton = newton(u -> (f(u, λ_grid[1]), Duf(u, λ_grid[1])), u_init)

u_grid[1] = u_init

# run continuation scheme

for j = 2:npts
    w = u_grid[j-1] # predictor

    u_bar, success_newton_j = newton(u -> (f(u, λ_grid[j]), Duf(u, λ_grid[j])), w; verbose = true)
    success_newton_j || error()

    u_grid[j] = u_bar
end

# construct the approximation

u_fft = [reverse(u_grid) ; u_grid[begin+1:end-1]]
u_cheb = real(to_seq(u_fft, Chebyshev(K)))

Sequence in Chebyshev(10) with coefficients Vector{Float64}:
  1.2410369149251503
  0.10907140054331423
 -0.009688555942083859
  0.0014382541523247298
 -0.0002564899825927712
  5.0342871589922137e-5
 -1.0483288779396366e-5
  2.2725722957282205e-6
 -5.072065108469026e-7
  1.1572783048626124e-7
 -2.6866343460769933e-8

Step 3: Approximate inverse (floating-point arithmetic)

We construct the approximate inverse $A(\lambda) \approx D_u f(\bar{u}(\lambda), \lambda)^{-1}$ across the continuation branch using standard floating-point arithmetic.

A_grid = inv.(Duf.(u_grid, λ_grid))
A_fft = [reverse(A_grid) ; A_grid[begin+1:end-1]]
A_cheb = real(to_seq(A_fft, Chebyshev(K)))

Sequence in Chebyshev(10) with coefficients Vector{Float64}:
  0.22730980307522403
 -0.040940634508741705
  0.009167001988608968
 -0.002186410740412038
  0.0005374770746072676
 -0.00013449087973091016
  3.404835870885278e-5
 -8.691414303604028e-6
  2.2323466863233277e-6
 -5.761101173660076e-7
  1.4924567912798963e-7

Step 4: Bounds estimation (interval arithmetic)

To apply the Radii Polynomial Theorem, we need to theoretically derive and rigorously evaluate the bounds $Y, Z_1, Z_2$. The computer-assisted proof is completed by evaluating these bounds with interval arithmetic:

λ_cheb_interval = interval(Sequence(Chebyshev(1), [0, 0.5]))
u_cheb_interval = interval(u_cheb)
A_cheb_interval = interval(A_cheb)

#- Y bound

Y = norm(A_cheb_interval * f(u_cheb_interval, λ_cheb_interval), 1)

#- Z₁ bound

Z₁ = norm(exact(1) - A_cheb_interval * Duf(u_cheb_interval, λ_cheb_interval), 1)

#- Z₂ bound

R = 10 * sup(Y)
Z₂ = exact(3) * norm(A_cheb_interval, 1) * (norm(u_cheb_interval, 1) + exact(2 * R))

# verify the contraction

ie, contraction_success = interval_of_existence(Y, Z₁, Z₂, R; verbose = true)

┌ Info: success: interval found
│ Y = [1.66179e-8, 1.66179e-8]_com
│ Z₁ = [4.19304e-7, 4.19304e-7]_com
│ Z₂ = [1.48207, 1.48207]_com
│ R = 1.6617944662661552e-7
│ Δ = [0.999999, 0.999999]_com
└ roots = ([1.66179e-8, 1.6618e-8]_com, [1.34946, 1.34946]_com)

inf(ie) # smallest error

1.6617954922205423e-8

The following figure^[1] shows the numerical approximation of the proven branch of the cubic root.

using GLMakie

fig = Figure()
ax = Axis(fig[1,1], xticks = -3:3)
lines!(ax, [Point2f(λ, cbrt(λ + 2)) for λ = LinRange(-3, 2, 501)];
    color = :black, label = L"(λ+2)^{1/3}")
lines!(ax, [Point2f(λ, u_cheb(λ)) for λ = LinRange(-1, 1, 501)];
    color = :blue, label = L"\bar{u}(λ)")
scatter!(ax, [Point2f(λ_grid[j], u_grid[j]) for j = 1:npts];
    color = :red)
axislegend(ax; position = :lt)
fig