Can you help with a powerful Sage code to solve the system in N-R method or other numerical methods ?

asked 2022-04-30 08:02:44 +0200

anonymous user


updated 2022-05-03 02:41:27 +0200

dan_fulea gravatar image

I want to approximate the zeros of the following nonlinear system using Newton-Raphson process/ any other method. $$F(x,y)=x+\frac{y^{3^2}}{3}+\frac{x^{3^5}}{3^2}+\frac{y^{3^7}}{3^3}+\frac{x^{3^{10}}}{3^4}=0,$$ $$G(x,y)=y+\frac{x^{3^3}}{3}+\frac{y^{3^5}}{3^2}+\frac{x^{3^{8}}}{3^3}+\frac{y^{3^{10}}}{3^4}=0.$$

For the following easier system $$f(x,y)=x+\frac{y^{2^2}}{2}+\frac{x^{2^5}}{2^2}+\frac{y^{2^7}}{2^3}=0 , $$ $$g(x,y)=y+\frac{x^{2^3}}{2}+\frac{y^{2^5}}{2^2}+\frac{x^{2^8}}{2^3}=0,$$ we can use the initial guess $(-1,-1)$, and checked that the N-R process converges at $8$th iteration.

IR = RealField(150)


f = x + y^4/2 + x^32/4 + y^128/8
g = y + x^4/2 + y^32/4 + x^256/8

a, b = IR(-1), IR(-1)

J = matrix(2, 2, [diff(f, x), diff(f, y), diff(g, x), diff(g, y)])
v = vector(IR, 2, [IR(a), IR(b)])

for k in range(9):

    a, b = v[0], v[1]
    fv, gv = f.subs({x : a, y : b}), g.subs({x : a, y : b})
    print(f'x{k} = {a}\ny{k} = {b}')
    print(f'f(x{k}, y{k}) = {fv}\ng(x{k}, y{k}) = {gv}\n')
    v = v - J.subs({x : a, y : b}).inverse() * vector(IR, 2, [fv, gv])

That is all about the code.

Unfortunately the same code seems to be inefficient to solve the above system

$$F(x,y)=0,$$ $$G(x,y)=0.$$

It seems this system is associated with higher degrees.

I think this code needs to be modified.

Can you help with a powerful Sage code to solve the system in N-R method or other numerical methods ?

edit retag flag offensive close merge delete


One can use Sage's infrastructure to Singular :

R1.<u,v> = QQ[]
Sys2 = [u + v^2^2/2 + u^2^5/2^2 + v^2^7/2^3,
        v + u^2^3/2 + v^2^5/2^2 + u^2^8/2^3]
Sys3 = [u + v^3^2/3 + u^3^5/3^3 + v^3^7/3^3 + u^3^10/3^4,
        v + u^3^3/3 + v^3^5/3^2 + u^3^8/3^3 + v^3^10/3^4]
J2 = R1.ideal(*Sys2)
J3 = R1.ideal(*Sys3)

but, if this allows to show that these systems have both at most a finite number of roots :

sage: J2.dimension()
sage: J3.dimension()

even Sys2seems too hard to compute tjhem : %time Sol2 = J2.variety(ring=RR) ran for more than 6 hours without result and used more than 13 GB of RAM and a ...(more)

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2022-05-01 23:44:49 +0200 )edit

Please always mention the source, if any, and the reason for the post. In this case, the posted question is some amalgamated result of: stackexchange question 4438130 and some other post on stackexchange that i lost. Your problem is not the code, but the mathematical part, the point to get started from. As in the case of 4428139, your point to start the iteration was leading to an explosion of the iterates, taking $(-1,-1)$ as starting point the Newton iterations converge. Now you have a solution, that was not so easily obtained. Immediately we have one more question. Please tell us why do you need to solve the above, if the motivation is there, i will try...

dan_fulea gravatar imagedan_fulea ( 2022-05-03 02:50:29 +0200 )edit