Investigating branches of the square-root function.
Requirements: ore_algebra package
# Define the rings of differential and shift operators
from ore_algebra import *
Ind.<n> = PolynomialRing(QQ); Shift.<Sn> = OreAlgebra(Ind)
Pols.<z> = PolynomialRing(QQ); Diff.<Dz> = OreAlgebra(Pols)
# Any branch of f(z) = sqrt(z) satisfies f'(z) = 1/2/sqrt(z) = f/2/z, so any branch of the square-root satisfies the D-finite equation
diff = 2*z*Dz - 1
# Starting with a fixed value of the square-root at a (non-zero) starting point we can locally compute branches of the square-root using numeric analytic continuation
# For instance, starting with f(1) = 1 we can analytically continue f in a line from 1 to 2 to get the value of sqrt(2)
# (The two branches of sqrt(z) must have f(1)=1 or f(1)=-1 -- other values of f(1) gives solutions of the D-finite equation which don't represent a square-root)
diff.numerical_solution(ini=[1], path=[1, 2], eps=1e-100)
[1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273 +/- 5.11e-102]
# This is the "usual" (principal branch) value of sqrt(2)
sqrt(2).n(333)
1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157
# We can continue this branch in a polygonal path 1 -> I -> -1 to get a branch of sqrt(z) near -1, which has sqrt(-1) = I
diff.numerical_solution(ini=[1], path=[1, I, -1], eps=1e-100)
[+/- 7.87e-103] + [1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 9.95e-103]*I
# If continue our starting branch in a polygonal path 1 -> -I -> -1 we get a different branch of sqrt(z) near -1, which has sqrt(-1) = -I
diff.numerical_solution(ini=[1], path=[1, -I, -1], eps=1e-100)
[+/- 7.87e-103] + [-1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 9.95e-103]*I
# This difference occurs because we must pass a branch cut going in a circle around the origin
# Going in a circle around the origin once gives a branch where sqrt(1)=-1
diff.numerical_solution(ini=[1], path=[1, I, -1, -I, 1], eps=1e-100)
[-1.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 1.12e-103] + [+/- 3.30e-112]*I
# Going in a circle around the origin twice returns us to the original branch where sqrt(1)=1
diff.numerical_solution(ini=[1], path=[1, I, -1, -I, 1, I, -1, -I, 1], eps=1e-100)
[1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 5.59e-104] + [+/- 4.01e-111]*I
# Similar things happen with other algebraic equations. For instance, branches of z^(1/3) satisfy the D-finite equation
diff2 = 3*z*Dz - 1
# Now starting with the branch 1^(1/3) = 1 and going around the origin in a circle multiplies by the root of unity e^(2*pi*I/3)
print(diff2.numerical_solution(ini=[1], path=[1, I, -1, -I, 1], eps=1e-50))
print(exp(2*pi*I/3).n(150))
[-0.50000000000000000000000000000000000000000000000000000 +/- 5.23e-54] + [0.86602540378443864676372317075293618347140262690519031 +/- 7.79e-54]*I -0.50000000000000000000000000000000000000000000 + 0.86602540378443864676372317075293618347140263*I
# Going around the origin in a circle twice multiplies by e^(4*pi*I/3)
print(diff2.numerical_solution(ini=[1], path=[1, I, -1, -I, 1, I, -1, -I, 1], eps=1e-50))
print(exp(4*pi*I/3).n(150))
[-0.50000000000000000000000000000000000000000000000000000 +/- 2.62e-54] + [-0.86602540378443864676372317075293618347140262690519031 +/- 5.18e-54]*I -0.50000000000000000000000000000000000000000000 - 0.86602540378443864676372317075293618347140263*I
# Going around the origin in a circle three times returns to the original branch
print(diff2.numerical_solution(ini=[1], path=[1, I, -1, -I, 1, I, -1, -I, 1, I, -1, -I, 1], eps=1e-50))
[1.0000000000000000000000000000000000000000000000000000 +/- 5.23e-54] + [+/- 2.41e-60]*I
# The function sqrt(z*(2-z)) satisfies the D-finite equation defined by
diff3 = z*(z-2) * Dz + (1-z)
# Analytically continue the branch of sqrt(z*(2-z)) defined by sqrt(1) = 1
# What happens if you circle exactly one of z=0 or z=2? What happens if you circle *both*?