Consider the combinatorial rational generating function $$\frac{1}{1 - x - y} = \sum_{j, k \geq 0} \binom{j+k}{j} x^j y^k.$$ The critical point equations are $$ \begin{split} 1 = x + y \\ x = y \end{split} $$ and there is only one solution, $(1/2, 1/2)$ which is also minimal.
using HomotopyContinuation
@polyvar x y z
H = 1 - x - y
r, s = 1, 1
sys = System(
[H, s*x*differentiate(H, x) - r*y*differentiate(H, y)],
variables=[x, y]
)
solve_result = solve(sys)
Result with 1 solution ====================== • 1 path tracked • 1 non-singular solution (1 real) • random_seed: 0xb43ea4d8 • start_system: :polyhedral
Homotopy continuation is a numerical method, but there are methods for ensuring that a solution exists using interval arithmetic. Once HC.jl solves the system, we may ask it to return intervals $(a, b)$ such that there is exactly one solution inside $(a, b)$. This is done by a generalization of Newton's method to intervals arithmetic.
sols = solutions(solve_result) # get the list of solutions
certify_result = certify(sys, sols) # certify the solutions
println(certify_result)
certs = certificates(certify_result) # datastructure containing certification result
intervals = certified_solution_interval.(certs)
CertificationResult =================== • 1 solution candidates given • 1 certified solution intervals (1 real, 0 complex) • 1 distinct certified solution intervals (1 real, 0 complex)
1-element Array{Arblib.AcbMatrix,1}:
[[0.500000000000 +/- 1.14e-13] + [+/- 1.14e-13]im; [0.500000000000 +/- 1.14e-13] + [+/- 1.14e-13]im]
To verify that the point $(1/2, 1/2)$ is minimal, we solve the extended critical point system with the minimality equation $H(tx, ty) = 0$. However, notice that $(1/2, 1/2, 1)$ is a solution. Since HC.jl returns intervals, it is impossible to verify that the $t$ coordinate is exactly 1. Therefore, we add an extra equation $u(1-t) - 1$ to ensure no solution with $t = 1$ exists.
@polyvar t u
sys = System([ # full system
H,
s*x*differentiate(H, x) - r*y*differentiate(H, y),
H([x, y] => t*[x, y]),
u*(1-t) - 1
], variables=[x, y, t, u])
# get the certified solution intervals
solve(sys) # we get 0 solutions
Result with 0 solutions ======================= • 1 path tracked • 0 non-singular solutions (0 real) • random_seed: 0x9fce91e7 • start_system: :polyhedral
@polyvar λ
H = (1-x-y)*(20-x-40*y)-1
sys = System(# get the critical points
[H; [x, y] .* differentiate(H, [x, y]) .- λ],
variables=[x, y, λ]
)
certs = certificates(certify(sys, solutions(solve(sys))))
crits = solution_approximation.(certs)
Tracking 4 paths... 100%|███████████████████████████████| Time: 0:00:01 # paths tracked: 4 # non-singular solutions (real): 4 (2) # singular endpoints (real): 0 (0) # total solutions (real): 4 (2)
4-element Array{Array{Complex{Float64},1},1}:
[9.997110519821023 + 0.0im, 0.2527749731647174 + 0.0im, 93.55290965306227 + 0.0im]
[0.490149016126495 - 0.28079160650357426im, 0.5808033325272964 + 0.16235478163384356im, 3.5706646451470894 + 1.9223315807277772im]
[0.490149016126495 + 0.28079160650357426im, 0.5808033325272964 - 0.16235478163384356im, 3.5706646451470894 - 1.9223315807277772im]
[0.548232473567013 + 0.0im, 0.30997733613966416 + 0.0im, -3.9442389433564395 + 0.0im]
sys = System([ # full system
H,
s*x*differentiate(H, x) - r*y*differentiate(H, y),
H([x, y] => t*[x, y]),
u*(1-t) - 1
], variables=[x, y, t, u])
# get the certified solution intervals
real_solutions(solve(sys))
Tracking 8 paths... 100%|███████████████████████████████| Time: 0:00:01 # paths tracked: 8 # non-singular solutions (real): 4 (2) # singular endpoints (real): 0 (0) # total solutions (real): 4 (2)
2-element Array{Array{Float64,1},1}:
[0.5482324735670131, 0.3099773361396641, 1.7099367491047748, -1.4085761883167662]
[9.997110519821023, 0.2527749731647175, 0.09218565523266332, 1.1015468149011123]
# This code can be entirely reused for the next examples
function check_minimality(H)
full = System([ # full system
H;
[x, y, z] .* differentiate(H, [x, y, z]) .- r*λ;
H([x, y, z] => t*[x, y, z]);
u*(1-t) - 1
], variables=[x, y, z, λ, t, u])
real_sols = real_solutions(solve(full))
if length(real_sols) > 1
certs = certificates(certify(full, real_sols))
sols = certified_solution_interval.(certs)
println(length(sols)) # this system has potentially many solutions
# here check the t coordinate
for sol in filter(sol -> 0 < abs(sol[5]) < 1, sols)
display(sol) # this mean the point is NOT minimal
end
end
end
check_minimality (generic function with 1 method)
The number of walks of length $n$ in $\mathbb{Z}^2$ which begin at the origin and take steps in the set $$S = \left\{ (1, 0), (-1, 0), (0, 1), (0, -1) \right\}$$ is the coefficient on the main diagonal of a rational generating function with denominator $$H = 1 - (zx^2y + zxy^2 + zx + zy).$$
H = 1 - (z*x^2*y + z*x*y^2 + z*x + z*y)
# get the critical points
sys = System(
[H; [x, y, z] .* differentiate(H, [x, y, z]) .- λ],
variables=[x, y, z, λ]
)
certs = certificates(certify(sys, solutions(solve(sys))))
crits = solution_approximation.(certs)
Tracking 4 paths... 100%|███████████████████████████████| Time: 0:00:01 # paths tracked: 4 # non-singular solutions (real): 2 (2) # singular endpoints (real): 0 (0) # total solutions (real): 2 (2)
2-element Array{Array{Complex{Float64},1},1}:
[-0.9999999999999998 + 0.0im, -0.9999999999999998 + 0.0im, -0.2500000000000001 + 0.0im, -1.0 + 0.0im]
[0.9999999999999998 + 0.0im, 0.9999999999999998 + 0.0im, 0.2500000000000001 + 0.0im, -1.0 + 0.0im]
check_minimality(H)
Tracking 16 paths... 100%|██████████████████████████████| Time: 0:00:01 # paths tracked: 16 # non-singular solutions (real): 6 (2) # singular endpoints (real): 0 (0) # total solutions (real): 6 (2) 2
Both points are minimal and the asymptotic formula is $$O(C + 1/n) \cdot \frac{4^n}{n}.$$
For the step set $S = \{ (1, 1), (1, -1), (-1, 1), (-1, -1) \}$ the denominator will be $$H = 1 - (zx^2y^2 + zx^2 + zy^2 + z)$$
H = 1 - (z*x^2*y^2 + z*x^2 + z*y^2 + z)
sys = System(
[H; [x, y, z] .* differentiate(H, [x, y, z]) .- λ],
variables=[x, y, z, λ]
)
certs = certificates(certify(sys, solutions(solve(sys))))
crits = solution_approximation.(certs)
Tracking 8 paths... 100%|███████████████████████████████| Time: 0:00:01 # paths tracked: 8 # non-singular solutions (real): 4 (4) # singular endpoints (real): 0 (0) # total solutions (real): 4 (4)
4-element Array{Array{Complex{Float64},1},1}:
[-1.0 + 0.0im, 1.0 + 0.0im, 0.25 + 0.0im, -1.0000000000000002 + 0.0im]
[-1.0 + 0.0im, -1.0 + 0.0im, 0.25 + 0.0im, -1.0000000000000002 + 0.0im]
[1.0 + 0.0im, 1.0 + 0.0im, 0.25 + 0.0im, -1.0000000000000002 + 0.0im]
[1.0 + 0.0im, -1.0 + 0.0im, 0.25 + 0.0im, -1.0000000000000002 + 0.0im]
check_minimality(H)
Tracking 40 paths... 100%|██████████████████████████████| Time: 0:00:01 # paths tracked: 40 # non-singular solutions (real): 16 (0) # singular endpoints (real): 0 (0) # total solutions (real): 16 (0)
Again, all critical points are minimal. The asymptotic formula is the same $$O(C + 1/n) \cdot \frac{4^n}{n}.$$
The denominator of the generating for bi-colored supertrees is $$H = x^5y^2 + 2x^2y - 2x^3y + 4y + x - 2.$$
H = x^5*y^2 + 2*x^2*y - 2*x^3*y + 4*y + x - 2
sys = System(
[H; x*differentiate(H, x) - y*differentiate(H, y)],
variables=[x, y]
)
certs = certificates(certify(sys, solutions(solve(sys))))
crits = certified_solution_interval.(certs)
Tracking 7 paths... 100%|███████████████████████████████| Time: 0:00:01 # paths tracked: 7 # non-singular solutions (real): 2 (2) # singular endpoints (real): 3 (3) # total solutions (real): 5 (5)
2-element Array{Arblib.AcbMatrix,1}:
[[3.2360679775 +/- 2.05e-11] + [+/- 2.03e-11]im; [0.047745751406 +/- 8.16e-13] + [+/- 5.53e-13]im]
[[-1.23606797750 +/- 1.22e-12] + [+/- 1.01e-12]im; [0.327254248594 +/- 6.70e-13] + [+/- 4.07e-13]im]
println(H(2, 1/8)) # this point is missed.
CE = x*differentiate(H, x) - y*differentiate(H, y)
println(CE(2, 1/8))
0.0 0.0
solve(sys)
Result with 4 solutions ======================= • 7 paths tracked • 2 non-singular solutions (2 real) • 2 singular solutions (2 real) • random_seed: 0x6182306a • start_system: :polyhedral • multiplicity table of singular solutions: ╭───────┬───────┬────────┬────────────╮ │ mult. │ total │ # real │ # non-real │ ├───────┼───────┼────────┼────────────┤ │ 1 │ 1 │ 1 │ 0 │ │ 2 │ 1 │ 1 │ 0 │ ╰───────┴───────┴────────┴────────────╯