Can you please help me whether $N-R$ method or other methods converges?

Question

Consider the system of equations \begin{align} &f(x,y)=x+\frac{y^4}{2}+\frac{x^{32}}{4}+\frac{y^{128}}{8}=0 \\ &g(x,y)=y+\frac{x^8}{2}+\frac{y^{32}}{4}+\frac{x^{256}}{8}=0. \end{align}

I want to solve it using Newton-Raphson process or any other methods

Consider the initial guess $(x_0,y_0)=((-2)^{5/31},(-2)^{9/31})$, which is either $( 1.118,-1.223 )$ or $(x_0,y_0)=(-1.118,1.223)$.

I want to see whether the Newton-Raphson method converges with this initial guess.

Note that the initial guess $(x_0,y_0)$ is a simultaneous zero of the following truncated system obtained from the original system \begin{align} &x+\frac{y^4}{2}=0\\ &y+\frac{x^8}{2}=0. \end{align}

By hand calculation, it is laborious. For the jacobian is $$J =\begin{pmatrix} 1+8x^{31} & 2y^3+16y^{127} \\ 2x^3+32x^{255} & 1+8y^{31} \end{pmatrix} \Rightarrow J((x_0,y_0)) \approx \begin{pmatrix} 255 & -2.03 \times 10^{12} \\ 7.21 \times 10^{13} & -4104\end{pmatrix}$$ So the 2nd iteration is given by \begin{align} \begin{pmatrix} x_1 \\ y_1 \end{pmatrix}&=\begin{pmatrix}1.118 \\ 1.223 \end{pmatrix} -J((x_0,y_0))^{-1} \begin{pmatrix}f((x_0,y_0)) \\ g((x_0,y_0)) \end{pmatrix} \end{align} which seems difficult to calculate because I can not invert the huge matrix.

The zeros of the truncated system should converge to the solutions of the original system.

Can you please help me whether N-R method or any other numerical methods converges ?

More, specifically, how to show the simultaneous zeroes of the truncated system $ x+\frac{y^4}{2}=0=y+\frac{x^8}{2}$ converges to the simultaneous zeroes of the original system $f(x,y)=0=g(x,y)$ ?

Thanks

Answer 1

Here is some sage code starting the iteration in $(-1, -1)$.

IR = RealField(150)
var('x,y');

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])

---

The printed results are following:

x0 = -1.0000000000000000000000000000000000000000000
y0 = -1.0000000000000000000000000000000000000000000
f(x0, y0) = -0.12500000000000000000000000000000000000000000
g(x0, y0) = -0.12500000000000000000000000000000000000000000

x1 = -1.0024422735346358792184724689165186500888099
y1 = -1.0059946714031971580817051509769094138543517
f(x1, y1) = 0.048583106252741437695070286532161290853457141
g(x1, y1) = 0.035005003661522385716279832064793355617823982

x2 = -1.0020504863163137765771033606257451944731308
y2 = -1.0047357187454832971083384658776119569301797
f(x2, y2) = 0.0032718331735723517740747794193782531073714280
g(x2, y2) = 0.0013588050706032233978301227450805778637631180

x3 = -1.0020413711644051490584615403424219543507335
y3 = -1.0046329944054482289164004688385874135513553
f(x3, y3) = 0.000019403097353192714808877039288506402587269048
g(x3, y3) = 2.0763587244920856969768249854671070482986312e-6

x4 = -1.0020414289196482740462473476044443285649078
y4 = -1.0046323504589185162361046727090015558601915
f(x4, y4) = 7.5595118968724701695743900080517263081615759e-10
g(x4, y4) = 8.1922112283336376999011329830106254188700964e-11

x5 = -1.0020414289218795422093499868226092493754357
y5 = -1.0046323504338328451030885177269370761065741
f(x5, y5) = 1.1471443918740692095841154251012109118025975e-18
g(x5, y5) = 1.2374987496689888382750646154304742963132461e-19

x6 = -1.0020414289218795422127464111069351864618923
y6 = -1.0046323504338328450650188283625352002492464
f(x6, y6) = 2.6419634883773767349590310873145265485723181e-36
g(x6, y6) = 2.8548246676615458259051746549294473950928253e-37

x7 = -1.0020414289218795422127464111069351864697057
y7 = -1.0046323504338328450650188283625352001615711
f(x7, y7) = -5.6051938572992682836949183331596645251210478e-45
g(x7, y7) = 1.4538471567369977110833694426632879862032718e-44

x8 = -1.0020414289218795422127464111069351864697057
y8 = -1.0046323504338328450650188283625352001615711
f(x8, y8) = -5.6051938572992682836949183331596645251210478e-45
g(x8, y8) = 1.4538471567369977110833694426632879862032718e-44

The value $(x_7,y_7)$ is kept at the next iteration step.

Answer 2

By performing the contour plot for some truncated approximations, we can infer for which the NR algorithm will work.

$$ \cases{ x + \frac{y^9}{3}+\frac{x^{243}}{9}=0 \\ y + \frac{x^{27}}{3}+\frac{y^{243}}{9}=0 \\ } $$

Here NR works with one solution.

$$ \cases{ x + \frac{y^9}{3}+\frac{x^{243}}{9} + \frac{y^{2187}}{27}=0 \\ y + \frac{x^{27}}{3}+\frac{y^{243}}{9} + \frac{x^{6561}}{27}=0 \\ } $$

Here NR succeeds with three solutions

$$ \cases{ x + \frac{y^9}{3}+\frac{x^{243}}{9} + \frac{y^{2187}}{27}+\frac{x^{59049}}{81}=0 \\ y + \frac{x^{27}}{3}+\frac{y^{243}}{9} + \frac{x^{6561}}{27}+\frac{y^{59049}}{81}=0 \\ } $$

Here NR only gives one solution. The two contours only cross one time.