Converting Homotopy continuation in Maple into SageMath

asked 2024-09-06 09:56:47 +0200

Rowing0914 gravatar image

updated 2024-09-23 10:30:49 +0200

FrédéricC gravatar image

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.

edit retag flag offensive close merge delete

Comments

maybe this could help

sage: A=PuiseuxSeriesRing(QQ,'t')
sage: t=A.gen()
sage: f=t**(-1/5)+4+t
sage: f**3
t^(-3/5) + 12*t^(-2/5) + 48*t^(-1/5) + 64 + 3*t^(3/5) + 24*t^(4/5) + 48*t + 3*t^(9/5) + 12*t^2 + t^3
FrédéricC gravatar imageFrédéricC ( 2024-09-06 15:52:09 +0200 )edit

@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)?

Rowing0914 gravatar imageRowing0914 ( 2024-09-07 14:45:25 +0200 )edit

The Puiseux ring handles the precision automatically :

sage: A=PuiseuxSeriesRing(QQ,'t', default_prec=15)
sage: t=A.gen()
sage: f=(t**(-1/3)+4+t**6).add_bigoh(15)
sage: f**2
t^(-2/3) + 8*t^(-1/3) + 16 + 2*t^(17/3) + 8*t^6 + t^12 + O(t^(44/3))
FrédéricC gravatar imageFrédéricC ( 2024-09-07 19:59:27 +0200 )edit