Iteraciones en $\mathbb{R}^2$

Dada una aplicación $f:U\subset\mathbb{R}^2\rightarrow\mathbb{R}^2$, suficientemente regular en el dominio $U$, vamos a utilizar la notación $f(x,y)=(u(x,y),v(x,y))$ con $u,v:U\subset\mathbb{R}^2\rightarrow\mathbb{R}$ suficientemente regulares. Denotemos por $Df(x_0,y_0)$ a la matriz jacobiana de $f$ en el punto $(x_0,y_0)$, ésto es $Df(x_0,y_0)=\left(\begin{array}{cc}u_x(x_0,y_0)&u_y(x_0,y_0)\\ v_x(x_0,y_0)&v_y(x_0,y_0)\end{array}\right)$. Supongamos que $Df$ puede invertirse en $U$, entonces el método de Newton se define mediante

$$\left(\begin{array}{cc}x_{k+1}\\ y_{k+1}\end{array}\right)=\... ...(x_k,y_k)\\ v(x_k,y_k)\end{array}\right),\quad (x_0,y_0)\in U,$$ (1)

donde hemos supuesto que los iterados $(x_k,y_k)$ se encuentran en $U$.

Figura 2: Curvas de nivel $u(x,y)=0$, $v(x,y)=0$ (punteada) y su intersección

Daremos varias interpretaciones geométricas de cómo se obtiene $(x_{k+1},y_{k+1})$ a partir de $(x_k,y_k)$. Tomemos para ello la aplicación ejemplo $f(x,y)=(x^3-3xy^2-1,3x^2y-y^3)$(¹), que tiene por ceros las raíces terceras de $1$. Gráficamente, los ceros pueden localizarse mediante la intersección de las curvas de nivel $u(x,y)=x^3-3xy^2-1=0$ y $v(x,y)=3x^2y-y^3=0$. La Figura 2 muestra estos detalles. Los comandos usados para estos dibujos son:

(* Definici'on de las funciones componentes *)
u[x_, y_] := x^3 - 3.*x*y^2 - 1.;
v[x_, y_] := 3.*x^2*y - y^3;

(* Dibujo de la curva de nivel u(x,y)==0 *)
CNu = ImplicitPlot[u[x, y] == 0., {x, -1.2, 1.2}, {y, -1.2, 1.2}]

(* Dibujo de la curva de nivel v(x,y)==0, punteada *)
CNv = ImplicitPlot[v[x, y] == 0., {x, -1.2, 1.2}, {y, -1.2, 1.2}, 
       PlotStyle -> {Dashing[{0.01, 0.01}]}, PlotPoints -> 75]

(* Gr'afico de las tres ra'ices de 1 *)
punts = Show[
    Graphics[Join[{PointSize[0.04]}, 
                  Table[Point[{Cos[k*2*Pi/3], Sin[k*2*Pi/3]}], {k, 0, 2}]
                 ]
            ]]

(* Dibujo uni'on de los tres gr'aficos anteriores *)
Show[CNu, CNv, punts, Axes -> None, Frame -> True]

La siguiente orden nos da el punto $N_f(x_0,y_0)$:

(* Siguiente iterado *)
Nf[u_, v_, x0_, y0_] := Module[{deter, ux, uy, vx, vy, uu, vv},
    uu = Function[{s, t}, u /. {x -> s, y -> t}];
    vv = Function[{s, t}, v /. {x -> s, y -> t}];
    ux = Function[{s, t}, D[u, x] /. {x -> s, y -> t}];
    uy = Function[{s, t}, D[u, y] /. {x -> s, y -> t}];
    vx = Function[{s, t}, D[v, x] /. {x -> s, y -> t}];
    vy = Function[{s, t}, D[v, y] /. {x -> s, y -> t}];
    deter = ux[x0, y0]*vy[x0, y0] - vx[x0, y0]*uy[x0, y0];
    {x0 - (vy[x0, y0]*uu[x0, y0] - uy[x0, y0]*vv[x0, y0])/deter, 
     y0 - (ux[x0, y0]*vv[x0, y0] - vx[x0, y0]*uu[x0, y0])/deter}
  ]

Las expresiones que aparecen en el código se han obtenido a partir de (1) o, equivalentemente, del sistema lineal compatible determinado con dos incógnitas $x_{k+1}$ y $y_{k+1}$:

$$Df(x_k,y_k)\left(\begin{array}{cc}x_{k+1}-x_k\\ y_{k+1}-y_k\... ...ft(\begin{array}{cc}u(x_k,y_k)\\ v(x_k,y_k)\end{array}\right).$$

Si denotamos $\Delta_k=\det Df(x_k,y_k)$ y $p_k=(x_k,y_k)$, entonces las expresiones de $x_{k+1}$ y $y_{k+1}$ vienen dadas por:

$$x_{k+1} = x_k-\frac1{\Delta_k}(v_y(p_k)u(p_k)-u_y(p_k)v(p_k)),$$

$$y_{k+1} = y_k-\frac1{\Delta_k}(u_x(p_k)v(p_k)-v_x(p_k)u(p_k)).$$ (2)

A partir de Nf[ ], que sólo nos da un iterado a partir del anterior, podemos obtener una lista de tantos iterados como queramos:

(* Lista de iterados, utiliza Nf[ ] *)
ListNf[u_, v_, x0_, y0_, N_] := Module[{llista, auxil},
    llista = {{x0, y0}};
    auxil = {x0, y0};
    Do[auxil = Nf[u, v, auxil[[1]], auxil[[2]]]; 
       llista = Join[llista, {auxil}], {N}];
    llista
    ]

Figura 3: Varios iterados de $N_f$ con $f(x,y)=(x^3-3xy^2-1,3x^2y-y^3)$

Figura 4: Iteraciones de $N_f$ para distintos puntos iniciales

Finalmente, podemos dibujar estos puntos mediante la orden

(* Dibuja los iterados, utiliza ListNf[ ] *)
DibListNf[u_, v_, x0_, y0_, N_, opt___] := Module[{dib1, dib2, llista},
    llista = ListNf[u, v, x0, y0, N];
    dib1 = ListPlot[llista, PlotJoined -> True, DisplayFunction -> Identity];
    dib2 = Show[Graphics[Map[Function[t, Disk[t, 0.03]], llista]], 
                DisplayFunction -> Identity];
    Show[dib1, dib2, opt, DisplayFunction -> $DisplayFunction]
    ]

La gráfica izquierda de la Figura 3, se ha obtenido mediante

DibListNf[u[x, y], v[x, y], 0.3, 0.5, 3, AspectRatio -> 1]

La gráfica derecha de la misma figura se ha generado usando la orden Show[ ] para unir el dibujo de DibListNf[ ] con el tercer dibujo de la Figura 2.

La Figura 4 muestra cinco iteraciones de $N_f$ empezando en los puntos $(0.5,0.3)$, $(-0.5,0.3)$ y $(-0.5,-0.3)$.