Converting Homotopy continuation in Maple into SageMath
I'm trying to replicate the behaviour of Maple code in [1] on Page 15. But I'm stuck at a bug described below. Any help?
Reference code in Maple
> h := t-x^5 + x^7 - x^(11)*t:
> hy := subs(x = y*t^(1/5),h):
> hyt := simplify(hy/t):
> newton := x -> x - subs(y=x,hyt/diff(hyt,y)):
> x[0] := 1:
> for k from 1 to 6 do
> x[k] := newton(x[k-1]):
> s[k] := series(x[k],t=0,15):
> lprint(op(1,s[k]-s[k-1]));
> end do:
My attempt
# Import necessary modules
from sage.all import *
# Define the polynomial ring
var("x y t")
h = t - x^5 + x^7 - x^11 * t
# Substitute x = y * t^(1/5) into h
hy = h.substitute(x = y * t^(1/5))
# Simplify hy/t
hyt = simplify(hy/t)
# Define Newton's method iteration function
def newton(x_val):
return x_val - (hyt.substitute(y = x_val) / diff(hyt, y).substitute(y = x_val))
# Initial guess for x: Start with a symbolic expression in terms of t
x0 = 1.0
# Perform Newton's method iterations
x_vals = [x0]
for k in range(1, 7):
x_next = newton(x_vals[k-1])
x_vals.append(x_next)
# Compute series expansion around t = 0, up to 15 terms
if k > 1:
s_prev = x_vals[k-1].series(t, 15)
s_curr = x_vals[k].series(t, 15)
# Print the difference between successive series approximations
print(s_curr - s_prev)
this code stucks with the division by zero error as follows;
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In [1], line 30
28 # Compute series expansion around t = 0, up to 15 terms
29 if k > Integer(1):
---> 30 s_prev = x_vals[k-Integer(1)].series(t, Integer(15))
31 s_curr = x_vals[k].series(t, Integer(15))
33 # Print the difference between successive series approximations
File /home/sc_serv/sage/src/sage/symbolic/expression.pyx:4923, in sage.symbolic.expression.Expression.series()
4921 sig_on()
4922 try:
-> 4923 x = self._gobj.expand(0).series(symbol0._gobj, cprec, options)
4924 nex = SymbolicSeries.__new__(SymbolicSeries)
4925 nex._parent = self._parent
ValueError: power::eval(): division by zero
Ref
[1] Sommese, Andrew J., Jan Verschelde, and Charles W. Wampler. "Introduction to numerical algebraic geometry." Graduate School on Systems of Polynomial Equations: From Algebraic Geometry to Industrial Applications (2003): 14-25.
maybe this could help
@FredericC Thank you for your comment. I checked the API of PuiseuxSeriesRing but couldn't find seeming ly relevant method..? How can I use PuiseuxSeriesRing to generate the series upto the 15th term done by
.series(t, 15)
?The Puiseux ring handles the precision automatically :