Experiment 2: Intersecting curves on a Hirzebruch surface¶
In this notebook we use the Cox homotopy to intersect curves on the Hirzebruch surface \(\mathscr{H}_2\).
We first load our routines.
[2]:
include("CoxHomotopy.jl")
[2]:
trackPathRandom (generic function with 1 method)
We generate two bivariate equations.
[21]:
n = 2;
@polyvar t[1:n];
f̂ = [1 + t[1] + t[2] + t[1]*t[2] + t[1]^2*t[2] + t[1]^3*t[2]; 1 + t[1] + 2*t[2] + 3*t[1]*t[2] + 4*t[1]^2*t[2] + 1*t[1]^3*t[2]];
f̂ = [1 + t[1] + t[2]*(t[1]-1)*(1+sum(map(i->randn()*t[1]^i,1:2))); 1 + t[1] + t[2]*(t[1]-1)*(1+sum(map(i->randn()*t[1]^i,1:2)))]
[21]:
2-element Array{DynamicPolynomials.Polynomial{true,Float64},1}:
-0.36515775180019705t₁³t₂ + 0.3626534954403268t₁²t₂ + 1.0025042563598703t₁t₂ + t₁ - t₂ + 1.0
0.4374666674815006t₁³t₂ + 0.0017765144844340552t₁²t₂ + 0.5607568180340654t₁t₂ + t₁ - t₂ + 1.0
We compute the polytope information.
[22]:
P, Fᵀ, a = computePolytope(f̂);
We solve the system using the Cox homotopy.
[23]:
@time homSols, f, toricSols, Fᵀ, x̂, K, B, data = coxHomotopy(f̂; fixSlice = true, details = true);
0.053850 seconds (67.02 k allocations: 4.872 MiB)
We compute the residuals of the toric solutions and the solutions in their homogeneous coordinates.
[24]:
display(get_residual(f̂,toricSols,t))
display(get_residual(f,homSols,variables(f)))
4-element Array{Float64,1}:
2.1506011421870678e-16
1.0816808319406293e-15
2.5502337173693015e-16
1.6503625390396785e-16
4-element Array{Float64,1}:
2.2484317241831537e-16
1.7830125319256865e-15
9.223413031523224e-17
2.08199467410775e-16
We plot the magnitude of the homogeneous coordinates of the solutions. Each solution corresponds to a column in the table and each row represents a Cox coordinate. Dark (blue) colors indicate small absolute values. The plot shows that only one solution is in the torus. The others are on the divisors \(D_1, D_2\) and \(D_3\).
[25]:
using PyPlot
imshow(log10.(abs.(hcat(homSols...))))

[25]:
PyObject <matplotlib.image.AxesImage object at 0x179c7b790>
To illustrate the orthogonal slicing method from Section 4.4 in our paper, we generate some equations of a higher degree on the same toric variety.
[26]:
Pbig = @pm polytope.Polytope(FACETS = hcat(5*a, Fᵀ))
exps = getLatticePoints(Pbig,t)
mons = [prod(t.^exps[ℓ,:]) for ℓ = 1:size(exps,1)]
cmatrix = randn(2,length(mons))
f̂ = [cmatrix[ℓ,:]'*mons for ℓ = 1:2];
We compute a generic start system and track the solution with index “solidx” (chosen randomly) using 10 different random, fixed linear spaces.
[27]:
ĝ,pt = computeStartSystem(f̂)
homStartsols, p, fullHomStartsols = solveGeneric(ĝ, pt, Fᵀ)
γ = randn(ComplexF64)
(k,n) = size(Fᵀ)
@polyvar x[1:4]
f, αf = homogenize(f̂, Fᵀ,x)
g, αg = homogenize(ĝ, Fᵀ,x)
Δt = 0.001 #discretization step for naive homotopy algorithm
numexps = 10
solidx = 37
κᵣ = fill([],numexps)
torsols = fill([],numexps)
linslices = fill([],numexps)
for i = 1:numexps
println("random slice number " *string(i))
zᵣ, κᵣ[i], linslices[i] = trackPathRandom(f, g, x, fullHomStartsols[solidx], n, k, Δt,γ)
if i == 1
semilogy(1:-Δt:Δt,κᵣ[i],color = "#14edc5", label="random slice")
else
semilogy(1:-Δt:Δt,κᵣ[i],color = "#14edc5")
end
println("residual: " * string(get_residual(f,[zᵣ],x)))
torsols[i] = dehomogenizeSolutions([zᵣ], Fᵀ)
println("computed solution: " * string(torsols[i]))
end
legend()
xlabel("τ")
ylabel("κ")
meanκ = 10.0.^(sum([log10.(ℓ) for ℓ ∈ κᵣ])/numexps);
Tracking 100 paths... 100%|█████████████████████████████| Time: 0:00:03
# paths tracked: 100
# non-singular solutions (real): 100 (0)
# singular endpoints (real): 0 (0)
# total solutions (real): 100 (0)
random slice number 1
residual: [6.998545843998793e-16]
computed solution: Any[Complex{Float64}[-0.6117540699006732 + 1.1173490185265766im, 0.10427879561967845 - 0.48453334544885374im]]
random slice number 2
residual: [7.913460749832348e-16]
computed solution: Any[Complex{Float64}[-0.6117540699006729 + 1.1173490185265764im, 0.10427879561967861 - 0.48453334544885385im]]
random slice number 3
residual: [1.272576505061576e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006734 + 1.1173490185265766im, 0.1042787956196785 - 0.4845333454488539im]]
random slice number 4
residual: [6.932547917400776e-16]
computed solution: Any[Complex{Float64}[-0.6117540699006732 + 1.1173490185265766im, 0.10427879561967854 - 0.48453334544885385im]]
random slice number 5
residual: [1.0122187521958066e-17]
computed solution: Any[Complex{Float64}[-0.6117540699006729 + 1.1173490185265766im, 0.10427879561967854 - 0.48453334544885374im]]
random slice number 6
residual: [1.338705067677462e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006732 + 1.1173490185265766im, 0.10427879561967861 - 0.4845333454488538im]]
random slice number 7
residual: [2.518462780271932e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006734 + 1.1173490185265766im, 0.10427879561967832 - 0.4845333454488536im]]
random slice number 8
residual: [2.7865124190615667e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006732 + 1.1173490185265764im, 0.10427879561967845 - 0.4845333454488537im]]
random slice number 9
residual: [2.333360640446741e-15]
computed solution: Any[Complex{Float64}[-0.6117540699006736 + 1.1173490185265766im, 0.1042787956196782 - 0.4845333454488536im]]
random slice number 10

residual: [1.4710097636599483e-16]
computed solution: Any[Complex{Float64}[-0.6117540699006733 + 1.1173490185265764im, 0.10427879561967829 - 0.48453334544885374im]]
We track the same path using orthogonal slicing.
[28]:
NF = smith(Fᵀ);
PP = round.(Int,inv(NF.S));
@time zₒ, κₒ = trackPathOrthogonal(f, g, x, fullHomStartsols[solidx], PP, n, k, Δt, γ);
torsolorth = dehomogenizeSolutions([zₒ], Fᵀ);
println("computed solution: " * string(torsolorth))
5.568242 seconds (54.64 M allocations: 4.059 GiB, 18.42% gc time)
computed solution: Array{Complex{Float64},1}[[-0.6117519112530363 + 1.1173471479487864im, 0.10427739646186379 - 0.4845344267016235im]]
We plot the results.
[29]:
semilogy(1:-Δt:Δt,κₒ,color = "#0f87bf", label = "orthogonal slice")
semilogy(1:-Δt:Δt,meanκ, color = "#ff0000", label = "mean random slice")
xlabel("τ")
ylabel("κ")
#legend(["orthogonal slice", "random slice"])
legend()

[29]:
PyObject <matplotlib.legend.Legend object at 0x179343c40>