Numerical Nonlinear Algebra

Numerical methods can also be used for computing the dimension of the vector spaces we present in the introduction by using Theorem 5.1. The following code in julia provides an example for this by using the package HomotopyContinuation.jl. We compute the Euler characteristic of the very affine surface \(X\) from Example 2.7, with \(f\) as in Equation (2.8) as the number of solution to a system of rational functions in the variables \(x,y\). The result in this case is \(6\).

using HomotopyContinuation
@var x y s ν[1:2]
f = -x*y^2 + 2*x*y^3 + 3*x^2*y - x^2*y^3 - 2*x^3*y + 3*x^3*y^2
L = s*log(f) + ν[1]*log(x) + ν[2]*log(y)
F = System(differentiate(L,[x;y]), parameters = [s;ν])
monodromy_solve(F)

<span class=”ansi-green-fg”>Solutions found: 6 Time: 0:00:02</span>

<span class=”ansi-blue-fg”> tracked loops (queued): 36 (0)</span> <span class=”ansi-blue-fg”> solutions in current (last) loop: 0 (0)</span> <span class=”ansi-blue-fg”> generated loops (no change): 6 (5)</span> </pre>

textcolor{ansi-green}{Solutions found: 6 Time: 0:00:02}

textcolor{ansi-blue}{ tracked loops (queued): 36 (0)} textcolor{ansi-blue}{ solutions in current (last) loop: 0 (0)} textcolor{ansi-blue}{ generated loops (no change): 6 (5)} end{sphinxVerbatim}

Solutions found: 6 Time: 0:00:02

 tracked loops (queued): 36 (0)  solutions in current (last) loop: 0 (0)  generated loops (no change): 6 (5)

MonodromyResult
===============
• return_code → :heuristic_stop
• 6 solutions
• 36 tracked loops
• random_seed → 0x287750da

The following \(\texttt{julia}\) code is the one used for the experiments in Section 5 of the article. It consists of six functions, copied below with a few lines of documentation. One can also find it in Appendix B in the paper.

# One step of Euler’s method on dy/dx = ω(x)y, with stepsize Δx.
function euler_step(x,y,Δx,ω)
newy = (1+ω(x)*Δx)*y
return (x+Δx,newy)
end

# One step of Newton iteration on y^k - prod(f(x).^(k*s))*x^(k*ν).
function newton_step(y, x, f, s, ν, k)
Fy = y^k - prod(f(x).^(k*s))*x^(k*ν)
Δy = -Fy/(k*y^(k-1))
return y + Δy
end

# Compute function values at N equidistant nodes on the line segment SxTx.
# The initial condition is given by Sy.

function track_line_segment(Sx,Sy,Tx,N,f,ω,s,ν,k)
xx = zeros(ComplexF64,N)
yy = zeros(ComplexF64,N)
xx[1] = Sx
yy[1] = Sy
Δx = (Tx - Sx)/(N-1)
for i = 2:N
(xx[i], yi_tilde) = euler_step(xx[i-1],yy[i-1],Δx,ω)
for j = 1:4
yi_tilde = newton_step(yi_tilde, xx[i], f, s, ν, k)
end
yy[i] = copy(yi_tilde)
end
return xx, yy
end


# Numerical integration based on function values yy and interval length h.
function integrate_trapezoidal(yy,h)
return (yy[1]/2 + sum(yy[2:end-1]) + yy[end]/2)*h
end
# Compute the integral over the line segment AB with N discretization nodes.
# The initial condition is given by phiAB_at_A.
# The integral is computed for the cocycles in aabb.
function integrate_line_segment(A,phiAB_at_A,B,N,f,ω,s,ν,k,aabb)
IAB = []
xxAB, yyAB = track_line_segment(A,phiAB_at_A,B,N,f,ω,s,ν,k)
for ab in aabb
a = ab[1]; b = ab[2];
single_valued = [x^(b-1)*prod(f(x).^a) for x in xxAB]
IAB = push!(IAB,integrate_trapezoidal(yyAB.*single_valued,(B-A)/(N-1)))
end
return IAB,yyAB
end
# Compute the integral over the twisted cycle defined by A, B, C,
# and the initial condition phiAB_at_A.
# The integral is computed for the cocycles in aabb.
function integrate_loop(A,B,C,phiAB_at_A,N,f,ω,s,ν,k,aabb)
IAB, yyAB = integrate_line_segment(A,phiAB_at_A,B,N,f,ω,s,ν,k,aabb)
phiBC_at_B = yyAB[end]
IBC, yyBC = integrate_line_segment(B,phiBC_at_B,C,N,f,ω,s,ν,k,aabb)
phiCA_at_C = yyBC[end]
ICA, yyCA = integrate_line_segment(C,phiCA_at_C,A,N,f,ω,s,ν,k,aabb)
return IAB+IBC+ICA
end

integrate_loop (generic function with 1 method)

Example 5.4

The following lines illustrate how to use the code presented above in the case \(f = (x-1, x-2),\) \(s = (\frac{1}{2}, \frac{1}{2})\), and \(\nu = \frac{1}{2}\)

f = x -> [x-1; x-2]
s = [1/2;1/2]; ν = 1/2; N = 1000; k = 2;
ω = x -> s[1]/(x-1) + s[2]/(x-2) + ν/x
cocycles = [[[-1;0],1], [[0,-1],1], [[0,0],0]]
A1 = 1/2+im; B1 = 1/2-im; C1 = 3+0im
phiA1B1_at_A1 = A1^ν*prod(f(A1).^s)
I1 = integrate_loop(A1,B1,C1,phiA1B1_at_A1,N,f,ω,s,ν,k,cocycles)
A2 = -1+0im; B2 = 3/2+im; C2 = 3/2-im
phiA2B2_at_A2 = A2^ν*prod(f(A2).^s)
I2 = integrate_loop(A2,B2,C2,phiA2B2_at_A2,N,f,ω,s,ν,k,cocycles)

3-element Vector{ComplexF64}:
   3.496076532931323 - 3.2862601528904634e-14im
  0.6482432272682694 - 4.529709940470639e-14im
 -4.1443176753269615 - 7.549516567451064e-15im