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!