Revision history [back]

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.