Help Debugging my first Julia Conversion of Maple Code

I’m trying to replicate the behavior of a Maple code snippet (from [1] on Page 15) in Julia, as I’m new to Julia I am not sure if this is being implemented correctly and a particular line of my code is running super slowly. Thank you for checking my post in advance!

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 in Julia


using SymPy

# Define the symbolic variables
@syms t x y

# Define the function h
h = t - x^5 + x^7 - x^11 * t

# Substitute x = y * t^(1/5) into h
hy = subs(h, x => y * t^(1/5))

# Simplify hy / t
hyt = simplify(hy / t)

# Define the Newton iteration function
function newton(x_val)
    diff_hyt = diff(hyt, y)
    return x_val - subs(hyt, y => x_val) / subs(diff_hyt, y => x_val)
end

# Initialize x_vals as symbolic list with the first element as a symbolic 1
x_vals = [Sym(1)]  # Start with symbolic 1

# Array to store series expansions
s_vals = []

# Perform 6 iterations of Newton's method
for k in 1:6
    # Perform Newton's method step
    x_next = newton(x_vals[end])
   push!(x_vals, simplify(x_next))  # Append the simplified result
    
    # Compute series expansion up to 15 terms
    s_current = series(x_vals[end], t, 0, 15)
    push!(s_vals, s_current)
    
    # Print the difference between successive series expansions
    if k > 1
        diff_series = s_vals[k] - s_vals[k-1]
        println("Difference in series at iteration $k: ", numerator(diff_series))
    end
end

Issue

The code takes quite long…

push!(x_vals, simplify(x_next))

References

[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.


Any insights or suggestions on how to resolve this issue would be greatly appreciated. Thanks!

1 Like

Isn’t SymPy.jl the Julia bindings for the Python package?

You may have better luck (and performance) using a native Julia library like Symbolics.jl.

3 Likes

It’s safe to assume that push! won’t take too long, so the question is why simplify takes so long. To diagnose the problem, you could print out x_next to see what’s being simplified.

Sorry for my late response. I should’ve mentioned but, yeah, I’d looked into it and had found that the expression was rather longer… But still compared to SageMath script that I’d written it was too long to iterate through so I aksed the question. Anyway, I’d try to rewrite this using Symbolics.jl!

Got it! I will try this!

1 Like