En este capítulo generalizaremos el problema de hallar un cero de una función real de variable real desarrollado en el Capítulo de ceros del curso Métodos numéricos con Python. Cálculo Numérico.

Concretamente, en lugar de considerar una función \(f:D\subseteq \mathbb{R}\longrightarrow \mathbb{R}\) y usar técnicas para hallar soluciones aproximadas \(\hat{x}\) de la ecuación \(f(x)=0\), consideraremos funciones de \(n\) variables, con \(n>1\), \(\mathbf{F}:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}^n\) y estudiaremos técnicas para hallar soluciones aproximadas \(\hat{\mathbf{x}}\) de la ecuación \(\mathbf{F}(\mathbf{x})=\mathbf{0}\).

Si escribimos \[ \mathbf{F}(\mathbf{x})=\begin{bmatrix}f_1(x_1,\ldots,x_n)\\ f_2(x_1,\ldots,x_n)\\\vdots \\ f_n(x_1,\ldots,x_n)\end{bmatrix}, \] hallar un cero de la ecuación \(\mathbf{F}(\mathbf{x})=\mathbf{0}\) equivale a resolver el siguiente sistema de ecuaciones no lineal: \[ \left. \begin{align*} f_1(x_1,\ldots,x_n) = & 0,\\ f_2(x_1,\ldots,x_n) = & 0,\\ &\vdots\\ f_n(x_1,\ldots,x_n) = & 0. \end{align*} \right\} \]

El sistema anterior muestra que resolver la ecuación \(\mathbf{F}(\mathbf{x})=\mathbf{0}\) puede interpretarse como una generalización de los capítulos del curso Métodos numéricos con Python. Álgebra Lineal Numérica donde se desarrollaban métodos para resolver sistemas lineales \(\mathbf{A}\mathbf{x}=\mathbf{b}\), donde \(\mathbf{A}\) es una matriz \(n\times n\), \(\mathbf{x}\), \(n\times 1\), es la solución buscada y \(\mathbf{b}\), \(n\times 1\) es el término independiente.

En este caso, la función \(\mathbf{F}(\mathbf{x})\) sería lineal de la forma \(\mathbf{F}(\mathbf{x})=\mathbf{A}\mathbf{x}-\mathbf{b}\): \[ \left. \begin{align*} f_1(x_1,\ldots,x_n)=a_{11}x_1+a_{12}x_2+\cdots +a_{1n}x_n- b_1=&0,\\ f_2(x_1,\ldots,x_n)=a_{21}x_1+a_{22}x_2+\cdots +a_{2n}x_n- b_2=&0,\\ \vdots & \\ f_n(x_1,\ldots,x_n)=a_{n1}x_1+a_{n2}x_2+\cdots +a_{nn}x_n-b_n=&0. \end{align*} \right\} \]

Entonces, resolver la ecuación \(\mathbf{F}(\mathbf{x})=\mathbf{0}\) equivale a hallar aproximaciones de ceros \(\hat{\mathbf{x}}\) cuando la función \(\mathbf{F}(\mathbf{x})\) no es lineal, por ejemplo, si consideramos \(n=2\) y la función siguiente: \[ \mathbf{F}(\mathbf{x})=\mathbf{F}(x_1,x_2)=\begin{bmatrix}\frac{1}{5} \sin (x_1 x_2)-\frac{4 x_1}{5}+1 \\ \frac{1}{5}\mathrm{e}^{-x_1 x_2}-x_2 \end{bmatrix}. \] Nos preguntamos cómo hallar una solución aproximada \(\hat{\mathbf{x}}=(\hat{x}_1,\hat{x}_2)\) de la ecuación: \[ \mathbf{F}(\mathbf{x})=\mathbf{F}(x_1,x_2)=\begin{bmatrix}\frac{1}{5} \sin (x_1 x_2)-\frac{4 x_1}{5}+1 \\ \frac{1}{5}\mathrm{e}^{-x_1 x_2}-x_2 \end{bmatrix}=\begin{bmatrix}0\\ 0\end{bmatrix}. \]

El cálculo de ceros de funciones de diversas variables aparecen sobre todo en problemas de optimización donde se tienen que hallar extremos relativos (máximos o mínimos) de funciones \(f\) reales de \(n\) variables \(x_1,x_2,\ldots,x_n\), es decir, \[ f:\mathbf{D}\subseteq \mathbb{R}^n\longrightarrow \mathbb{R}, \] para cada \(\mathbf{x}=(x_1,x_2,\ldots,x_n)\in\mathbf{D}\), la función \(f\) nos da \(f(\mathbf{x})=f(x_1,x_2,\ldots,x_n)\in\mathbb{R}\).

En este caso, para hallar dichos extremos, tenemos que resolver \(\mathbf{Df}(\mathbf{x})=\mathbf{0}\), donde \(\mathbf{Df}(\mathbf{x})\) es la función que nos da las parciales respecto cada variable: \[ \mathbf{Df}(\mathbf{x})=\begin{bmatrix}\frac{\partial f}{\partial x_1}(\mathbf{x})\\ \frac{\partial f}{\partial x_2}(\mathbf{x})\\\vdots \\ \frac{\partial f}{\partial x_n}(\mathbf{x})\end{bmatrix}=\begin{bmatrix}0\\ 0\\\vdots\\ 0\end{bmatrix}. \]

En este caso, la función \(\mathbf{F}(\mathbf{x})\) de la cual tenemos que hallar un cero sería \(\mathbf{F}(\mathbf{x})=\mathbf{Df}(\mathbf{x})\).

Recordemos unas definiciones sobre límites y continuidad de funciones de varias variables:

Definición de límite.

Sea \(f\) una función real de varias variables, es decir, \(f:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}\). Sea \(\mathbf{x}_0\in\mathbf{D}\). Diremos que el límite de \(f\) en \(\mathbf{x}_0\) vale \(L\) y escribiremos, \[ \lim_{\mathbf{x}\to\mathbf{x}_0}f(\mathbf{x})=L, \] si, dado un valor \(\epsilon >0\), existe un valor \(\delta >0\) tal que si \(\mathbf{x}\in\mathbf{D}\) con \(\|\mathbf{x}-\mathbf{x}_0\|<\delta\), entonces \(|f(\mathbf{x})-L|<\epsilon\).

Observación.

Cuando escribimos \(\|\mathbf{x}-\mathbf{x}_0\|\), \(\|\cdot\|\) representa cualquier norma vectorial definida sobre \(\mathbb{R}^n\), como ya vimos en la parte de Álgebra lineal numérica del curso Métodos numéricos con Python. Álgebra Lineal Numérica.

Definición de continuidad en un punto.

Sea \(f\) una función real de varias variables, es decir, \(f:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}\). Sea \(\mathbf{x}_0\in\mathbf{D}\). Diremos que la función \(f\) es continua en \(\mathbf{x}_0\) si

• el límite de \(f(\mathbf{x})\) en \(\mathbf{x}_0\) existe, \(\displaystyle\exists\lim_{\mathbf{x}\to\mathbf{x}_0}f(\mathbf{x})\) y,
• \(\displaystyle\lim_{\mathbf{x}\to\mathbf{x}_0}f(\mathbf{x})=f(\mathbf{x}_0)\).

Diremos que \(f\) es continua en todo su dominio \(\mathbf{D}\) si es continua en todo punto \(\mathbf{x}_0\in\mathbf{D}\).

Las definiciones de límite de una función en un punto y continuidad se pueden extender a funciones \(\mathbf{F}:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}^n\).

Definición de límite.

Sea \(\mathbf{F}\) una función de varias variables, es decir, \(\mathbf{F}:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}^n\) donde escribimos: \[ \mathbf{F}(\mathbf{x})=\begin{bmatrix}f_1(x_1,\ldots,x_n)\\ f_2(x_1,\ldots,x_n)\\\vdots \\ f_n(x_1,\ldots,x_n)\end{bmatrix}. \]

Sea \(\mathbf{x}_0\in\mathbf{D}\). Diremos que el límite de \(\mathbf{F}\) en \(\mathbf{x}_0\) vale \(\mathbf{L}=\begin{bmatrix}L_1\\ \vdots\\ L_n\end{bmatrix}\) si, y sólo sí, el límite de cada componente \(f_i(\mathbf{x})\) en \(\mathbf{x}_0\) vale \(L_i\), para \(i=1,2,\ldots,n\).

Definición de continuidad en un punto.

Sea \(\mathbf{F}\) una función de varias variables, es decir, \(\mathbf{F}:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}^n\) donde escribimos: \[ \mathbf{F}(\mathbf{x})=\begin{bmatrix}f_1(x_1,\ldots,x_n)\\ f_2(x_1,\ldots,x_n)\\\vdots \\ f_n(x_1,\ldots,x_n)\end{bmatrix}. \]

Sea \(\mathbf{x}_0\in\mathbf{D}\). Diremos que \(\mathbf{F}\) es continua en \(\mathbf{x}_0\) si, y sólo sí, cada componente \(f_i(\mathbf{x})\) es continua en \(\mathbf{x}_0\), para \(i=1,2,\ldots,n\).

El siguiente resultado relaciona las derivadas parciales con la continuidad para funciones de varias variables sin entrar en detalle en el concepto de diferenciabilidad al ser éste un concepto difícil de manejar en la práctica.

Teorema: relación derivadas parciales y continuidad para funciones de varias variables.

Sea \(f\) una función real de varias variables, es decir, \(f:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}\). Sea \(\mathbf{x}_0\in\mathbf{D}\). Supongamos que todas las derivadas parciales \(\frac{\partial f}{\partial x_i}(\mathbf{x}_0)\) existen en \(\mathbf{x}_0\) que está localmente acotadas en \(\mathbf{x}_0\), esto es, existen constantes \(\delta >0\), \(K>0\) tal que para todo \(\mathbf{x}\in\mathbf{D}\) con \(\|\mathbf{x}-\mathbf{x}_0\|<\delta\) se cumple \[ \left|\frac{\partial f}{\partial x_i}(\mathbf{x}_0)\right|\leq K,\ i=1,2,\ldots,n. \] Entonces \(f\) es continua en \(\mathbf{x}_0\).

Uno de los métodos que está desarrollado en el Capítulo para hallar ceros de funciones reales de variable real del curso Métodos numéricos con Python. Cálculo Numérico es el método del punto fijo.

De forma concisa, dicho método consiste en transformar la ecuación \(f(x)=0\) de la cual queremos hallar una solución aproximada \(\hat{x}\) en otra equivalente de la forma \(x=g(x)\) y, a partir de esta ecuación, si la función \(g\) verifica que \(|g'(x)|<k<1\) en un entorno de la solución \(p\) de la ecuación \(f(x)=0\), entonces la sucesión \((x_k)_{k\geq 0}\) definida por \(x_{k+1}=g(x_k)\) converge hacia dicha solución \(p\).

Vamos a generalizar este método para funciones \(\mathbf{F}:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}^n\) de varias variables.

Primero, definamos el concepto de punto fijo:

Definición de punto fijo para funciones de varias variables.

Sea \(\mathbf{F}:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow\mathbb{R}^n\) una función de varias variables. Diremos que la función \(\mathbf{F}\) tiene un punto fijo \(\mathbf{p}\in\mathbf{D}\) si \(\mathbf{F}(\mathbf{p})=\mathbf{p}\).

Observación.

La condición \(\mathbf{F}(\mathbf{p})=\mathbf{p}\) son \(n\) condiciones que se tienen que cumplir si usamos las componentes de la función \(\mathbf{F}(\mathbf{x})=\begin{bmatrix}f_1(x_1,\ldots,x_n)\\ f_2(x_1,\ldots,x_n)\\\vdots \\ f_n(x_1,\ldots,x_n)\end{bmatrix}\),

\[ \mathbf{F}(\mathbf{p})=\begin{bmatrix}f_1(p_1,\ldots,p_n)\\ f_2(p_1,\ldots,p_n)\\\vdots \\ f_n(p_1,\ldots,p_n)\end{bmatrix}=\begin{bmatrix}p_1\\ p_2\\\vdots\\ p_n\end{bmatrix}=\mathbf{p},\ \Rightarrow \begin{align*} f_1(p_1,\ldots,p_n) = & p_1,\\ f_2(p_1,\ldots,p_n) = & p_2,\\ &\vdots\\ f_n(p_1,\ldots,p_n) = & p_n. \end{align*} \]

Supongamos que somos capaces de escribir nuestra ecuación a resolver \(\mathbf{F}(\mathbf{x})=\mathbf{0}\) de una forma equivalente a \(\mathbf{G}(\mathbf{x})=\mathbf{x}\).

Veamos bajo qué condiciones podemos afirmar que la función \(\mathbf{G}(\mathbf{x})\) tiene un punto fijo.

Teorema del punto fijo.

Sea \(\mathbf{D}=\{\mathbf{x}=(x_1,\ldots,x_n)\in\mathbb{R}^n,\ a_i\leq x_i\leq b_i,\ i=1,\ldots,n\}\) para un conjunto de valores \(a_1,\ldots,a_n\), \(b_1,\ldots,b_n\). Sea \(\mathbf{G}(\mathbf{x})\) una función de varias variables continua definida sobre \(\mathbf{D}\), \(\mathbf{G}:\mathbf{D}\subset\mathbb{R}^n\longrightarrow\mathbb{R}^n\) tal que \(\mathbf{G}(\mathbf{x})\in\mathbf{D}\) para todo valor \(\mathbf{x}\in \mathbf{D}\). Entonces, \(\mathbf{G}(\mathbf{x})\) tiene un punto fijo en \(\mathbf{D}\), es decir, existe un valor \(\mathbf{p}\in\mathbf{D}\) tal que \(\mathbf{G}(\mathbf{p})=\mathbf{p}\).

Observación.

El Teorema anterior nos dice que si una función \(\mathbf{G}(\mathbf{x})\) definida sobre un rectángulo \(n\)-dimensional cumple que la imagen de todo punto del rectángulo pertenece a dicho rectángulo, entonces la función \(\mathbf{G}(\mathbf{x})\) tiene un punto fijo en dicho rectángulo.

La idea es que si \(\mathbf{p}\) es un punto fijo de la función de varias variables \(\mathbf{G}(\mathbf{x})\), entonces \(\mathbf{p}\) necesariamente ha de ser un cero de nuestra función de varias variables \(\mathbf{F}(\mathbf{x})\).

Vamos a ver qué condiciones debe cumplir una función de varias variables \(\mathbf{G}(\mathbf{x})\) para que la sucesión en varias variables \((\mathbf{x}^{(k)})_{k\geq 0}\) definida de forma recurrente como \(\mathbf{x}^{(k)}=\mathbf{G}(\mathbf{x}^{(k-1)})\), con \(\mathbf{x}^{(0)}\) en el dominio de definición de dicha función \(\mathbf{G}(\mathbf{x})\) converja hacia un punto fijo de la función \(\mathbf{G}(\mathbf{x})\) y por tanto, hacia un cero de la función \(\mathbf{F}(\mathbf{x})\):

Teorema. Condición de convergencia.

Sea \(\mathbf{D}=\{\mathbf{x}=(x_1,\ldots,x_n)\in\mathbb{R}^n,\ a_i\leq x_i\leq b_i,\ i=1,\ldots,n\}\) para un conjunto de valores \(a_1,\ldots,a_n\), \(b_1,\ldots,b_n\). Sea \(\mathbf{G}(\mathbf{x})\) una función de varias variables continua definida sobre \(\mathbf{D}\), \(\mathbf{G}:\mathbf{D}\subset\mathbb{R}^n\longrightarrow\mathbb{R}^n\) tal que \(\mathbf{G}(\mathbf{x})\in\mathbf{D}\) para todo valor \(\mathbf{x}\in \mathbf{D}\).

Sean \(g_1,\ldots,g_n\) las componentes de la función \(\mathbf{G}\), \(\mathbf{G}(\mathbf{x})=(g_1(\mathbf{x}),\ldots,g_n(\mathbf{x}))^\top\).

Suponemos además que \(\mathbf{G}\) es diferenciable y que existe una constante \(C<1\) tal que: \[ \left|\frac{\partial g_i(\mathbf{x})}{\partial x_j}\right|\leq \frac{C}{n}, \] para todo \(i,j=1,\ldots,n\) y para todo \(\mathbf{x}\in\mathbf{D}\).

Teorema. Condición de convergencia (continuación).

Entonces la sucesión en \(\mathbb{R}^n\) definida por: \[ \mathbf{x}^{(k)}=\mathbf{G}(\mathbf{x}^{(k-1)}),\ \mbox{ con }\mathbf{x}^{(0)}\in\mathbf{D}, \] converge hacia el único punto fijo \(\mathbf{p}\) de la función \(\mathbf{G}(\mathbf{x})\) y además, \[ \|\mathbf{x}^{(k)}-\mathbf{p}\|\leq \frac{C^k}{1-C}\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\|. \]

• INPUT \(n\), \(g_1(x_1,\ldots,x_n),\ldots,g_n(x_1,\ldots,x_n)\),componentes de la función de iteración \(\mathbf{G}(\mathbf{x})\), \(\mathbf{x}^{(0)}\) (valor inicial), \(TOL\) (tolerancia), \(Nmax\) (número máximo de iteraciones)
• Set \(k=1\) (inicializamos el contador de las iteraciones)
• While \(k\leq Nmax\)
  • For i=1,...,n (calculamos el siguiente valor \(\mathbf{x}\) de la sucesión)
    • Set \(x_i=g_i(\mathbf{x}^{0})\)
  • Set \(ERROR = \|\mathbf{x}-\mathbf{x}^{(0)}\|_2\) (calculamos el error entre los dos términos consecutivos de la sucesión)

• • If \(ERROR < TOL\)
    • Print La solución aproximada vale \(\mathbf{x}\).
    • STOP
  • Set \(\mathbf{x}^{(0)}=\mathbf{x}\) (actualizamos los términos de la sucesión)
  • Set \(k=k+1\)
• Print El método no converge
• STOP

Para acelerar la convergencia de la sucesión podemos usar una metodología similar al método de Gauss-Seidel introducido en el capítulo de métodos iterativos del curso Métodos numéricos con Python. Álgebra Lineal Numérica.

Concretamente, para calcular la componente \(x_i^{(k)}\), podemos usar las componentes calculadas \(x_j^{(k)}\), \(j=1,\ldots,i-1\). Es decir la recurrencia \(\mathbf{x}^{(k)}=\mathbf{G}(\mathbf{x}^{(k-1)})\) sería la siguiente: \[ \begin{align*} x_1^{(k)}= & g_1(x_1^{(k-1)},\ldots,x_n^{(k-1)}),\\ x_2^{(k)}= & g_2(x_1^{(k)},x_2^{(k-1)}\ldots,x_n^{(k-1)}),\\ &\vdots\\ x_i^{(k)}=& g_i(x_1^{(k)},\ldots,x_{i-1}^{(k)},x_i^{(k-1)},\ldots,x_n^{(k-1)}),\\ &\vdots\\ x_n^{(k)}= & g_n(x_1^{(k)},\ldots,x_n^{(k-1)}). \end{align*} \]

Ahora bien, no siempre se acelera la convergencia usando la metodología anterior.

Recordemos el método de Newton-Raphson introducido en el capítulo de ceros de funciones del curso Métodos numéricos con Python. Cálculo Numérico:

Queremos hallar una aproximación de la ecuación \(f(x)=0\), donde \(f\) es derivable. Para ello definimos la sucesión \((x_k)_k\) de forma recurrente de la forma siguiente: \[ x_{k}=x_{k-1}-\frac{f(x_{k-1})}{f'(x_{k-1})}, \] dando un valor \(x_0\) inicial.

En el curso anterior se estudian cuestiones referentes a la convergencia y a la velocidad de convergencia.

En primer lugar, observemos que el método de Newton-Raphson puede considerarse un caso particular de aplicar el método del punto fijo usando como función de iteración \(g(x)=x-\frac{f(x)}{f'(x)}\).

De hecho, y de cara a generalizar dicho método a hallar ceros de funciones de varias variables \(\mathbf{F}(\mathbf{x})=\mathbf{0}\), podemos pensar que el método de Newton-Raphson es aplicar el método del punto fijo cuando la función de iteración tiene la forma siguiente \(g(x)=x-\varphi(x)f(x)\), donde en el caso del método de Newton-Raphson, la función \(\varphi(x)=\frac{1}{f'(x)}\).

La generalización natural a las funciones de diversas variables es aplicar el método del punto fijo usando la función de iteración siguiente: \(\mathbf{G}({\mathbf{x}})=\mathbf{x}-\mathbf{A}(\mathbf{x})^{-1}\mathbf{F}(\mathbf{x})\), donde \(\mathbf{A}(\mathbf{x})\) sería la matriz Jacobiana o de las derivadas parciales de la función \(\mathbf{F}(x)\) de la que queremos hallar un cero.

Como el método de Newton-Raphson es un caso particular de aplicación del método del punto fijo, antes de entrar en detalles en el método de Newton-Raphson, veamos el resultado siguiente que nos da condiciones de convergencia para el método del punto fijo:

Teorema: convergencia del método del punto fijo.

Sea \(\mathbf{p}\) un punto fijo o una solución de la ecuación \(\mathbf{G}(\mathbf{x})=\mathbf{x}\). Supongamos que existe un número \(\delta >0\) que cumple lo siguiente:

1. \(\frac{\partial g_i}{\partial x_j}\) es una función continua en el dominio \(N_\delta = \{\mathbf{x},\ |\ \|\mathbf{x}-\mathbf{p}\| <\delta\}\) para \(i,j=1,\ldots,n\).
2. \(\frac{\partial^2 g_i}{\partial x_j\partial x_k}\) es una función continua y además existe una constante \(M\) tal que \(\left|\frac{\partial^2 g_i}{\partial x_j\partial x_k}\right|\leq M,\) si \(\mathbf{x}\in N_\delta\), para \(i,j,k=1,\ldots,n\).
3. \(\frac{\partial g_i(\mathbf{p})}{\partial x_j}=0,\) para \(i,j=1,\ldots,n\).

Teorema: convergencia del método del punto fijo (continuación).

Entonces existe un valor \(0<\delta'\leq \delta\) tal que la sucesión generada por \(\mathbf{x}^{(k)}=\mathbf{G}(\mathbf{x}^{(k-1)})\) converge con velocidad cuadrática al punto fijo \(\mathbf{p}\) para cualquier valor inicial \(\mathbf{x}^{(0)}\). Además, si \(\|\mathbf{x}^{(0)}-\mathbf{p}\|<\delta'\) se cumple lo siguiente: \[ \|\mathbf{x}^{(k)}-\mathbf{p}\|_\infty\leq \frac{n^2M}{2}\|\mathbf{x}^{(k-1)}-\mathbf{p}\|_\infty^2,\ k\geq 1. \]

Observación.

La velocidad de convergencia es cuadrática ya que usando la desigualdad anterior tenemos que: \[ \frac{\|\mathbf{x}^{(k)}-\mathbf{p}\|_\infty}{\|\mathbf{x}^{(k-1)}-\mathbf{p}\|_\infty^2}\leq \frac{n^2M}{2}, \] lo que significa que \[ \lim_{k\to\infty} \frac{\|\mathbf{x}^{(k)}-\mathbf{p}\|_\infty}{\|\mathbf{x}^{(k-1)}-\mathbf{p}\|_\infty^2}=\mbox{constante.} \]

Observación.

La demostración del Teorema anterior se basa en el desarrollo de Taylor de la función de varias variables \(\mathbf{G}(x)\) alrededor del punto \(\mathbf{p}\).

Usando las condiciones del Teorema, no es complicado demostrar su tesis. Sin embargo, la notación es muy técnica, creemos que no aporta ningún avance en la comprensión del tema y por dicho motivo, se omite.

A continuación, si escribimos \(\mathbf{G}(x)=\mathbf{x}-\mathbf{A}(\mathbf{x})^{-1}\mathbf{F}(\mathbf{x})\) veamos qué condiciones debe verificar la matriz-función de varias variables \(\mathbf{A}(\mathbf{x})\) para verificar las condiciones del teorema anterior.

Sean \(b_{ij}(\mathbf{x})\) la componente \((i,j)\) de la matriz \(\mathbf{A}(\mathbf{x})^{-1}\), es decir:

\[ \mathbf{A}(\mathbf{x})^{-1}=\begin{bmatrix}b_{11}(\mathbf{x})&b_{12}(\mathbf{x})&\ldots&b_{1n}(\mathbf{x}) \\b_{21}(\mathbf{x})&b_{22}(\mathbf{x})&\ldots&b_{2n}(\mathbf{x}) \\\vdots&\vdots&\ddots&\vdots \\b_{n1}(\mathbf{x})&b_{n2}(\mathbf{x})&\ldots&b_{nn}(\mathbf{x}) \\\end{bmatrix}. \]

Entonces la componente \(g_i(\mathbf{x})\) de la función \(\mathbf{G}(\mathbf{x})\) será: \[ g_i(\mathbf{x})=x_i-\sum_{k=1}^n b_{ik}(\mathbf{x})f_k(\mathbf{x}),\ i=1,\ldots,n, \] siendo \(f_k(\mathbf{x})\) la componente \(k\)-ésima de la función \(\mathbf{F}(\mathbf{x})\).

La condición 3. del Teorema nos dice que \(\frac{\partial g_i(\mathbf{p})}{\partial x_j}=0\), \(i,j=1,\ldots,n\). Hallemos, pues, \(\frac{\partial g_i(\mathbf{x})}{\partial x_j}\) en un valor cualquiera \(\mathbf{x}\): \[ \frac{\partial g_i(\mathbf{x})}{\partial x_j}=\begin{cases}1-\sum_{k=1}^n \left(b_{ik}(\mathbf{x})\frac{\partial f_k(\mathbf{x})}{\partial x_j}+\frac{\partial b_{ik}(\mathbf{x})}{\partial x_j}f_k(\mathbf{x})\right),& \mbox{ si } i=j,\\ -\sum_{k=1}^n \left(b_{ik}(\mathbf{x})\frac{\partial f_k(\mathbf{x})}{\partial x_j}+\frac{\partial b_{ik}(\mathbf{x})}{\partial x_j}f_k(\mathbf{x})\right),& \mbox{ si } i\neq j, \end{cases} \]

Como el punto fijo de la función \(\mathbf{G}(\mathbf{x})\), \(\mathbf{G}(\mathbf{p})=\mathbf{p}\) debe ser un cero de la función \(\mathbf{F}(\mathbf{x})\), tenemos que \(\mathbf{F}(\mathbf{p})=0\) o, escrito en componentes \(f_k(\mathbf{p})=0\), \(k=1,\ldots,n\).

Si calculamos \(\frac{\partial g_i(\mathbf{p})}{\partial x_j}\) usando la expresión anterior obtenemos: \[ \frac{\partial g_i(\mathbf{p})}{\partial x_j}=\begin{cases}1-\sum_{k=1}^n b_{ik}(\mathbf{p})\frac{\partial f_k(\mathbf{p})}{\partial x_j},& \mbox{ si } i=j,\\ -\sum_{k=1}^n b_{ik}(\mathbf{p})\frac{\partial f_k(\mathbf{x})}{\partial x_j},& \mbox{ si } i\neq j, \end{cases} \] Como \(\frac{\partial g_i(\mathbf{p})}{\partial x_j}=0\), \(i,j=1,\ldots,n\), podemos escribir matricialmente que: \[ \mathbf{A}^{-1}(\mathbf{p})\mathbf{J}(\mathbf{p})=\mathbf{Id}, \]

donde \(\mathbf{J}(\mathbf{x})\) es la matriz jacobiana de la función \(\mathbf{F}(\mathbf{x})\):

\[ \mathbf{J}(\mathbf{x})=\begin{bmatrix}\frac{\partial f_1(\mathbf{x})}{\partial x_1}&\frac{\partial f_1(\mathbf{x})}{\partial x_2}&\ldots&\frac{\partial f_1(\mathbf{x})}{\partial x_n} \\\frac{\partial f_2(\mathbf{x})}{\partial x_1}&\frac{\partial f_2(\mathbf{x})}{\partial x_2}&\ldots&\frac{\partial f_2(\mathbf{x})}{\partial x_n} \\\vdots&\vdots&\ddots&\vdots \\\frac{\partial f_n(\mathbf{x})}{\partial x_1}&\frac{\partial f_n(\mathbf{x})}{\partial x_2}&\ldots&\frac{\partial f_n(\mathbf{x})}{\partial x_n} \\\end{bmatrix}. \]

Entonces, como \(\mathbf{A}^{-1}(\mathbf{p})\mathbf{J}(\mathbf{p})=\mathbf{Id},\) tenemos que \(\mathbf{A}(\mathbf{p})=\mathbf{J}(\mathbf{p})\) y la función de iteración será, pues: \[ \mathbf{G}(\mathbf{x})=\mathbf{x}-\mathbf{J}(\mathbf{x})^{-1}\mathbf{F}({\mathbf{x}}), \] y la sucesión \(\mathbf{x}^{(k)}\) queda definida de la forma siguiente: \[ \mathbf{x}^{(k)}=\mathbf{G}(\mathbf{x}^{(k-1)})=\mathbf{x}^{(k-1)}-\mathbf{J}(\mathbf{x}^{(k-1)})^{-1}\mathbf{F}({\mathbf{x}^{(k-1)}}), \] eligiendo un valor inicial \(\mathbf{x}^{(0)}\).

El método anterior se denomina método de Newton-Raphson para resolver sistemas de ecuaciones no lineales \(\mathbf{F}(\mathbf{x})=\mathbf{0}\).

La parte positiva del método anterior es que se espera que tenga convergencia cuadrática en caso que se cumplan las condiciones del Teorema anterior.

La parte negativa es que se tiene que calcular una inversa de una matriz \(n\times n\), \(\mathbf{J}(\mathbf{x})\) en cada iteración.

En la práctica, no se calcula ninguna inversa sino que en cada paso se resuelve el siguiente sistema de ecuaciones lineal: \[ \begin{align*} \mathbf{J}(\mathbf{x}^{(k-1)})\mathbf{y}^{(k)}= & -\mathbf{F}({\mathbf{x}^{(k-1)}}),\ \Rightarrow \mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+\mathbf{y}^{(k)}. \end{align*} \] Es decir, en el paso \(k\)-ésimo se tiene que resolver un sistema lineal de la forma \(\mathbf{A}\mathbf{x}=\mathbf{b}\), donde la matriz del sistema \(\mathbf{A}\) vale \(\mathbf{A}=\mathbf{J}(\mathbf{x}^{(k-1)})\) y el término independiente \(\mathbf{b}\) vale \(\mathbf{b}=-\mathbf{F}({\mathbf{x}^{(k-1)}})\) y cuya solución hemos denominado \(\mathbf{y}^{(k)}\).

Notemos como en cada paso de la sucesión, la matriz del sistema cambia y debe ser calculada nuevamente antes de resolver numéricamente el sistema correspondiente.

Podemos usar un método directo o un método iterativo para resolver dicho sistema. Ver los temas de Métodos directos y Métodos iterativos del curso Métodos numéricos con Python. Álgebra Lineal Numérica para más detalles de cuando usar un método u otro.

El término siguiente de la sucesión será \(\mathbf{x}^{(k)}=\mathbf{x}^{(k-1)}+\mathbf{y}^{(k)}\).

• INPUT \(n\) número de ecuaciones, \(f_1(x_1,\ldots,x_n),\ldots,f_n(x_1,\ldots,x_n)\),componentes de la función de iteración \(\mathbf{F}(\mathbf{x})\), \(\mathbf{J}(\mathbf{x})\), cálculo de la matriz jacobiana en un punto \(\mathbf{x}\), \(\mathbf{x}^{(0)}\) aproximación inicial, \(Nmax\) (número máximo de iteraciones)
• Set \(k=1\) (contador para el número de iteraciones)
• While \(k\leq Nmax\)
  • Compute \(\mathbf{F}(\mathbf{x}^{(0)})\) y \(\mathbf{J}(\mathbf{x}^{(0)})\) (calculamos todas las componentes de la función \(\mathbf{F}(\mathbf{x})\) en \(\mathbf{x}^{(0)}\) y la matriz jacobiana \(\mathbf{J}(\mathbf{x}^{(0)})\) en \(\mathbf{x}^{(0)}\).
  • Solve \(\mathbf{J}(\mathbf{x}^{(0)})(\mathbf{y})=-\mathbf{F}(\mathbf{x}^{(0)})\) (resolvemos el sistema lineal \(\mathbf{J}(\mathbf{x}^{(0)})\mathbf{y}=-\mathbf{F}(\mathbf{x}^{(0)})\))

• • Set \(\mathbf{x}^{(0)}=\mathbf{x}^{(0)}+\mathbf{y}\) (calculamos el siguiente término de la sucesión)
  • If \(\|\mathbf{y}\| < TOL\) (como el error entre dos términos consecutivos vale la solución del sistema de ecuaciones \(\|\mathbf{y}\|\), miramos si dicho error es menor que la tolerancia)
    • Print La solución aproximada vale \(\mathbf{x}^{(0)}\).
    • STOP.
  • Set \(k=k+1\) (aumentamos el contador de iteraciones en uno)
• Print El método no converge
• STOP.

El método de Newton-Raphson, como hemos comentado

• tiene como ventaja la velocidad de convergencia cuadrática pero
• tiene como inconveniente, aparte que se tiene que resolver un sistema lineal en cada iteración, que no tenemos asegurada la convergencia en el sentido de que si el valor inicial no está suficientemente cercano a la solución, puede que el método no converja.

En esta sección, vamos a introducir el método del gradiente descendente que sólo tiene convergencia lineal pero usualmente converge aunque el valor inicial sea una aproximación pobre de la solución buscada.

Dicho método está basado en el resultado siguiente:

Proposición: conexión entre el cero de funciones de varias variables y el mínimo de una función real de varias variables.

Sea \(\mathbf{F}(\mathbf{x})\) una función de varias variables, \(\mathbf{F}:\mathbf{D}\subseteq\mathbb{R}^n\longrightarrow \mathbb{R}^n\) de componentes \(f_1(\mathbf{x}),\ldots,f_n(\mathbf{x})\).

Sea \({\mathbf{p}}=(p_1,\ldots,p_n)\) una solución del sistema de ecuaciones no lineal \(\mathbf{F}(\mathbf{x})=\mathbf{0}\): \[ \left. \begin{align*} f_1(x_1,\ldots,x_n) = & 0,\\ f_2(x_1,\ldots,x_n) = & 0,\\ &\vdots\\ f_n(x_1,\ldots,x_n) = & 0. \end{align*} \right\} \]

Proposición: conexión entre el cero de funciones de varias variables y el mínimo de una función real de varias variables. (continuación)

Entonces \(\mathbf{p}\) es un mínimo de la siguiente función real de varias variables: \[ g(\mathbf{x})=g(x_1,\ldots,x_n)=\sum_{i=1}^n (f_i(x_1,\ldots,x_n))^2, \] y además \(g(\mathbf{p})=0\).

El resultado anterior nos transforma el problema de hallar un cero de una función de diversas variables \(\mathbf{F}(\mathbf{x})\) en un problema de hallar un mínimo de una función \(\mathbf{g}(\mathbf{x})\).

Para hallar dicho mínimo, vamos a obtener una sucesión de aproximaciones \((\mathbf{x}^{(k)})_{k\geq 0}\) de la forma siguiente:

1. Se considera una aproximación inicial \(\mathbf{x}^{(0)}=(x_1^{(0)},\ldots,x_n^{(0)})^\top\).
2. Se halla la dirección \(\mathbf{v}\) que maximiza el decrecimiento de la función \(g(\mathbf{x})\) a partir de la aproximación inicial \(\mathbf{x}^{(0)}\).
3. Nos movemos hacia la dirección hallada en el paso anterior para hallar la nueva aproximación \(\mathbf{x}^{(1)}\).
4. Volvemos al primer paso reemplazando el valor de \(\mathbf{x}^{(0)}\) por \(\mathbf{x}^{(1)}\) hasta que la diferencia \(\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\|_2\) sea menor que una cierta tolerancia permitida.

Los pasos 2. y 3. han quedado un poco difusos y se tienen que concretar.

En el paso 2. se habla de hallar la dirección \(\mathbf{v}\) que maximiza el decrecimiento de la función \(g(\mathbf{x})\) a partir de la aproximación inicial \(\mathbf{x}^{(0)}\).

Es decir, dada una dirección \(\mathbf{v}\), que suponemos unitaria o de norma euclídea unidad, (\(\displaystyle\|\mathbf{v}\|_2^2=\sum_{i=1}^n v_i^2=1\)) consideramos la función real de una variable \(h_{\mathbf{v}}(t)=g(\mathbf{x}^{(0)}+t\mathbf{v})\) que dice cómo se mueve la función \(g(\mathbf{x})\) en dicha dirección.

Fijémonos que \(h_{\mathbf{v}}(0)=g(\mathbf{x}^{(0)})\). La idea es hallar la dirección \(\mathbf{v}\) tal que \(h'_{\mathbf{v}}(0)\) sea mínima ya que la derivada de la función \(h_{\mathbf{v}}(t)\) nos da la tasa de crecimiento o decrecimiento de dicha función.

El valor \(h'_{\mathbf{v}}(0)\) se denomina derivada direccional de la función \(g(\mathbf{x})\) en \(\mathbf{x}^{(0)}\) en la dirección \(\mathbf{v}\): \[ h'_{\mathbf{v}}(0)=D_{\mathbf{v}}(g)(\mathbf{x}^{(0)})=\lim_{u\to 0}\frac{h_{\mathbf{v}}(u)-h_{\mathbf{v}}(0)}{u}=\lim_{u\to 0}\frac{g(\mathbf{x}^{(0)}+u\mathbf{v})-g(\mathbf{x}^{(0)})}{u}. \]

Para calcular \(h'_{\mathbf{v}}(t)\) podemos usar la regla de la cadena: \[ \begin{align*} h'_{\mathbf{v}}(t) = & \nabla g(\mathbf{x}^{(0)}+t\mathbf{v})\cdot\mathbf{v}=\left(\frac{\partial g}{\partial x_1}(\mathbf{x}^{(0)}+t\mathbf{v}),\ldots,\frac{\partial g}{\partial x_n}(\mathbf{x}^{(0)}+t\mathbf{v})\right)\cdot\mathbf{v} \\ = & \sum_{i=1}^n \frac{\partial g}{\partial x_i}(\mathbf{x}^{(0)}+t\mathbf{v})\cdot v_i. \end{align*} \]

Usando la expresión anterior, podemos hallar la derivada direccional de la función \(g(\mathbf{x})\) en \(\mathbf{x}^{(0)}\) en la dirección \(\mathbf{v}\): \[ h'_{\mathbf{v}}(0)=\nabla g(\mathbf{x}^{(0)})\cdot\mathbf{v}=\sum_{i=1}^n \frac{\partial g}{\partial x_i}(\mathbf{x}^{(0)})\cdot v_i. \]

Fijémonos que \(h'_{\mathbf{v}}\) es el producto escalar del vector gradiente de la función \(g(\mathbf{x})\) por la dirección \(\mathbf{v}\). (\(h'_{\mathbf{v}}(0)=\nabla g(\mathbf{x}^{(0)})\cdot\mathbf{v}\))

Ahora, ¿cuál es la dirección \(\mathbf{v}\) que minimiza el producto escalar anterior \(\nabla g(\mathbf{x}^{(0)})\cdot\mathbf{v}\)?

En general, dado un producto escalar de dos vectores, en nuestro caso \(\nabla g(\mathbf{x}^{(0)})\cdot\mathbf{v}\) donde el primero de ellos es fijo \(\nabla g(\mathbf{x}^{(0)})\), dicho producto escalar se minimiza cuando el vector \(\mathbf{v}\) se escoge en la dirección opuesta de la del vector fijo \(\nabla g(\mathbf{x}^{(0)})\) ya que recordemos que el producto escalar de dos vectores vale: \[ \nabla g(\mathbf{x}^{(0)})\cdot\mathbf{v}=\|\nabla g(\mathbf{x}^{(0)})\|_2\cdot\|\mathbf{v}\|_2\cdot\cos(\alpha)=\|\nabla g(\mathbf{x}^{(0)})\|_2\cdot\cos(\alpha), \] donde \(\alpha\) es el ángulo que forman. Entonces, para minimizar \(\cos(\alpha)\), tenemos que escoger \(\alpha=180^o=\pi\mathrm{\ rad.}\) que equivale a escoger \(\mathbf{v}=-\frac{\nabla g(\mathbf{x}^{(0)})}{\|\nabla g(\mathbf{x}^{(0)})\|_2}\).

En resumen, ya hemos “resuelto” el paso 2: tenemos que escoger como dirección \(\mathbf{v}\) la opuesta de la del gradiente de la función \(g(\mathbf{x})\) en \(\mathbf{x}^{(0)}\): \(\mathbf{v}=-\frac{\nabla g(\mathbf{x}^{(0)})}{\|\nabla g(\mathbf{x}^{(0)})\|_2}\).

Vamos a por el tercer paso.

Tenemos que movernos de la forma siguiente una cantidad \(u>0\) a determinar tal que la función siguiente: \[ \mathbf{x}^{(1)}=h_{\mathbf{v}}(u)=g\left(\mathbf{x}^{(0)}-u\frac{\nabla g(\mathbf{x}^{(0)})}{\|\nabla g(\mathbf{x}^{(0)})\|_2}\right), \] donde \(\mathbf{v}=-\frac{\nabla g(\mathbf{x}^{(0)})}{\|\nabla g(\mathbf{x}^{(0)})\|_2}\), tenga un minimo en \(u>0\).

Hallar el valor \(u>0\) que minimize la función anterior \(h_{\mathbf{v}}(u)\) es complicado ya que tenemos que hallar un cero de la siguiente función no lineal: \[ h'_{\mathbf{v}}(u)=\nabla g(\mathbf{x}^{(0)}-u\nabla g(\mathbf{x}^{(0)}))\cdot \left(-\nabla g(\mathbf{x}^{(0)}\right)=0. \]

El método del gradiente descendente simplifica el paso anterior considerando tres valores \(u_1 <u_2<u_3\) que esperamos que estén cerca del mínimo a calcular, hallamos la función cuadrática \(P_2(u)\) que interpola la función \(h_{\mathbf{v}}(u)\) en los tres valores anteriores y consideramos que el valor mínimo a computar es el mínimo de la función cuadrática \(P_2(u)\).

Hallar el mínimo de una función cuadrática es una tarea relativamente sencilla de calcular.

Sea \(\hat{u}\) el mínimo de la función cuadrática \(P_2(u)\). El nuevo valor de la sucesión \(\mathbf{x}^{(1)}\) será pues: \[ \mathbf{x}^{(1)}=h_{\mathbf{v}}(\hat{u})=g(\mathbf{x}^{(0)}-\hat{u}\nabla g(\mathbf{x}^{(0)})). \]

En la práctica, se escoge \(u_1 =0\) ya que \(h_{\mathbf{v}}(u_1=0)=g(\mathbf{x}^{(0)})\).

Como \(u_1\) no es el mínimo de la función \(h_{\mathbf{v}}(u)\), existe un número \(u_3>u_1=0\) tal que \(h_{\mathbf{v}}(u_3)<h_{\mathbf{v}}(u_1=0)\). Por último, escogemos \(u_2=\frac{u_3}{2}\).

El mínimo de la función cuadrática \(P_2(u)\) en el intervalo cerrado \([u_1=0,u_3]\) está asegurado ya que toda función cuadrática alcanza un mínimo en un intervalo cerrado, o en el interior o en los extremos.

Hallemos a continuación una expresión del mínimo de la función cuadrática \(P_2(u)\).

Sean \(g_1=h_{\mathbf{v}}(0)=g(x^{(0)})\), \(g_2=h_{\mathbf{v}}(u_2)=g(\mathbf{x}^{(0)}-u_2\nabla g(\mathbf{x}^{(0)}))\), \(g_3=h_{\mathbf{v}}(u_3)=g(\mathbf{x}^{(0)}-u_3\nabla g(\mathbf{x}^{(0)}))\). Para hallar el polinomio \(P_2(u)\) que interpola los puntos \((0,g_1), (u_2,g_2)\) y \((u_3,g_3)\) podemos usar el método de las diferencias divididas de Newton (ver el tema de interpolación del curso Métodos Numéricos con Python. Cálculo Numérico):

El polinomio de interpolación \(P_2(u)\) será, pues, \[ P_2(u)=g_1+h_1 u+h_3 u (u-u_2). \] Para hallar el mínimo del polinomio anterior, hemos de derivar e igualar a cero: \[ P_2'(u)=h_1+h_3(2u-u_2)=0,\ \Rightarrow u =\frac{1}{2}\left(u_2-\frac{h_1}{h_3}\right). \]

Como para aplicar el método del gradiente descendente, necesitamos conocer el gradiente de la función \(g(\mathbf{x})\), \(\nabla g(\mathbf{x})=\left(\frac{\partial g(\mathbf{x})}{\partial x_1},\ldots,\frac{\partial g(\mathbf{x})}{\partial x_n}\right)^\top\), vamos a dar una expresión de dicho gradiente en función de la matriz jacobiana de la función \(\mathbf{F}(\mathbf{x})\).

Recordemos la expresión de la función \(g(\mathbf{x})\): \[ g(\mathbf{x})=g(x_1,\ldots,x_n)=\sum_{i=1}^n (f_i(x_1,\ldots,x_n))^2. \] Entonces, \[ \frac{\partial g(\mathbf{x})}{\partial x_k}=2\sum_{i=1}^n f_i(x_1,\ldots,x_n)\frac{\partial f_i(\mathbf{x})}{\partial x_k},\ k=1,\ldots,n. \]

Matricialmente, la expresión anterior se traduce en: \[ \nabla g(\mathbf{x})=2\mathbf{J}(\mathbf{x})^\top\cdot\mathbf{F}(\mathbf{x}), \] donde \(\mathbf{J}(\mathbf{x})\) sería la matriz jacobiana de la función \(\mathbf{F}(\mathbf{x})\) en \(\mathbf{x}\) y \(\mathbf{F}(\mathbf{x})=(f_1(\mathbf{x}),\ldots,f_n(\mathbf{x}))^\top\).

Suponemos que tenemos definidas rutinas que dan las componentes \(f_i(\mathbf{x})\), \(i=1,\ldots,n\) y las componentes de la matriz jacobiana de la función \(\mathbf{F}(\mathbf{x})\): \(\frac{\partial f_i(\mathbf{x})}{\partial x_k}\), \(i,k=1,\ldots,n\)

En primer lugar definimos la función \(g\) de la cual queremos hallar un mínimo:

• g=rutina(x1,...,xn)
  • Set \(y=0\) (inicializamos el valor de \(g(\mathbf{x}\)))
  • For i=1,...,n (calculamos el valor \(g(\mathbf{x})\))
    • Set \(y=y+f_i(x_1,\ldots,x_n)^2\)
  • Return \(y\).

En segundo lugar definimos la función \(\nabla g\) o el gradiente de \(g\) en un punto \(\mathbf{x}\):

• grad_g = rutina(x1,...,xn)
  • Set \(grad=\mathbf{0}\) (inicializamos el vector gradiente a \(\mathbf{0}\))
  • For k=1,...,n
    • For i=1,...,n
      • Set \(grad[k]=grad[k]+\frac{\partial f_i(x_1,\ldots,x_n)}{\partial x_k}\cdot f_i(x_1,\ldots,x_n)\)
  • Return \(2\cdot grad\).

Programa principal:

• INPUT \(n\), \(f_1(x_1,\ldots,x_n),\ldots,f_n(x_1,\ldots,x_n)\),componentes de la función \(\mathbf{F}(\mathbf{x})\), \(\mathbf{x}^{(0)}\) (valor inicial), \(TOL\) (tolerancia), \(Nmax\) (número máximo de iteraciones)
• Set \(k=1\) (inicializamos el contador de las iteraciones)
• While \(k\leq Nmax\)
  • Set \(g_1=g(x_1^{(0)},\ldots,x_n^{(0)})\).
  • Set \(gradx = grad_g(x_1^{(0)},\ldots,x_n^{(0)})\).
  • Set \(normx = \|gradx\|_2\).

• • If \(normx=0\)
    • Print: hemos encontrado un valor de gradiente cero. Comprobar si \(\mathbf{x}^{(0)}\) es una aproximación de un cero de la función \(\mathbf{F}(\mathbf{x})\)
    • Print \(\mathbf{x}^{(0)}\)
    • STOP.
  • Set \(gradx=\frac{gradx}{normx}\) (normalizamos el gradiente)
  • Set \(u_1=0\) (inicializamos las \(u_i\))
  • Set \(u_3=1\)
  • Set \(g_3=g(\mathbf{x}^{(0)}-u_3\cdot gradx)\)

• • While \(g_3\geq g_1\) (refinamos el valor \(u_3\))
    • Set \(u_3=\frac{u_3}{2}\)
    • Set \(g_3=g(\mathbf{x}^{(0)}-u_3\cdot gradx)\)
    • If \(u_3 <\frac{TOL}{2}\)
      • Print: no hay mejora significativa. Comprobar si \(\mathbf{x}^{(0)}\) es una aproximación de un cero de la función \(\mathbf{F}(\mathbf{x})\)
    • Set \(u_2=\frac{u_3}{2}\).
    • Set \(g_2=g(\mathbf{x}^{(0)}-u_2\cdot gradx)\).

• • Set \(h_1=\frac{g_2-g_1}{u_2}\). (hallamos el mínimo del polinomio de interpolación \(P_2(u)\))
  • Set \(h_2=\frac{g_3-g_2}{u_3-u_2}\).
  • Set \(h_3=\frac{h_2-h_1}{u_3}\).
  • Set \(u_m = \frac{1}{2}\left(u_2-\frac{h_1}{h_3}\right)\).
  • Set \(g_m = g(\mathbf{x}^{(0)}-u_m\cdot gradx)\).
  • If \(g_3 < g_m\) (calculamos el valor mínimo entre \(g_3\) y \(g_m\))
    • Set \(u_m=u_3\).
    • Set \(g_m=g_3\).

• • Set \(\mathbf{x}^{(1)}=\mathbf{x}^{(0)}-u_m\cdot gradx\). (hallamos el siguiente término de la sucesión)
  • If \(\|\mathbf{x}^{(1)}-\mathbf{x}^{(0)}\| < TOL\) (miramos si hemos llegado a la tolerancia permitida)
    • Print solución aproximada:
    • Print \(\mathbf{x}^{(1)}\).
    • Set \(k=k+1\).
    • STOP.
  • Set \(\mathbf{x}^{(0)}=\mathbf{x}^{(1)}\).
  • Set \(k=k+1\).
• Print: hemos alcanzado el número máximo de iteraciones. El método no converge.
• STOP.