Un método cuasi-Newton para funciones reales de dos variables basado en la minimización del número de condición de la matriz de actualización .

==Un método cuasi-Newton para funciones reales de dos variables basado en la minimización del número de condición de la matriz de actualización .==

'''José Jacobo Oliveros Oliveros, Julio Andrés Acevedo Vázquez, Juan Alberto Escamilla Reyna '''

==Resumen==

In this work, a quasi-Newton method is proposed to solve unconstrained non-linear equations based on the minimization of the condition number of the updating matrix considering the Frobenius norm. The convergence of the method is proved using fixed point theory. Some numerical examples are presented to show the performance of the method, and it is compared with the classical DFP and BFGS methods. The results show that the method proposed here is feasible, and for certain kinds of problems, the solution is obtained using fewer iterations and less computing time than the other ones.  The method is applied to solve systems of ordinary differential equations. The results obtained are compared with the ones obtained with the BFGS method. These results show that the two methods have a similar performance.

'''Keywords:''' cuasi-Newton, condition number.

En este trabajo, se propone un método cuasi-Newton para resolver ecuaciones no lineales sin restricciones, que se basa en la minimización del número de condición de la matriz de actualización, considerando la norma de Frobenius. La convergencia del método es probada usando teoría del punto fijo. Se presentan algunos ejemplos numéricos para mostrar el desempeño del método y es comparado con los métodos clásicos DFP y BFGS. Los resultados muestran que el método propuesto es factible y que para cierta clase de problemas, obtiene la solución utilizando menos iteraciones que los otros. El método es aplicado para resolver sistemas de ecuaciones diferenciales ordinarias. Los resultados obtenidos son comparados con aquellos obtenidos con el método BFGS. Estos resultados muestran que los dos métodos tiene un desempeño similar.

'''Palabras clave:''' cuasi-Newton, número de condición.

==1 Introducción==

Los problemas de minimización de funciones tienen gran importancia por sus múltiples aplicaciones y llevan a la resolución de  sistemas de ecuaciones  ya que tenemos que encontrar los puntos críticos de la función que se minimiza. Resolver ecuaciones tiene gran importancia por si mismo, ya que aparecen con frecuencia  en muchas aplicaciones. En particular, los sistemas de ecuaciones aparecen en los métodos numéricos para resolver ecuaciones diferenciales ordinarias, como el método de Euler implícito que es un caso de los llamados métodos BFD (por sus siglas en inglés Backward Differentiation Formulas) <span id='citeF-1'></span>[[#cite-1|[1]]]. Por lo anterior, se ha puesto especial atención en el desarrollo de métodos para la resolución de estas ecuaciones tanto lineales como no lineales. En este trabajo, se propone un método cuasi-Newton para resolver ecuaciones no lineales y se muestra que para cierta clase de funciones produce mejores resultados que los conocidos métodos DFP y BFGS. También se implementa este método y el  BFGS en el Método Implícito de Euler y se halla que producen resultados similares.

A continuación, se presenta el planteamiento del problema de minimización que interesa en este trabajo.

Sea <math display="inline">\boldsymbol {f}:\mathbb{R}^n\to \mathbb{R}</math> una función suave. Consideremos el siguiente problema

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\min  \boldsymbol {f}(\boldsymbol {x});</math>
|-
| style="text-align: center;" | <math>     \hbox{ s. a}  \boldsymbol {x}\in \mathbb{R}^n.      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

El primer método cuasi-Newton (denotado por DFP) fue creado por Davidon, Fletcher y Powell <span id='citeF-2'></span><span id='citeF-3'></span>[[#cite-2|[2,3]]], demostrando que este método era más eficaz que los métodos existentes. El método DFP es un ''método de búsqueda en la linea'' que considera una aproximación cuadrática de la función objetivo en la iteración <math display="inline">\boldsymbol {x}^k</math>. Tenemos la siguiente notación: <math display="inline">\boldsymbol {f}^k=\boldsymbol {f}(\boldsymbol {x}^k)</math>, <math display="inline">\nabla \boldsymbol {f}^k=\nabla \boldsymbol {f}(\boldsymbol {x}^k)</math> y <math display="inline">\nabla ^2\boldsymbol {f}^k=\nabla ^2\boldsymbol {f}(\boldsymbol {x}^k)</math> donde <math display="inline">\nabla ^2\boldsymbol {f}</math> denota el hessiano de la función <math display="inline">\boldsymbol {f}</math>, por lo que la aproximación está dada por

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {m}_k(\boldsymbol {p})=\boldsymbol {f}^k+(\nabla \boldsymbol {f}^k)^t \boldsymbol {p}+\frac{1}{2}\boldsymbol {p}^t\boldsymbol {B}_k\boldsymbol {p},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

donde <math display="inline">\boldsymbol {B_k}</math> es una matriz simétrica definida positiva de tamaño <math display="inline">n\times n</math>, que se actualizará en cada iteración. La matriz <math display="inline">\boldsymbol {B_K}</math> no es el hessiano de la función <math display="inline">\boldsymbol {f}</math>, pero es una aproximación de ella. Notemos que <math display="inline">\boldsymbol {m}_k(\boldsymbol {0})=\boldsymbol {f}^k</math> y dado que <math display="inline">\boldsymbol {m}_k</math> es una función convexa, alcanza el mínimo en

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {p}^k=-\boldsymbol {B}_k^{-1}\nabla \boldsymbol {f}^k, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

el cual será usado en la nueva iteración

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {x}^{k+1}=\boldsymbol {x}^k+\alpha _k\boldsymbol {p}^k, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

donde la longitud de paso <math display="inline">\alpha _k</math> se escoge de tal forma que cumpla las condiciones de Wolfe.

Un requerimiento que parece razonable pedir en <math display="inline">\boldsymbol {B}_{k+1}</math> es que el gradiente de <math display="inline">\boldsymbol {m}_{k+1}</math> debe coincidir con el gradiente de la función objetivo <math display="inline">\boldsymbol {f}</math> en al menos dos iteraciones <math display="inline">\boldsymbol {x}^k</math> y <math display="inline">\boldsymbol {x}^{k+1}</math>. Obsérvese que <math display="inline">\nabla \boldsymbol {m}_{k+1}(\boldsymbol {0})=\nabla \boldsymbol {f}^{k+1}</math>, así que la segunda de estas condiciones se satisface. La primera condición puede ser escrita como

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {B}_{k+1}(\boldsymbol {x}^{k+1}-\boldsymbol {x}^k) = \nabla \boldsymbol {f}^{k+1}-\nabla \boldsymbol {f}^k. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

Para simplificar la ecuación anterior se definen los vectores <span id="eq-6"></span>

<span id="eq-6.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {s}^k = \boldsymbol {x}^{k+1}-\boldsymbol {x}^k = \alpha _k\boldsymbol {p}^k; </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.a)
|}

<span id="eq-6.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}^k = \nabla \boldsymbol {f}^{k+1}-\nabla \boldsymbol {f}^k. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.b)
|}

Así, la ecuación [[#eq-5|(5)]] se reescribe como

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {B}_{k+1}\boldsymbol {s}^k=\boldsymbol {y}^k, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

la cual es conocida como la ''ecuación secante''.

Del hecho de que la matriz <math display="inline">\boldsymbol {B}_{k+1}</math> es definida positiva y por la ecuación secante, se tiene que

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(\boldsymbol {s}^k)^t\boldsymbol {y}^k>0. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

La desigualdad anterior es conocida como la ''condición de curvatura''.

Cuando la función es fuertemente convexa, la condición de curvatura se cumplirá, para cualesquiera dos puntos <math display="inline">\boldsymbol {x}^{k+1}</math> y <math display="inline">\boldsymbol {x}^k</math>, pero no siempre se cumplirá cuando la función no es convexa. En este último caso, se tiene que forzar a que se cumpla [[#eq-8|(8)]] imponiendo las condiciones de Wolfe o las condiciones fuertes de Wolfe en <math display="inline">\alpha _k</math> <span id='citeF-4'></span>[[#cite-4|[4]]].

Cuando la condición de curvatura se cumple, la ecuación secante [[#eq-7|(7)]] tiene ''infinitas soluciones'', ya que hay <math display="inline">\frac{n(n+1)}{2}</math> grados de libertad en una matriz simétrica y la ecuación secante representa <math display="inline">n</math> condiciones y hay <math display="inline">n</math> condiciones más como consecuencia de que la matriz es definida positiva, pero estas condiciones no absorben los grados de libertad restantes. Así, ''hay infinitas maneras de elegir la matriz de actualización <math>\boldsymbol {B}_{k+1}</math>'', lo que da la posibilidad de proponer diferentes métodos que se adapten al problema que se está estudiando.

==2 Descripción de los métodos DFP y BFGS==

En esta sección se explicará brevemente la manera en que se construyeron los métodos DFP y BFGS usuales. Notamos que, aunque los dos métodos se construyen de manera muy similar, éstos poseen propiedades muy distintas.

===2.1 El método DFP===

Una manera de determinar <math display="inline">\boldsymbol {B}_{k+1}</math> es pidiendo que, entre todas las matrices simétricas que cumplen la ecuación secante, <math display="inline">\boldsymbol {B}_{k+1}</math> debe ser la más cercana a la matriz actual <math display="inline">\boldsymbol {B}_k</math>, en otras palabras, tenemos que resolver el siguiente subproblema <span id="eq-9"></span>

<span id="eq-9.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\min _{\boldsymbol {B}} \left\Vert \boldsymbol {B}-\boldsymbol {B}_k\right\Vert ; </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9.a)
|}

<span id="eq-9.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hbox{s. a } \boldsymbol {B} =\boldsymbol {B}^t;</math>
|-
| style="text-align: center;" | <math>         \boldsymbol {Bs}^k =\boldsymbol {y}^k,     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9.b)
|}

donde <math display="inline">\boldsymbol {s}^k</math> y <math display="inline">\boldsymbol {y}^k</math> satisfacen la condición de curvatura y <math display="inline">\boldsymbol {B}_k</math> es una matriz simétrica definida positiva. Muchas normas pueden ser usadas en [[#eq-9.a|(9.a)]] y cada una de ellas da lugar a un método cuasi-Newton distinto. Una norma que permite una fácil resolución del problema es la ''norma pesada de Frobenius''

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\Vert \boldsymbol {A}\right\Vert _{\boldsymbol {W}}^2=tr\left(\boldsymbol {WA}^t\boldsymbol {WA}\right). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

El peso <math display="inline">\boldsymbol {W}</math> puede ser escogido como cualquier matriz que satisfaga la relación <math display="inline">\boldsymbol {Wy}^k=\boldsymbol {s}^k</math>. Se asume <math display="inline">\boldsymbol {W}=\overline{\boldsymbol {G}}_k^{-1}</math>, donde <math display="inline">\overline{\boldsymbol {G}}_k</math> es el hessiano promedio, dado por

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline{\boldsymbol {G}}_k=\int _0^1 \nabla ^2 \boldsymbol {f}(\boldsymbol {x}^k+\tau \alpha _k \boldsymbol {p}^k)d\tau{.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

Fácilmente puede verificarse que <math display="inline">\overline{\boldsymbol {G}}_k^{-1}\boldsymbol {y}^k=\boldsymbol {s}^k</math>.

Con la norma y matriz de peso descritas arriba, la única solución de [[#eq-9|(9)]] es

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {B}_{k+1}^{DFP}=(\boldsymbol {I}-\gamma _k\boldsymbol {y}^k(\boldsymbol {s}^k)^t)\boldsymbol {B}_k(\boldsymbol {I}-\gamma _k\boldsymbol {s}^k(\boldsymbol {y}^k)^t)+\gamma _k\boldsymbol {y}^k(\boldsymbol {y}^k)^t, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

donde <math display="inline">\gamma _k=\dfrac{1}{(\boldsymbol {y}^k)^t\boldsymbol {s}^k}.</math>

Esta fórmula es llamada la ''fórmula de actualización DFP'', ya que fue propuesta por Davidon en 1959 <span id='citeF-2'></span>[[#cite-2|[2]]], y posteriormente estudiada, implementada y popularizada por Fletcher y Powell <span id='citeF-3'></span>[[#cite-3|[3]]].

Sea <math display="inline">\boldsymbol {H}_k=\boldsymbol {B}_k^{-1}</math>. Usando la fórmula de Sherman-Morrison-Woodbury <span id='citeF-4'></span>[[#cite-4|[4]]] y la ecuación [[#eq-12|(12)]] se sigue que

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {H}_{k+1}^{DFP}=\boldsymbol {H}_k-\frac{\boldsymbol {H}_k\boldsymbol {y}^k(\boldsymbol {y}^k)^t\boldsymbol {H}_k}{(\boldsymbol {y}^k)^t\boldsymbol {H}_k\boldsymbol {y}^k}+ \frac{\boldsymbol {s}^k(\boldsymbol {s}^k)^t}{(\boldsymbol {y}^k)^t\boldsymbol {s}^k}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

Esto resulta útil, ya que nos permite calcular la dirección de búsqueda simplemente con una multiplicación matriz-vector.

Nótese que los últimos dos términos del lado derecho de [[#eq-13|(13)]] son matrices de rango 1, así <math display="inline">\boldsymbol {H}_k</math> es una modificación de rango 2. Esa es la idea fundamental de la actualización cuasi-Newton, en lugar de recalcular la actualización desde cero, se aplica una simple modificación que combina la información recientemente observada de la función objetivo con el conocimiento existente en nuestra aproximación de la hessiana actual.


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:Punto inicial <math display="inline">\boldsymbol {x}^0</math>, tolerancia <math display="inline">\varepsilon </math>, aproximación del inverso del hessiano <math display="inline">\boldsymbol {H}_0</math>. <math display="inline">\boldsymbol {x}^*</math> mínimo de <math display="inline">\boldsymbol {f}</math>. <math display="inline">k\leftarrow 0.</math> <math>\left\Vert \nabla \boldsymbol {f}(\boldsymbol {x}^k)\right\Vert  > \varepsilon </math>Calcular la dirección de búsqueda

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {p}^k=-\boldsymbol {H}_k\nabla \boldsymbol {f}(\boldsymbol {x}^k). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Calcular <math display="inline">\boldsymbol {x}^{k+1}=\boldsymbol {x}^k+\alpha _k\boldsymbol {p}^k</math>, donde <math display="inline">\alpha _k</math> es calculado mediante un procedimiento de búsqueda en la línea que satisfaga las condiciones de Wolfe, <span id="eq-15"></span>


<span id="eq-15.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {f}(\boldsymbol {x}_k+\alpha \boldsymbol {p}_k) \leq \boldsymbol {f}(\boldsymbol {x}_k)+c_1\alpha \nabla \boldsymbol {f}_k^t\boldsymbol {p}_k, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15.a)
|}



<span id="eq-15.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla \boldsymbol {f}(\boldsymbol {x}_k+\alpha _k \boldsymbol {p}_k)^T\boldsymbol {p}_k \geq  c_2 \nabla \boldsymbol {f}_k^t\boldsymbol {p}_k, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15.b)
|}

con <math display="inline">0<c_1<c_2<1.</math> Definir <math display="inline">\boldsymbol {s}^k=\boldsymbol {x}^{k+1}-\boldsymbol {x}^k</math> y <math display="inline">\boldsymbol {y}^k=\nabla \boldsymbol {f}(\boldsymbol {x}^{k+1})-\nabla \boldsymbol {f}(\boldsymbol {x}^k)</math>. Calcular <math display="inline">\boldsymbol {H}_{k+1}</math> por medio de [[#eq-13|(13)]].  <math display="inline">k\leftarrow k+1.</math> 


|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-1'></span>'''Algoritmo. 1''' Algoritmo DFP
|}

===2.2 El método BFGS===

La fórmula del método BFGS puede ser derivado haciendo un simple cambio en el argumento que lleva a [[#eq-12|(12)]]. Las condiciones impuestas en las aproximaciones del hessiano <math display="inline">\boldsymbol {B}_k</math>, pueden ser impuestas en sus inversas <math display="inline">\boldsymbol {H}_k</math>. Por lo tanto, la actualización <math display="inline">\boldsymbol {H}_{k+1}</math> debe ser simétrica, definida positiva, y satisfacer la ecuación secante [[#eq-7|(7)]], ahora escrita como <math display="inline">\boldsymbol {H}_{k+1}\boldsymbol {y}^k=\boldsymbol {s}^k.</math>

La condición de la cercanía a <math display="inline">\boldsymbol {H}_k</math> es ahora especificada de forma análoga a [[#eq-9|(9)]]

<span id="eq-16"></span>

<span id="eq-16.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\min _{\boldsymbol {H}} \left\Vert \boldsymbol {H}-\boldsymbol {H}_k\right\Vert ; </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16.a)
|}

<span id="eq-16.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hbox{s. a } \boldsymbol {H} =\boldsymbol {H}^t;</math>
|-
| style="text-align: center;" | <math>         \boldsymbol {Hy}^k =\boldsymbol {s}^k.     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16.b)
|}

La norma considerada es de nuevo la norma pesada de Frobenius descrita anteriormente, donde la matriz de peso <math display="inline">\boldsymbol {W}</math> es ahora cualquier matriz que satisfaga <math display="inline">\boldsymbol {Ws}^k=\boldsymbol {y}^k</math>. (Por concreción, se asume de nuevo que <math display="inline">\boldsymbol {W}</math> está dado por el inverso del hessiano promedio <math display="inline">\overline{\boldsymbol {G}}_k</math> definida en [[#eq-11|(11)]]). La única solución <math display="inline">\boldsymbol {H}_{k+1}</math> de [[#eq-16|(16)]] está dada por

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {H}_{k+1}^{BFGS}=\left(\boldsymbol {I}-\rho _k\boldsymbol {s}^k(\boldsymbol {y}^k)^t\right)\boldsymbol {H}_k\left(\boldsymbol {I}-\rho _k\boldsymbol {y}^k(\boldsymbol {s}^k)^t\right)+\rho _k\boldsymbol {s}^k(\boldsymbol {s}^k)^t, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

donde

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho _k=\frac{1}{(\boldsymbol {y}^k)^t\boldsymbol {s}^k}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

o bien, desarrollando las multiplicaciones

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {H}_{k+1}^{BFGS}=\boldsymbol {H}_k+\left(1+\frac{(\boldsymbol {y}^k)^t\boldsymbol {H}_k\boldsymbol {y}^k}{(\boldsymbol {s}^k)^t\boldsymbol {y}^k}\right)\frac{\boldsymbol {s}^k(\boldsymbol {s}^k)^t}{(\boldsymbol {s}^k)^t\boldsymbol {y}^k}-\left(\frac{\boldsymbol {s}^k(\boldsymbol {y}^k)^t\boldsymbol {H}_k+\boldsymbol {H}_k\boldsymbol {y}^k(\boldsymbol {s}^k)^t}{(\boldsymbol {s}^k)^t\boldsymbol {y}^k}\right). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

La demostración del siguiente teorema puede encontrarse en <span id='citeF-5'></span>[[#cite-5|[5]]].

<span id='theorem-3.3.2'></span>Teorema 1:  Si la fórmula de actualización <math display="inline">BFGS</math> [[#eq-19|(19)]], es ahora escrita como <math display="inline">\boldsymbol {H}_{k+1}=\boldsymbol {H}_k+\boldsymbol {E}</math>, entonces <math display="inline">\boldsymbol {E}</math> resuelve el problema <span id="eq-20"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\min _{\boldsymbol {E}} \left\Vert \boldsymbol {E}\right\Vert _{\boldsymbol {W}};</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (20.a)
|-
| style="text-align: center;" | <math> s. a   \boldsymbol {E} = \boldsymbol {E}^t,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (20.b)
|-
| style="text-align: center;" | <math>             \boldsymbol {Ey}^k = \boldsymbol {\eta }, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (20.c)
|}
|}

donde <math display="inline">\boldsymbol {W}</math> satisface <math display="inline">\boldsymbol {Ws}^k=\boldsymbol {y}^k</math> &nbsp;y &nbsp;<math display="inline">\boldsymbol {\eta }=\boldsymbol {s}^k-\boldsymbol {H}_k\boldsymbol {y}^k</math>.

De acuerdo con Nocedal <span id='citeF-4'></span>[[#cite-4|[4]]], si <math display="inline">\boldsymbol {H}_k</math> es definida positiva, entonces <math display="inline">\boldsymbol {H}_{k+1}</math> es definida positiva.


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:Punto inicial <math display="inline">\boldsymbol {x}_0</math>, tolerancia <math display="inline">\varepsilon </math>, aproximación del inverso del hessiano <math display="inline">\boldsymbol {H}_0</math>. <math display="inline">\boldsymbol {x}^*</math> mínimo de <math display="inline">\boldsymbol {f}</math>. <math display="inline">k\leftarrow 0.</math> <math>\left\Vert \nabla \boldsymbol {f}(\boldsymbol {x}^k)\right\Vert  > \varepsilon </math>Calcular la dirección de búsqueda

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {p}^k=-\boldsymbol {H}_k\nabla \boldsymbol {f}(\boldsymbol {x}^k). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

Calcular <math display="inline">\boldsymbol {x}^{k+1}=\boldsymbol {x}^k+\alpha _k\boldsymbol {p}^k</math>, donde <math display="inline">\alpha _k</math> es calculado mediante un procedimiento de búsqueda en la línea que satisfaga las condiciones de Wolfe [[#algorithm-1|(1)]]. Definir <math display="inline">\boldsymbol {s}^k=\boldsymbol {x}^{k+1}-\boldsymbol {x}^k</math> y <math display="inline">\boldsymbol {y}^k=\nabla \boldsymbol {f}(\boldsymbol {x}^{k+1})-\nabla \boldsymbol {f}(\boldsymbol {x}^k)</math>. Calcular <math display="inline">\boldsymbol {H}_{k+1}</math> por medio de [[#eq-19|(19)]].  <math display="inline">k\leftarrow k+1.</math> 


|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-2'></span>'''Algoritmo. 2''' Algoritmo BFGS
|}

<!-- comment

===Algoritmo BFGS===

Sólo un problema tiene que ser resuelto antes de que podamos definir un algoritmo BFGS completo ¿Cómo debemos escoger la aproximación inicial <math display="inline">H_0</math>? Desafortunadamente, no hay una fórmula que funcione bien en todos los casos. Podemos usar información especifica acerca del problema, por ejemplo, configurándolo como el inverso de una aproximación del hessiano calculado por diferencias finitas en <math display="inline">x_0</math>. De otra manera, podemos simplemente configurarla para que sea la matriz identidad, o un múltiplo de la matriz identidad, donde el múltiplo es escogido para reflejar el escalamiento de las variables.


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:Punto inicial <math display="inline">x_0</math>, tolerancia <math display="inline">\varepsilon </math>, aproximación del inverso de el hessiano <math display="inline">H_0</math>. <math display="inline">x_*</math> <math display="inline">k\leftarrow 0;</math> <math>\Vert \nabla f_k\Vert > \varepsilon </math> Calcular la dirección de búsqueda

<span id=""></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_k=-H_k\nabla f_k; </math>
|}
|}

Configurar <math display="inline">x_{k+1}=x_k+\alpha _kp_k</math>, donde <math display="inline">\alpha _k</math> es calculado mediante un procedimiento de búsqueda en la línea que satisfaga las condiciones de Wolf; Definir <math display="inline">s^k=x_{k+1}-x_k</math> y <math display="inline">y_k=\nabla f_{k+1}-\nabla f_k</math>; Calcule <math display="inline">H_{k+1}</math> por medio de [[#eq-17|(17)]]; <math display="inline">k\leftarrow k+1;</math> 


|-
| style="text-align: center; font-size: 75%;"|
'''Algoritmo. ''' Algoritmo BFGS
|}

Cada iteración puede ser ejecutado con un costo de <math display="inline">O(n^2)</math> operaciones aritméticas (más el costo de las evaluaciones en la función y el gradiente); no hay operaciones de <math display="inline">O(n^3)</math>, tales como solución de un sistema de ecuaciones lineal u operaciones matriz-matriz. El algoritmo es robusto, y su rapidez de convergencia es superlineal, lo cual es suficientemente rápido para la mayoría de los propósitos prácticos. Incluso aunque el método de Newton converge más rápido a la solución (es decir, cuadráticamente), su costo por iteración es mayor, ya que requiere la solución de un sistema lineal. Una ventaja muy importante del método BFGS es, por supuesto, que no requiere el cálculo de segundas derivadas. 

-->

==3 Propuesta del nuevo método==

Como se vio en la sección anterior, se puede conseguir distintos métodos cuasi-Newton si se cambia el subproblema que determina la matriz de actualización, tal subproblema puede adaptarse a la naturaleza misma de la función. Cuando se considera funciones con error, los métodos anteriores pueden presentar problemas de convergencia o de velocidad, por ello, en esta sección se cambia el subproblema y se le pide a la matriz de actualización que tenga el menor número de condición. En <span id='citeF-6'></span>[[#cite-6|[6]]], se utiliza también un criterio que se basa en el menor número de condición para el método de gradiente conjugado.

===3.1 Motivación===

Considere la función <math display="inline">\boldsymbol {\bar{f}}:\mathbb{R}^2\to \mathbb{R}</math> definida por

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\bar{f}}(x,y)=\boldsymbol {f}(x,y)+\boldsymbol {e}(x,y),     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

donde <math display="inline">\boldsymbol {f}(x,y)=0.00001(x^4+y^4)</math>&nbsp;y&nbsp;<math display="inline">\boldsymbol {e}(x,y)= 0.1\boldsymbol {f}(x,y)sin(  x^2 + y^2 )</math>, el punto inicial <math display="inline">\boldsymbol {x}^0=(3,1)^t</math>, una tolerancia  de <math display="inline">10^{-8}</math> y la matriz identidad como matriz inicial. En este ejemplo, <math display="inline">\boldsymbol {e}(x,y)</math> representa un error que emula los diferentes errores que se pueden presentar, tales como error de mediciones, truncamientos, redondeo de la máquina, etc. Los resultados de la implementación se muestran a continuación

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-PP.png|390px|Resultados de los métodos DFP y BFGS para el punto inicial x⁰=(3,1)<sup>t</sup>.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 1:''' Resultados de los métodos DFP y BFGS para el punto inicial <math>\boldsymbol {x}^0=(3,1)^t</math>.
|}

Se observa que los métodos DFP y BFGS no convergen a la solución del problema, lo cual es debido a la presencia de errores. Esto nos lleva a plantear la posibilidad de elegir a la matriz <math display="inline">\boldsymbol {B}_{k+1}</math> de una manera diferente.

===3.2 Planteamiento del problema===

Como se vio anteriormente, el método DFP puede ser modificado si se cambia el problema que lleva a la obtención de la actualización <math display="inline">\boldsymbol {B} _{k+1}</math>. Tomando en cuenta que la presencia de errores provocó que los métodos DFP y BFGS no convergieran, un requerimiento que parece razonable, es pedir que la siguiente matriz tenga el menor número de condición, es decir, que en lugar de resolver [[#eq-9|(9)]] se resuelve <span id="eq-23"></span>

<span id="eq-23.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\min _{\boldsymbol {B} } cond(\boldsymbol {B} ); </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23.a)
|}

<span id="eq-23.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hbox{s. a }\boldsymbol {B} =\boldsymbol {B} ^t;</math>
|-
| style="text-align: center;" | <math>          \boldsymbol {{B}s}^k=\boldsymbol {y}^k,     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23.b)
|}

donde <math display="inline">cond(\boldsymbol {B} )</math> está dado por <math display="inline">\frac{\lambda _n}{\lambda _1}</math>, con <math display="inline">0<\lambda _1<\cdots{<\lambda}_n</math> valores propios de la matriz <math display="inline">\boldsymbol {B} </math> (se verá que la matriz obtenida es definida positiva y por lo tanto la dirección de paso será de descenso). Se enfatiza cuál es la razón de usar este método:

''Al pedir el menor número de condición, se trata de que la búsqueda de la dirección de descenso por esta vía, sea menos sensible a los errores.''

===3.3 Caso bidimensional===

Supongamos que <math display="inline">\boldsymbol {f}</math> es una función de <math display="inline">\mathbb{R}^2</math> a <math display="inline">\mathbb{R}</math>. Consideramos <math display="inline">\boldsymbol {y}^k=(y_1^k,y_2^k)</math>, <math display="inline">\boldsymbol {s}^k=(s_1^k,s_2^k)</math> y

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {B} = \begin{bmatrix}b_{11} & b_{12}\\ b_{12} & b_{22}  \end{bmatrix}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

En este caso el número de condición está dado por

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>cond(\boldsymbol {B} )=\frac{\lambda _2}{\lambda _1}=\frac{tr(\boldsymbol {B} )+\sqrt{tr(\boldsymbol {B} )^2-4det(\boldsymbol {B} )}}{tr(\boldsymbol {B} )-\sqrt{tr(\boldsymbol {B} )^2-4det(\boldsymbol {B} )}}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

donde <math display="inline">tr(\boldsymbol {B} )</math> y <math display="inline">det(\boldsymbol {B} )</math> denotan la traza y el determinante de <math display="inline">\boldsymbol {B} </math>, respectivamente. Ahora bien, dado que la matriz tiene que ser simétrica y cumplir la ecuación secante, tenemos que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b_{11}s_1^k+b_{12}s_2^k=y_1^k, </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b_{12}s_1^k+b_{22}s_2^k=y_2^k, </math>
|}
|}

de donde se obtiene que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b_{11}=\frac{y_1^k-b_{12}s_2^k}{s_1^k}, </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b_{22}=\frac{y_2^k-b_{12}s_1^k}{s_2^k}, </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>tr(\boldsymbol {B} )=\frac{y_1^k}{s_1^k}+\frac{y_2^k}{s_2^k}-b_{12}\left(\frac{s_2^k}{s_1^k}+\frac{s_1^k}{s_2^k}\right), </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>det(\boldsymbol {B} )=\left(\frac{y_1^k-b_{12}s_2^k}{s_1^k}\right)\left(\frac{y_2^k-b_{12}s_1^k}{s_2^k}\right)-b_{12}^2. </math>
|}
|}

Sustituyendo las ecuaciones anteriores en [[#eq-25|(25)]], se obtiene una función de <math display="inline">\mathbb{R}</math> a <math display="inline">\mathbb{R}</math>, la cual se minimiza calculando la primera derivada (respecto a <math display="inline">b_{12}</math>) e igualándola a 0. Con lo cual se obtiene que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b_{12}=\frac{y_1^ky_2^k\left\Vert \boldsymbol {s}^k\right\Vert ^2-s_1^ks_2^k\left\Vert \boldsymbol {y}^k\right\Vert ^2}{\left\Vert \boldsymbol {s}^k\right\Vert ^2(\boldsymbol {s}^k)^t\boldsymbol {y}^k}, </math>
|}
|}

que al sustituirlo en las expresiones de <math display="inline">b_{11}</math> y <math display="inline">b_{22}</math> se llega a

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b_{11}=\frac{(y_1^k)^2\left\Vert \boldsymbol {s}^k\right\Vert ^2 + (s_2^k)^2\left\Vert \boldsymbol {y}^k\right\Vert ^2}{\left\Vert \boldsymbol {s}^k\right\Vert ^2(\boldsymbol {s}^k)^t\boldsymbol {y}^k},     b_{22}=\frac{(y_2^k)^2\left\Vert \boldsymbol {s}^k\right\Vert ^2 + (s_1^k)^2\left\Vert \boldsymbol {y}^k\right\Vert ^2}{\left\Vert \boldsymbol {s}^k\right\Vert ^2(\boldsymbol {s}^k)^t\boldsymbol {y}^k}. </math>
|}
|}

Con lo cual se obtienen las entradas de la matriz a la cual llamaremos <math display="inline">\boldsymbol {B}_{k+1}^{MCN}</math>, ya que es la que tiene el ''menor número de condición (minor condition number)'', la cual resuelve el subproblema [[#eq-23|(23)]]. Además, sustituyendo el valor de <math display="inline">b_{12}</math> en la expresión del determinante, se obtiene que <math display="inline">\boldsymbol {B} _{k+1}^{MCN}</math> es definida positiva. Haciendo los cálculos y simplificaciones necesarias se sigue que

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {B} _{k+1}^{MCN}=\frac{\boldsymbol {y}^k(\boldsymbol {y}^k)^t}{(\boldsymbol {s}^k)^t\boldsymbol {y}^k}+\frac{\left\Vert \boldsymbol {y}^k\right\Vert ^2\boldsymbol {z}^k (\boldsymbol {z}^k)^t }{\left\Vert \boldsymbol {s}^k\right\Vert ^2(\boldsymbol {s}^k)^t\boldsymbol {y}^k}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

donde <math display="inline">(\boldsymbol {z}^k)^t=((\boldsymbol {s}^k)^\perp )^t=(s_2^k,-s_1^k)</math>. Notemos que, tanto el primer sumando como el segundo son matrices de rango 1.

Con el fin de que en la dirección de búsqueda no se tenga que calcular el inverso de la matriz anterior, denotamos por <math display="inline">\boldsymbol {H}_k=\boldsymbol {B}_k^{-1}</math>, teniendo así una expresión del inverso de [[#eq-29|(29)]]

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {H}_{k+1}^{MCN}=\frac{\boldsymbol {s}^k(\boldsymbol {s}^k)^t}{(\boldsymbol {s}^k)^t\boldsymbol {y}^k}+\frac{\left\Vert \boldsymbol {s}^k\right\Vert ^2\boldsymbol {w}^k (\boldsymbol {w}^k)^t }{\left\Vert \boldsymbol {y}^k\right\Vert ^2(\boldsymbol {s}^k)^t\boldsymbol {y}^k}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

donde <math display="inline">(\boldsymbol {w}^k)^t=((\boldsymbol {y}^k)^\perp )^t=(y_2^k,-y_1^k)</math>.

====3.3.1 Existencia de la matriz de peso W.====

Siguiendo la idea de la construcción del método DFP, lo que se requiere es que <math display="inline">\boldsymbol {B} _{k+1}^{MCN}</math> sea la matriz que resuelva el problema [[#eq-9|(9)]] con la misma norma <math display="inline">\left\Vert \boldsymbol {A}\right\Vert _{\boldsymbol {W}}=tr(\boldsymbol {WA}^t\boldsymbol {WA})</math>, pero con una matriz de peso adecuada. Así, el problema ahora es encontrar <math display="inline">\boldsymbol {W}</math> de tal forma que <math display="inline">\boldsymbol {Wy}^k=\boldsymbol {s}^k</math> y

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {B} _{k+1}^{MCN} =\arg \min _{\boldsymbol {B} } \left\Vert \boldsymbol {B} -\boldsymbol {B} _k\right\Vert _{\boldsymbol {W}};</math>
|-
| style="text-align: center;" | <math> \hbox{ s. a }\boldsymbol {B} =\boldsymbol {B} ^t;</math>
|-
| style="text-align: center;" | <math>          \boldsymbol {{B}s}^k=\boldsymbol {y}^k. </math>
|}
|}

Notando que <math display="inline">\boldsymbol {B} </math> y <math display="inline">\boldsymbol {B} _k</math> son matrices simétricas se tiene que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\Vert \boldsymbol {B} -\boldsymbol {B} _k\right\Vert _{\boldsymbol {W}}=tr\left(\left(\boldsymbol {W}(\boldsymbol {B} -\boldsymbol {B} _k)\right)^2 \right).</math>
|}
|}

Además, de la restricción <math display="inline">\boldsymbol {Wy}^k=\boldsymbol {s}^k</math> y haciendo un análisis análogo al realizado para <math display="inline">\boldsymbol {B} _{k+1}^{MCN}</math> se obtiene que

<span id="eq-31"></span>
<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w_{11}=\frac{s_1^k-w_{12}y_2^k}{y_1^k},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|-
| style="text-align: center;" | <math>     w_{22}=\frac{s_2^k-w_{12}y_1^k}{y_2^k}.  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}
|}

Si se denota por <math display="inline">c_1=b_{11}-b_{11}^k</math>, <math display="inline">c_2=b_{12}-b_{12}^k</math> y <math display="inline">c_3=b_{22}-b_{22}^k</math>,  se tiene que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {g}(w_{12}) = \left\Vert \boldsymbol {B} -\boldsymbol {B} _k\right\Vert _{\boldsymbol {W}}</math>
|-
| style="text-align: center;" | <math>     = 2\left(c_1w_{12} +\frac{c_2(s_2^k - w_{12}y_1^k)}{y_2^k}\right)\left(c_3w_{12} +\frac{c_2(s_1^k - w_{12}y_2^k)}{y_1^k}\right)</math>
|-
| style="text-align: center;" | <math>     + \left(c_2w_{12} + \frac{c_1(s_1^k - w_{12}y_2^k)}{y_1^k}\right)^2 + \left(c_2w_{12} +\frac{c_3(s_2^k - w_{12}y_1^k)}{y_2^k}\right)^2. </math>
|}
|}

Obteniendo los puntos estacionarios de la expresión anterior (con respecto a <math display="inline">b_{ij}</math>), obtenemos que

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w_{12}=\frac{b^k_{22}(s_2^k)^2(y_1^k)^2 - 2b^k_{12}(s_2^k)^2(y_1^k)(y_2^k) + b^k_{11}(s_2^k)^2(y_2^k)^2 + (s_2^k)(y_1^k)^2(y_2^k) + (s_1^k)(y_1^k)^3}{((\boldsymbol {s}^k)^t\boldsymbol {y}^k)^2}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

Sustituyendo ([[#eq-33|33]]) en ([[#eq-31|31]]) y ([[#eq-32|32]]), se halla a la matriz de peso <math display="inline">\boldsymbol {W}</math>.

==4 Algoritmo==

Definimos ahora, un algoritmo que encuentre la solución de [[#eq-1|(1)]]


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:Punto inicial <math display="inline">\boldsymbol {x}^0</math>, tolerancia <math display="inline">\varepsilon </math>, aproximación del inverso del hessiano <math display="inline">\boldsymbol {H} _0</math>. <math display="inline">\boldsymbol {x}^*</math> mínimo de <math display="inline">\boldsymbol {f}</math>. <math display="inline">k\leftarrow 0.</math> <math>\left\Vert \nabla \boldsymbol {f}(\boldsymbol {x}^k)\right\Vert  > \varepsilon </math>Calcular la dirección de búsqueda

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {p}^k=-\boldsymbol {H}_k\nabla \boldsymbol {f}(\boldsymbol {x}^k). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

Calcular <math display="inline">\boldsymbol {x}^{k+1}=\boldsymbol {x}^k+\alpha _k\boldsymbol {p}^k</math>, donde <math display="inline">\alpha _k</math> es calculado mediante un procedimiento de búsqueda en la línea que satisfaga las condiciones de Wolfe [[#algorithm-1|(1)]]. Definir <math display="inline">\boldsymbol {s}^k=\boldsymbol {x}^{k+1}-\boldsymbol {x}^k</math> y <math display="inline">\boldsymbol {y}^k=\nabla \boldsymbol {f}(\boldsymbol {x}^{k+1})-\nabla \boldsymbol {f}(\boldsymbol {x}^k)</math>. Calcular <math display="inline">\boldsymbol {H} _{k+1}</math> por medio de [[#eq-30|(30)]].  <math display="inline">k\leftarrow k+1.</math> 


|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-3'></span>'''Algoritmo. 3''' Algoritmo MCN
|}

==5 Convergencia==

Para analizar la convergencia del método MCN, usaremos los resultados dados en  <span id='citeF-7'></span>[[#cite-7|[7]]], en donde los métodos cuasi-Newton son analizados utilizando técnicas de punto fijo, se demuestran resultados de convergencia y se hallan tasas de convergencia.

En lo que sigue, se dará un resumen de las definiciones y resultados que están en  <span id='citeF-7'></span>[[#cite-7|[7]]]. Con esto, podremos demostrar la convergencia que nos interesa.  Aunque los resultados mostrados en esta sección, son utilizados para resolver sistemas de ecuaciones no lineales, estos pueden ser utilizados en la optimización de funciones ya que nos interesa resolver el sistema <math display="inline">\nabla \boldsymbol {f}(\boldsymbol {x})=\boldsymbol {0}</math>.

Sean <math display="inline">X</math> un espacio lineal (vectorial) de dimensión finita, <math display="inline">D\subset X</math>, <math display="inline">\Omega \subset \mathbb{R}^n</math>, <math display="inline">\boldsymbol {\phi }:\Omega \times D\to \mathbb{R}^n</math>. Se consideran los métodos iterativos definidos por

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {x}^{k+1}=\boldsymbol {\phi }(\boldsymbol {x}^k,E_k). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

<span id='theorem-punto fijo'></span>Definición 1: Decimos que <math display="inline">\boldsymbol {x}\in \Omega </math> es un punto fijo de <math display="inline">\boldsymbol {\phi }</math> si y sólo si,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {x}=\boldsymbol {\phi }(\boldsymbol {x},E),\quad \hbox{para todo }E\in D.     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

El objetivo de los algoritmos del tipo [[#eq-35|(35)]] es el de aproximar los puntos fijos de <math display="inline">\boldsymbol {\phi }</math>. Notemos que los métodos cuasi-Newton pertenecen a esta clase de métodos, ya que si consideramos <math display="inline">D=\lbrace \boldsymbol {B}\in \mathbb{R}^{n\times n}:\boldsymbol {B}\hbox{ es no singular}\rbrace </math>, <math display="inline">\Omega =\mathbb{R}^n</math>, <math display="inline">\boldsymbol {F}:\Omega \to \mathbb{R}^n</math>, definimos

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\phi }(\boldsymbol {x},\boldsymbol {B})=\boldsymbol {x}-\boldsymbol {B} ^{-1}\boldsymbol {F}(\boldsymbol {x}). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

Claramente, el conjunto de puntos <math display="inline">\boldsymbol {x} \in \mathbb{R}^n</math> tales que <math display="inline">\boldsymbol {x}=\boldsymbol {\phi }(\boldsymbol {x},\boldsymbol {B})</math> para todo <math display="inline">\boldsymbol {B}\in D</math>, es el conjunto de soluciones de <math display="inline">\boldsymbol {F}(\boldsymbol {x})=\boldsymbol {0}</math>. Los métodos definidos por <math display="inline">\boldsymbol {x}^{k+1}=\boldsymbol {\phi }(\boldsymbol {x}^k,\boldsymbol {B}_k)</math> forman la familia de métodos cuasi-Newton para resolver sistemas de ecuaciones no lineales.

Se establecen condiciones suficientes para que la sucesión definida por [[#eq-35|(35)]] esté bien definida y converja a un punto fijo de <math display="inline">\boldsymbol {\phi }</math> con una convergencia lineal. Se consideran dos suposiciones, A1 y A2 que se establecen a continuación.

'''Suposición A1.'''  Sea <math display="inline">\Omega </math> un conjunto abierto y convexo, <math display="inline">\left\vert \cdot \right\vert </math> una norma en <math display="inline">\mathbb{R}^n</math>, <math display="inline">X</math> un espacio lineal de dimensión finita, <math display="inline">\left\Vert \cdot \right\Vert </math> una norma en <math display="inline">X</math>, <math display="inline">D</math> un subconjunto abierto de <math display="inline">X</math>. Supongamos que <math display="inline">\boldsymbol {\phi }:\Omega \times D\to \mathbb{R}^n</math> es continua, y que la derivada de <math display="inline">\boldsymbol {\phi }</math> con respecto a <math display="inline">\boldsymbol {x}</math> existe y es continua, para todo <math display="inline">\boldsymbol {x}\in \Omega </math>, <math display="inline">E\in D</math>. Denotamos por <math display="inline">\boldsymbol {\phi }'(\boldsymbol {x},\boldsymbol {E})\in \mathbb{R}^{n\times n}</math> la matriz Jacobiana de <math display="inline">\boldsymbol {\phi }</math> con respecto a <math display="inline">\boldsymbol {x}</math>. Supongamos que:

<ol>

<li> Existe <math display="inline">\boldsymbol {x}^*\in \Omega </math>, <math display="inline">L,p>0</math> tales que

<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\left\vert \boldsymbol {\phi }'(\boldsymbol {x},E)-\boldsymbol {\phi }'(\boldsymbol {x}^*,E)\right\vert \leq L\left\vert \boldsymbol {x}-\boldsymbol {x}^*\right\vert ^p,         </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}</li>

para todo <math display="inline">\boldsymbol {x}\in \Omega </math>, <math display="inline">E\in D</math>. Esto implica que para todo <math display="inline">\boldsymbol {x},\boldsymbol {z}\in \Omega </math>, <math display="inline">E\in D</math>,

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\left\vert \boldsymbol {\phi }(\boldsymbol {z},E)-\boldsymbol {\phi }(\boldsymbol {x},E)-\boldsymbol {\phi }'(\boldsymbol {x}^*,E)(\boldsymbol {z}-\boldsymbol {x})\right\vert \leq L\left\vert \boldsymbol {z}-\boldsymbol {x}\right\vert \sigma (\boldsymbol {x},\boldsymbol {z})^p,         </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

donde

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\sigma (\boldsymbol {x},\boldsymbol {z})=\max \left\{\left\vert \boldsymbol {x}-\boldsymbol {x}^*\right\vert ,\left\vert \boldsymbol {z}-\boldsymbol {x}^*\right\vert \right\}.         </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}
<li>Para todo <math display="inline">E\in D</math>,

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\boldsymbol {x}^*=\boldsymbol {\phi }(\boldsymbol {x}^*,E).         </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}</li>

</ol>

'''Suposición A2.''' Existe <math display="inline">E_*\in D</math> tal que

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\vert \boldsymbol {\phi }'(\boldsymbol {x}^*,E_*)\right\vert =r^*<1. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

<span id='theorem-3.6'></span>Lema 1: Suponga que  &nbsp;<math display="inline">\boldsymbol {\phi },\boldsymbol {x}^*,E_*</math> satisfacen las Suposiciones A1 y A2. Sea <math display="inline">r\in (r^*,1)</math>. Entonces, existen <math display="inline">\overline{\varepsilon }=\overline{\varepsilon }(r), \overline{\delta }=\overline{\delta }(r)</math> tales que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\vert \boldsymbol {\phi }(\boldsymbol {x},E)-\boldsymbol {x}^*\right\vert \leq \left\vert \boldsymbol {x}-\boldsymbol {x}^*\right\vert , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

siempre que <math display="inline">\left\vert \boldsymbol {x}-\boldsymbol {x}^*\right\vert \leq \overline{\varepsilon }</math>, <math display="inline">\left\Vert E-E_*\right\Vert \leq \overline{\delta }</math>.

<span id='theorem-T3.1'></span>Teorema 2:  Suponga que <math display="inline">\boldsymbol {\phi },\boldsymbol {x}^*,E_*</math> satisfacen las Suposiciones A1 y A2. Sea <math display="inline">r\in (r^*,1)</math>. Entonces, existen <math display="inline">\varepsilon =\varepsilon (r), \delta =\delta (r)</math> tales que, si <math display="inline">\left\vert \boldsymbol {x}^0-\boldsymbol {x}^*\right\vert \leq \varepsilon </math> y <math display="inline">\left\Vert E_k-E_*\right\Vert \leq \delta </math> para todo <math display="inline">k=0,1,\ldots </math>, entonces la sucesión generada por [[#eq-35|(35)]] está bien definida, converge a <math display="inline">\boldsymbol {x}^*</math> y satisface

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\vert \boldsymbol {x}^{k+1}-\boldsymbol {x}^*\right\vert \leq \left\vert \boldsymbol {x}^k-\boldsymbol {x}^*\right\vert , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

para todo <math display="inline">k=0,1,2,\ldots </math>.

Para el caso de los métodos cuasi-Newton, a fin  de verificar la condición [[#eq-42|(42)]], se debe calcular <math display="inline">\boldsymbol {\phi }'(\boldsymbol {x},E)</math> que en este caso está dada por:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\phi }'(\boldsymbol {x},\boldsymbol {B})=\boldsymbol {I}-\boldsymbol {B}^{-1}\boldsymbol {F}'(\boldsymbol {x}), </math>
|}
|}

de donde se tiene que cumplir

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\vert \boldsymbol {I}-\boldsymbol {B}^{-1}\boldsymbol {F}'(\boldsymbol {x})\right\vert =r_*<1. </math>
|}
|}

Si se toma <math display="inline">\boldsymbol {B}_*=\boldsymbol {F}'(\boldsymbol {x}^*)</math>, se tiene <math display="inline">r_*=0</math>.

<span id='theorem-T4.1'></span>Teorema 3:      Suponga que &nbsp;<math display="inline">\boldsymbol {\phi },\boldsymbol {x}^*,E_*</math> y la sucesión <math display="inline">\{ \boldsymbol {x}^k\} </math> satisfacen las hipótesis del Teorema [[#theorem-T3.1|2]]. Suponga además que <math display="inline">\boldsymbol {x}^k\neq \boldsymbol {x}^*</math> para todo <math display="inline">k=0,1,2,\ldots </math> y que

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lim _{k\to \infty }\frac{\left\vert \left[(\boldsymbol {I}-\boldsymbol {\phi }'(\boldsymbol {x}^*,\boldsymbol {E}_k))^{-1}-(\boldsymbol {I}-\boldsymbol {\phi }'(\boldsymbol {x}^*,E_*))^{-1}\right](\boldsymbol {x}^{k+1}-\boldsymbol {x}^k)\right\vert }{\left\vert \boldsymbol {x}^{k+1}-\boldsymbol {x}^k\right\vert }=0.     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

Entonces

<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline{\lim }\frac{\left\vert \boldsymbol {x}^{k+1}-\boldsymbol {x}^*\right\vert }{\left\vert \boldsymbol {x}^k-\boldsymbol {x}^*\right\vert }\leq r_*.     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

<span id='theorem-C4.1'></span>Corolario 1:      Bajo las hipótesis del Teorema [[#theorem-T4.1|3]], si <math display="inline">\displaystyle \lim _{k\to \infty }\left\Vert E^k-E^*\right\Vert =0</math>, entonces

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline{\lim }\frac{\left\vert \boldsymbol {x}^{k+1}-\boldsymbol {x}^*\right\vert }{\left\vert \boldsymbol {x}^k-\boldsymbol {x}^*\right\vert }\leq r_*.     </math>
|}
|}

<span id='theorem-T4.3'></span>Teorema 4:      Suponga que <math display="inline">\boldsymbol {\phi },\boldsymbol {x}^*,E_*</math> satisfacen las Suposiciones A1 y A2, y que <math display="inline">r^*=0</math>. Suponga que una sucesión es generada por [[#eq-35|(35)]] y que, para todo <math display="inline">k=0,1,2,\ldots </math>

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\vert \boldsymbol {\phi }'(\boldsymbol {x}^*,\boldsymbol {E}_k)-\boldsymbol {\phi }'(\boldsymbol {x}^*,E_*)\right\vert \leq M\left\vert \boldsymbol {x}^k-\boldsymbol {x}^*\right\vert ^p,     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

para algún <math display="inline">M>0</math>. Entonces, existe <math display="inline">\varepsilon{>0}</math> tal que si <math display="inline">\left\vert \boldsymbol {x}^0-\boldsymbol {x}^*\right\vert \leq \varepsilon </math>, la sucesión <math display="inline">\{ \boldsymbol {x}^k\} </math> está bien definida, converge a <math display="inline">\boldsymbol {x}^*</math> y satisface

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\vert \boldsymbol {x}^{k+1}-\boldsymbol {x}^*\right\vert \leq (L+M)\left\vert \boldsymbol {x}^k-\boldsymbol {x}^*\right\vert ^{p+1}.     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

Habíamos visto que en el caso de los métodos cuasi-Newton, <math display="inline">\boldsymbol {\phi }'(\boldsymbol {x},\boldsymbol {B})=\boldsymbol {I}-\boldsymbol {B}^{-1}\boldsymbol {F}'(\boldsymbol {x})</math>. Así,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\phi }'(\boldsymbol {x}^*,\boldsymbol {B} _k)=\boldsymbol {I}-\boldsymbol {B} _k^{-1}\boldsymbol {F}'(\boldsymbol {x}^*) </math>
|}
|}

y

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[\boldsymbol {I}-\boldsymbol {\phi }'(\boldsymbol {x}^*,\boldsymbol {B} _k)\right]^{-1} =\left[\boldsymbol {F}'(\boldsymbol {x}^*)\right]^{-1}\boldsymbol {B} _k,</math>
|-
| style="text-align: center;" | <math>     \left[\boldsymbol {I}-\boldsymbol {\phi }'(\boldsymbol {x}^*,\boldsymbol {B}_*)\right]^{-1} =\left[\boldsymbol {F}'(\boldsymbol {x}^*)\right]^{-1}\boldsymbol {B}_*. </math>
|}
|}

Por lo tanto, la condición [[#eq-45|(45)]] se convierte en

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lim _{k\to \infty }\frac{\left\vert \boldsymbol {F}'(\boldsymbol {x}^*)(\boldsymbol {B} _k-\boldsymbol {B}_*)(\boldsymbol {x}^{k+1}-\boldsymbol {x}^k)\right\vert }{\left\vert \boldsymbol {x}^{k+1}-\boldsymbol {x}^k\right\vert }=0. </math>
|}
|}

Esto es equivalente a

<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lim _{k\to \infty }\frac{\left\vert (\boldsymbol {B} _k-\boldsymbol {B}_*)(\boldsymbol {x}^{k+1}-\boldsymbol {x}^k)\right\vert }{\left\vert \boldsymbol {x}^{k+1}-\boldsymbol {x}^k\right\vert }=0. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

Esta es la condición de Dennis-Walker de convergencia de los métodos cuasi-Newton a una tasa ideal (ver <span id='citeF-8'></span>[[#cite-8|[8]]]). Ya que el método MCN es un método cuasi-Newton, se deduce su convergencia. Más aún, la convergencia del MCN está dada por ([[#eq-49|49]]).

==6 Ejemplos numéricos==

En esta sección, se compara numéricamente el método MCN con los métodos DFP y BFGS para dos tipos de funciones. Para ello, se implementaron los tres métodos en el lenguaje de programación MATLAB. Se considera dos casos: con error y sin error en los datos de la función. Para la función del ejemplo 2, el MCN obtuvo mejores resultados que los otros dos, lo que justifica la propuesta.

===6.1 Caso sin error===

Se consideró la siguiente función para la primera prueba

<span id="eq-50"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {f}: \mathbb{R}^2 \to \mathbb{R},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|-
| style="text-align: center;" | <math>     (x_1,x_2) \longmapsto x_1^4+(x_1+x_2)^2+\left(e^{x_2}-1\right)^2, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}
|}

con el punto inicial <math display="inline">\boldsymbol {x}^0=(-10,17)</math>, una tolerancia de <math display="inline">10^{-8}</math> y la matriz identidad como matriz inicial. Se obtuvieron los  resultados mostrados en la Figura [[#img-2|2]]  (que corresponde a una ventana emergente del lenguaje de programación):

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Resultados sin error.png|420px|Resultados de los métodos.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 2:''' Resultados de los métodos.
|}

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Grafica sin error.png|300px|Gráfica de la función.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 3:''' Gráfica de la función.
|}

Se puede observar que el método DFP no converge a la solución, mientras que el MCN y el BFGS si lo hacen. Además, se observa que el método BFGS utiliza un menor número de iteraciones. Una posible explicación de este  hecho, es que la función de estudio tiene una ''buena  inclinación'' lo que hace que el Hessiano promedio  para el método BFGS arroje  mejores resultados que el método MCN. La Tabla [[#table-1|1]] muestra los resultados obtenidos usando los programas elaborados para diferentes puntos iniciales. En algunos  de estos puntos, el método DFP no converge. Notar que en los puntos en los que converge el método DFP, lo hace en menos iteraciones que el MCN.  El Método BFGS fue el que generó los mejores resultados.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Tabla. 1 Número de iteraciones y tiempo, medido en segundos, de los distintos métodos usando distintos puntos iniciales. Se observa que el método MCN converge en todos los puntos. El método DFP no converge en todos los puntos y en los que lo hace, emplea menos iteraciones que el método MCN. El método BFGS converge en todos los puntos y  en menos iteraciones que los otros dos.
|- style="border-top: 2px solid;"
| <math display="inline">\boldsymbol {x}^0</math> 
| BFGS [tiempo] 
| DFP [tiempo]  
| MCN [tiempo]
|- style="border-top: 2px solid;"
|  (-10,17) 
| 61 [10.55] 
| - 
| 76 [14.86] 
|-
| (10,-17) 
| 21 [3.04] 
| 27 [3.68] 
| 31 [4.45]
|-
| (-10,-17) 
| 17 [2.7]
| 18 [2.63] 
| 37 [5.59] 
|-
| (10,17) 
| 56 [9.92] 
| - 
| 62 [12.94]
|-
| (2,1) 
| 12 [1.75] 
| 15 [2.05] 
| 17 [2.38]
|-
| (-2,1) 
| 10 [1.63] 
| 11 [1.72] 
| 28 [4.1]
|-
| (-2,-1) 
| 11 [1.71] 
| 11 [1.69]
| 39 [5.94]
|- style="border-bottom: 2px solid;"
| (2,-1) 
| 14 [2.11] 
| 24 [3.33] 
| 26 [3.61] 

|}

A continuación veremos un ejemplo en la que la función no tiene la ''buena  inclinación'' mencionada arriba y el método  MCN obtiene mejores resultados que los otros dos.

Consideremos ahora la función

<span id="eq-52"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {f}: \mathbb{R}^2 \to \mathbb{R},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|-
| style="text-align: center;" | <math>     (x_1,x_2) \longmapsto 0.00001(x_1^4+x_2^4), </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}
|}

con el punto inicial <math display="inline">\boldsymbol {x}^0=(2,1)^T</math>, una tolerancia  de <math display="inline">10^{-8}</math> y la matriz identidad como matriz inicial. Los resultados de la implementación se muestran a continuación

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Felix sin error.png|420px|Resultados de los tres métodos utilizados. El MCN utiliza un menor número de iteraciones y consume menos tiempo de cómputo.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 4:''' Resultados de los tres métodos utilizados. El MCN utiliza un menor número de iteraciones y consume menos tiempo de cómputo.
|}

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Felix sin error grafica.png|300px|Gráfica de la función f(x,y)=0.00001(x⁴+y⁴).]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 5:''' Gráfica de la función <math>\boldsymbol {f}(x,y)=0.00001(x^4+y^4)</math>.
|}


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Tabla. 2 Número de iteraciones y tiempo, medido en segundos, de los tres métodos usando distintos puntos iniciales.
|- style="border-top: 2px solid;"
| <math display="inline">\boldsymbol {x}^0</math> 
| BFGS [tiempo] 
| DFP [tiempo]  
| MCN [tiempo]
|- style="border-top: 2px solid;"
|  (2,1) 
| 32 [4.39] 
| 526 [73] 
| 13 [2.15]
|-
| (-2,1) 
| 31 [3.95] 
| 936 [129.07] 
| 13 [1.75]
|-
| (-2,-1) 
| 32 [4.49] 
| 526 [70.42] 
| 13 [2.29]
|-
| (2,-1) 
| 31 [4.02] 
| 936 [132.79] 
| 13 [1.88]
|- style="border-bottom: 2px solid;"
| (3,1) 
| 24 [3.69] 
| 977 [137.59] 
| 15 [2.94] 

|}

===6.2 Caso con error===

Consideremos la función [[#eq-52|(52)]] del ejemplo anterior, pero le agregaremos un error, por lo que la función a minimizar ahora es

<span id="eq-54"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\tilde{f}}(x,y)=\boldsymbol {f}(x,y)+\boldsymbol {e}(x,y), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

donde <math display="inline">\boldsymbol {e}(x,y)</math> está definido por

<span id="eq-55"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {e}(x,y)=0.1\boldsymbol {f}(x,y)\sin{(x+y)}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

Se hizo uso del mismo punto inicial <math display="inline">\boldsymbol {x}^0=(2,1)^t</math>, tolerancia <math display="inline">10^{-8}</math> y matriz inicial. Los resultados obtenidos fueron los siguientes

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Resultados error simple.png|420px|Resultados de los métodos.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 6:''' Resultados de los métodos.
|}

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Grafica error simple.png|300px|Gráfica de la función ̃f(x,y)=f(x,y)+e(x,y).]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 7:''' Gráfica de la función <math>\boldsymbol {\tilde{f}}(x,y)=\boldsymbol {f}(x,y)+\boldsymbol {e}(x,y)</math>.
|}


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>Tabla. 3 Número de iteraciones y tiempo, medido en segundos, de los tres métodos usando distintos puntos iniciales.
|- style="border-top: 2px solid;"
| <math display="inline">\boldsymbol {x}^0</math> 
| BFGS [tiempo] 
| DFP [tiempo]  
| MCN [tiempo]
|- style="border-top: 2px solid;"
|  (2,1) 
| 30 [4.8] 
| - 
| 13 [2.36]
|-
| (-2,1) 
| 31 [4.81]
| - 
| 13 [2.35]
|-
| (-2,-1) 
| 29 [4.41] 
| - 
| 14 [3.09]
|-
| (2,-1) 
| 28 [4.71] 
| 372 [90.37] 
| 14 [2.41]
|- style="border-bottom: 2px solid;"
| (3,1) 
| 33 [5.42] 
| - 
| 15 [2.71] 

|}

En este caso, puede observarse que el método MCN fue que que obtuvo mejores resultados, ya que converge en aproximadamente la mitad de iteraciones   que el método BFGS, mientras que el método DFP en la mayoría de los casos no converge y en caso de converger, lo hace en número   de iteraciones mucho mayor. Además, consume aproximadamente la mitad de tiempo de cómputo.

Consideramos la misma función [[#eq-54|(54)]], donde <math display="inline">\boldsymbol {f}(x,y)</math> está definido por [[#eq-52|(52)]] y ahora

<span id="eq-56"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {e}(x,y)=0.1\boldsymbol {f}(x,y)\sin{(x^2+y^2)}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

Se consideró el punto inicial <math display="inline">\boldsymbol {x}^0=(2,1)^t</math>, una tolerancia de <math display="inline">10^{-8}</math> y la matriz identidad como matriz inicial. En la Figura [[#img-8|8]] se muestran los resultados obtenidos.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Resultados ejemplo 2.png|390px|Resultados de los métodos considerando error.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 8:''' Resultados de los métodos considerando error.
|}

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Grafica Felix con error.png|300px|Gráfica de la función considerando errores.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 9:''' Gráfica de la función considerando errores.
|}


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-4'></span>Tabla. 4 Número de iteraciones y tiempo, medido en segundos, de los tres métodos considerando el error descrito y usando distintos puntos iniciales.
|- style="border-top: 2px solid;"
| <math display="inline">\boldsymbol {x}^0</math>  
| BFGS [tiempo (seg)] 
| DFP [tiempo (seg)]  
| MCN [tiempo (seg)]
|- style="border-top: 2px solid;"
|  (2,1)  
| 32 [5.47] 
| 729 [125.09]  
| 14 [2.72]
|-
| (1,2)  
| 32 [5.53] 
| 371 [70.48]  
| 14 [3.14]
|-
| (3,1)  
| - 
| -  
| 15 [2.63] 
|-
| (-3,1)  
| 29 [4.99] 
| 218 [51] 
| 13 [3.04] 
|-
| (3,-1)  
| 29 [4.91] 
| 218 [54.7] 
| 13 [6.33] 
|- style="border-bottom: 2px solid;"
| (-3,-1)  
| - 
| - 
| 15 [2.75]

|}

De los ejemplos vistos, podemos notar que el método MCN tiene un mejor comportamiento (incluso mejor que el BFGS) en los casos donde hay presencia de errores o la función objetivo es "muy plana''  cerca del mínimo.  El tipo de error considerado tiene por objetivo ilustrar el desempeño de los métodos analizados aquí ya que puede considerarse una buena representación de los errores cometidos en la práctica.  En <span id='citeF-9'></span>[[#cite-9|[9]]], se considera que tanto la función como su gradiente tienen error acotado en la norma uniforme. Este caso, podría considerarse para probar el desempeño del MCN.

==7 Aplicación a la solución numérica de sistemas de ecuaciones diferenciales ordinarias==

Como se comentó anteriormente, los sistemas no lineales aparecen en la solución numérica de sistemas de ecuaciones diferenciales no lineales. En esta sección, vamos a considerar el caso de un sistema de la forma:

<span id="eq-57"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{d\boldsymbol {y}_1}{dt}=\boldsymbol {f}_1(t,\boldsymbol {y}_1,\boldsymbol {y}_2), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

<span id="eq-58"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{d\boldsymbol {y}_2}{dt}=\boldsymbol {f}_2(t,\boldsymbol {y}_1,\boldsymbol {y}_2), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

con condiciones iniciales

<span id="eq-59"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}_1(t_0)=y_1^0,\quad \boldsymbol {y}_2(t_0)=y_2^0. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

Las funciones <math display="inline">\boldsymbol {f}_1</math> y <math display="inline">\boldsymbol {f}_2</math> son suficientemente suaves. Supongamos que las funciones <math display="inline">\boldsymbol {f}_1</math> y <math display="inline">\boldsymbol {f}_2</math> cumplen las condiciones del Teorema de Existencia y Unicidad dado en <span id='citeF-10'></span>[[#cite-10|[10]]], por lo cual el problema de valor inicial [[#eq-57|(57)]]-[[#eq-59|(59)]] tiene una única solución en algún intervalo del eje <math display="inline">t</math> que contiene a <math display="inline">t_0</math>. Deseamos determinar valores aproximados <math display="inline">\boldsymbol {y}_1^1,\boldsymbol {y}_1^2,\ldots ,\boldsymbol {y}_1^n,\ldots </math>&nbsp;y&nbsp;<math display="inline">\boldsymbol {y}_2^1,\boldsymbol {y}_2^2,\ldots ,\boldsymbol {y}_2^n,\ldots </math> de las soluciones <math display="inline">\boldsymbol {y}_1=\boldsymbol {\phi }(t)</math>, <math display="inline">\boldsymbol {y}_2=\boldsymbol {\Phi }(t)</math> en los puntos <math display="inline">t_n=t_0+nh</math>, con <math display="inline">n=1,2,\ldots </math>.

En notación vectorial, el problema de valor inicial [[#eq-57|(57)]]-[[#eq-59|(59)]] puede ser escrito como

<span id="eq-60"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}'=\boldsymbol {f}(t,\boldsymbol {y}),\quad \boldsymbol {y}(t_0)=\boldsymbol {y}_0, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
|}

donde <math display="inline">\boldsymbol {y}</math> es un vector con componentes <math display="inline">\boldsymbol {y}_1</math> y <math display="inline">\boldsymbol {y}_2</math>, <math display="inline">\boldsymbol {f}</math> es el vector de funciones con componentes <math display="inline">\boldsymbol {f}_1</math> y <math display="inline">\boldsymbol {f}_2</math> y <math display="inline">\boldsymbol {y}_0</math> es el vector con componentes <math display="inline">y_1^0,y_2^0</math>.

Uno de los métodos que se utiliza para hallar la solución numérica del sistema [[#eq-60|(60)]], es el método de Euler que consiste en considerar la fórmula

<span id="eq-61"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}_{n+1}=\boldsymbol {y}_n+h\boldsymbol {f}_n. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
|}

Las condiciones iniciales son usadas para determinar <math display="inline">\boldsymbol {f}_0</math>, el cual es el vector tangente a la gráfica de la solución <math display="inline">\boldsymbol {y}=\boldsymbol {\Phi }(t)</math> en el plano <math display="inline">y_1y_2</math>. Nos movemos en dirección de este vector tangente por un paso de tiempo <math display="inline">h</math> con el fin de encontrar el siguiente punto <math display="inline">\boldsymbol {y}_1</math>. Después, calculamos un nuevo vector tangente <math display="inline">\boldsymbol {f}_1</math>, nos movemos a lo largo de éste por un paso de tiempo <math display="inline">h</math> para encontrar <math display="inline">\boldsymbol {y}_2</math>, y así sucesivamente.

Una variación del método de Euler puede ser obtenido de la siguiente manera. Dado que <math display="inline">\boldsymbol {y}=\boldsymbol {\Phi }(t)</math> es una solución del problema de valor inicial [[#eq-60|(60)]], obtenemos

<span id="eq-62"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\Phi }'(t)=\boldsymbol {f}(t,\boldsymbol {y}(t)). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
|}

Aproximando la derivada en [[#eq-62|(62)]] por el cociente de diferencias hacia atrás <math display="inline">\dfrac{\boldsymbol {\phi }(t_n)-\boldsymbol {\phi }(t_{n-1})}{h},</math> obteniendo

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {\phi }(t_n)-\boldsymbol {\phi }(t_{n-1})\approx h\boldsymbol {f}(t_n,\boldsymbol {y}_n)), </math>
|}
|}

o bien

<span id="eq-63"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}_n=\boldsymbol {y}_{n-1}+h\boldsymbol {f}(t_n,\boldsymbol {y}_n)). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
|}


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:Valores iniciales <math display="inline">t_0</math> y <math display="inline">\boldsymbol {y}_0</math>, longitud de paso <math display="inline">h</math> y número de pasos <math display="inline">m</math>. Aproximación <math display="inline">\boldsymbol {y}</math>. <math>n</math> desde 1 hasta <math>m</math><math display="inline">t_n=t_{n-1}+h</math>. Resolver [[#eq-63|(63)]]. 


|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-4'></span>'''Algoritmo. 4''' El Método de Euler Implícito
|}

Aunque la implementación del Algoritmo [[#algorithm-4|4]] parece sencilla, tenemos el problema de resolver el sistema de ecuaciones [[#eq-63|(63)]]. Es en este momento donde hacemos uso del método MCN. El sistema de ecuaciones podemos plantearlo como un problema de optimización haciendo

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {g}(\boldsymbol {a})=\boldsymbol {y}_{n-1}+h\boldsymbol {f}(t_n,\boldsymbol {a})-\boldsymbol {a},</math>
|}
|}

y trabajando con el problema de optimización

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\min _{\boldsymbol {a}\in \mathbb{R}^2} \frac{\left\Vert g(\boldsymbol {a})\right\Vert ^2}{2}. </math>
|}
|}

Se compara el desempeño del MCN y del BFGS en el Método de Euler Implícito (MEI). Para ello, se elaboraron programas en el sistema MATLAB para resolver problemas de valores iniciales por medio del MEI. En todos los casos se usó <math display="inline">\boldsymbol {y}_{n-1}</math> como punto inicial, la matriz identidad como matriz inicial y una tolerancia de <math display="inline">10^{-8}</math>.

Consideremos el problema dado en <span id='citeF-1'></span>[[#cite-1|[1]]], el cual involucra un sistema lineal de coeficientes constantes no homogéneo de dimensión 2.

<span id="eq-64"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{bmatrix}\boldsymbol {y}_1'\\   \boldsymbol {y}_2'   \end{bmatrix}   =   \begin{bmatrix}-2 & 1\\   1 & -2   \end{bmatrix}   \begin{bmatrix}\boldsymbol {y}_1\\   \boldsymbol {y}_2   \end{bmatrix}   +   \begin{bmatrix}2\sin{t}\\   2(\cos{t}-\sin{t})   \end{bmatrix}   ;\quad    \begin{bmatrix}\boldsymbol {y}_1(0)\\   \boldsymbol {y}_2(0)   \end{bmatrix}   =   \begin{bmatrix}2\\   3   \end{bmatrix}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
|}

<!-- comment

'''Problema 2'''

<span id=""></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{bmatrix}\boldsymbol {y}_1'\\   \boldsymbol {y}_2'   \end{bmatrix}   =   \begin{bmatrix}-2 & 1\\   998 & -999   \end{bmatrix}   \begin{bmatrix}\boldsymbol {y}_1\\   \boldsymbol {y}_2   \end{bmatrix}   +   \begin{bmatrix}2\sin{t}\\   999(\cos{t}-\sin{t})   \end{bmatrix}   ;\quad    \begin{bmatrix}\boldsymbol {y}_1(0)\\   \boldsymbol {y}_2(0)   \end{bmatrix}   =   \begin{bmatrix}2\\   3   \end{bmatrix} </math>
|}
|}

--> El problema tienen la siguiente solución exacta

<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{bmatrix}\boldsymbol {y}_1(t)\\   \boldsymbol {y}_2(t)   \end{bmatrix}   =   2e^{-t}   \begin{bmatrix}1\\   1   \end{bmatrix}   +   \begin{bmatrix}\sin{t}\\   \cos{t}   \end{bmatrix}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
|}

Se tomaron un total de 1 000 puntos de la partición del intervalo temporal <math display="inline">[0,10]</math>. En la Figura [[#img-10a|10a]], se muestran la solución exacta del sistema y la obtenida usando el método BFGS en el MEI. En la Figura [[#img-10b|10b]], se muestran la solución exacta del sistema y la obtenida usando el MCN. Obsérvese que las aproximaciones son tan cercanas a la solución, que las gráficas prácticamente están sobrepuestas.

En este ejemplo, el método MCN hizo un promedio de 3.9310 iteraciones, tardando en promedio 0.2775 segundos en cada paso para la resolución del sistema, mientras que el método BFGS hace un promedio de 3.9670 iteraciones, tardando en promedio 0.2783 segundos en cada paso para la resolución del sistema. Así, el MCN tiene un mejor desempeño que el BFGS. Nótese que el error relativo es mayor en los puntos en los que se alcanzan máximo y mínimos.

<div id='img-10a'></div>
<div id='img-10b'></div>
<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Solucion BFGS.png|600px|Solución con el método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Solucion MCN.png|600px|Solución con el método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (a) Solución con el método BFGS.
| (b) Solución con el método MCN.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 10:''' Solución numérica encontrada usando los dos métodos.
|}

<div id='img-11a'></div>
<div id='img-11b'></div>
<div id='img-11c'></div>
<div id='img-11d'></div>
<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Error BFGS.png|600px|Error absoluto del método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Error MCN.png|600px|Error absoluto del método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (a) Error absoluto del método BFGS.
| (b) Error absoluto del método MCN.
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Error relativo BFGS.png|600px|Error relativo del método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Error relativo MCN.png|600px|Error relativo del método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (c) Error relativo del método BFGS.
| (d) Error relativo del método MCN.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 11:''' Error de los métodos en cada componente de la solución.
|}

Ahora, consideremos el problema no lineal dado en <span id='citeF-1'></span>[[#cite-1|[1]]].

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}_1'=\dfrac{1}{\boldsymbol {y}_1}-\boldsymbol {y}_2\frac{e^{t^2}}{t^2}-t,\quad \boldsymbol {y}_1(1.5)=2/3,</math>
|-
| style="text-align: center;" | <math>     \boldsymbol {y}_2'=\dfrac{1}{\boldsymbol {y}_2}-e^{t^2}-2te^{-t^2},\quad \boldsymbol {y}_2(1.5)=e^{-9/4}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
|}

El problema tienen la siguiente solución exacta

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}_1(t)=\dfrac{1}{t},</math>
|-
| style="text-align: center;" | <math>     \boldsymbol {y}_2(t)=e^{-t^2}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
|}

Se tomaron un total de 2 000 puntos de la partición del intervalo temporal <math display="inline">[1.5,1.7]</math>. En la Figura [[#img-12a|12a]], se muestran la solución exacta del sistema y la obtenida usando el método BFGS en el MEI. En la Figura [[#img-12b|12b]], se muestran la solución exacta del sistema y la obtenida usando el MCN.

<div id='img-12a'></div>
<div id='img-12b'></div>
<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Solucion E2 BFGS.png|600px|Solución con el método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Solucion E2 MCN.png|600px|Solución con el método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (a) Solución con el método BFGS.
| (b) Solución con el método MCN.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 12:''' Solución numérica encontrada usando los dos métodos.
|}

<div id='img-13a'></div>
<div id='img-13b'></div>
<div id='img-13c'></div>
<div id='img-13d'></div>
<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Error E2 BFGS.png|600px|Error absoluto del método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Error E2 MCN.png|600px|Error absoluto del método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (a) Error absoluto del método BFGS.
| (b) Error absoluto del método MCN.
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Error relativo E2 BFGS.png|600px|Error relativo del método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Error relativo E2 MCN.png|600px|Error relativo del método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (c) Error relativo del método BFGS.
| (d) Error relativo del método MCN.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 13:''' Error de los métodos en cada componente de la solución.
|}

En este ejemplo, el método MCN hizo un promedio de 2.4345 iteraciones, tardando en promedio 0.1846 segundos en cada paso para la resolución del sistema, mientras que el método BFGS hace un promedio de 2.7170 iteraciones, tardando en promedio 0.1855 segundos en cada paso para la resolución del sistema. Así, el MCN tuvo un mejor desempeño que el BFGS.

Ahora, consideremos el problema no lineal

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}_1'=\dfrac{4\boldsymbol {y}_1}{t},\quad \boldsymbol {y}_1(1)=10^{-4},</math>
|-
| style="text-align: center;" | <math>     \boldsymbol {y}_2'=\dfrac{4\boldsymbol {y}_2}{t},\quad \boldsymbol {y}_2(1)=10^{-4},   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
|}

con <math display="inline">a=10^{-4}</math>. El problema tienen la siguiente solución exacta

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}_1(t)=10^{-4}t^4,</math>
|-
| style="text-align: center;" | <math>     \boldsymbol {y}_2(t)=10^{-4}t^4.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
|}

Se tomaron un total de 1 000 puntos de la partición del intervalo temporal <math display="inline">[1,2]</math>. En la Figura [[#img-14a|14a]], se muestran la solución exacta del sistema y la obtenida usando el método BFGS en el MEI. En la Figura [[#img-14b|14b]], se muestran la solución exacta del sistema y la obtenida usando el MCN.

<div id='img-14a'></div>
<div id='img-14b'></div>
<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Solución E3 BFGS.png|600px|Solución con el método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Solución E3 MCN.png|600px|Solución con el método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (a) Solución con el método BFGS.
| (b) Solución con el método MCN.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 14:''' Solución numérica encontrada usando los dos métodos.
|}

En este ejemplo, el método MCN hizo un promedio de 1 iteración, tardando en promedio 0.3006 segundos en cada paso para la resolución del sistema, mientras que el método BFGS hace un promedio de 1 iteración, tardando en promedio 0.2907 segundos en cada paso para la resolución del sistema. Obsérvese que, como en los casos anteriores, las aproximaciones son tan cercanas a la solución, que las gráficas prácticamente están sobrepuestas.

En la Figura [[#img-15|15]], se muestran los errores, tanto absolutos como relativos, cometidos por los dos métodos en cada entrada de la solución del sistema. Se observa que el desempeño es similar en ambos métodos.

<div id='img-15a'></div>
<div id='img-15b'></div>
<div id='img-15c'></div>
<div id='img-15d'></div>
<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Error absoluto E3 BFGS.png|600px|Error absoluto del método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Error absoluto E3 MCN.png|600px|Error absoluto del método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (a) Error absoluto del método BFGS.
| (b) Error absoluto del método MCN.
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Error relativo E3 BFGS.png|600px|Error relativo del método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Error relativo E3 MCN.png|600px|Error relativo del método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (c) Error relativo del método BFGS.
| (d) Error relativo del método MCN.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 15:''' Error de los métodos en cada componente de la solución.
|}

Ahora, consideremos el ejemplo dado en <span id='citeF-10'></span>[[#cite-10|[10]]] dado por un modelo de Lotka-Volterra

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol {y}_1'=\boldsymbol {y}_1-0.5\boldsymbol {y}_1\boldsymbol {y}_2,\quad \boldsymbol {y}_1(0)=2.5,</math>
|-
| style="text-align: center;" | <math>     \boldsymbol {y}_2'=-0.75\boldsymbol {y}_2+0.25\boldsymbol {y}_1\boldsymbol {y}_2,\quad \boldsymbol {y}_2(0)=1,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
|}

donde <math display="inline">\boldsymbol {y}_1</math> y <math display="inline">\boldsymbol {y}_2</math> denotan la población de presa y depredador, respectivamente. Se tomaron un total de 7 000 puntos de la partición del intervalo temporal <math display="inline">[0,30]</math>. En la Figura [[#img-16a|16a]], se muestra la solución obtenida usando el método BFGS en el MEI. En la Figura [[#img-16b|16b]], se muestra la solución obtenida usando el MCN.

<div id='img-16a'></div>
<div id='img-16b'></div>
<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Acevedo Vazquez_103892577-Solucion LV 7000 BFGS.png|600px|Solución con el método BFGS.]]
|[[Image:Draft_Acevedo Vazquez_103892577-Solucion LV 7000 MCN.png|600px|Solución con el método MCN.]]
|- style="text-align: center; font-size: 75%;"
| (a) Solución con el método BFGS.
| (b) Solución con el método MCN.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 16:''' Solución numérica encontrada usando los dos métodos.
|}

En este ejemplo, el método MCN hizo un promedio de 2.8163 iteraciones, tardando en promedio 0.4433s en cada paso para la resolución del sistema, mientras que el método BFGS hace un promedio de 2.6409 iteraciones, tardando en promedio 0.4064s en cada paso para la resolución del sistema.

==8 Conclusiones==

Se propone un método cuasi-Newton que tiene como criterio minimizar el número de condición de la matriz de actualización. Esto con la intención de manejar la sensibilidad ante errores en los datos. Se da una fórmula explícita para la actualización en el caso bidimensional. Se ilustra con ejemplos la factibilidad del método. Se distingue el caso sin error y con error. En el primero de ellos, puede observarse que, aunque todos los métodos convergen a la solución, el método DFP lo hace en un número mayor de iteraciones.  Esto es acorde con lo reportado en la literatura. Por otra parte, los métodos BFGS y MCN intercambian su rol de mejor opción, dependiendo del tipo de función que se esté minimizando. Para el Ejemplo 1, el método BFGS requiere menos iteraciones que el MCN pero para el Ejemplo 2, sucede lo contrario. Pasando ahora al caso con error, el método MCN converge con menos iteraciones en ambos casos presentados. Aún más, con algunos puntos iniciales no convergen ni el método DFP ni el BFGS. Así, se muestra que el método puede mejorar los resultados tanto en iteraciones como en convergencia para cierta clase de funciones, a saber, aquellas en las que la matriz del sistema para la actualización es mal condicionada.  

Se implementó el  MCN   en el método implícito de Euler para resolver sistema de dos ecuaciones diferenciales ordinarias no lineales. Se comparó su desempeño con el BFGS en el mismo método implícito de Euler con cuatro ejemplos encontrándose que los dos métodos obtienen  resultados similares tanto en número de iteraciones como tiempo de cómputo.

===BIBLIOGRAFÍA===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' Lambert, J. D. (1991) "Numerical Methods for Ordinary Differential Systems: The Initial Value Problem". John Wiley & Sons, Inc.

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' Davidon, W C. (1959) "VARIABLE METRIC METHOD FOR MINIMIZATION", Volume. Technical Report ANL–5990 (revised),     Argonne National Laboratory

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]''' Fletcher, R. and Powell, M. J. D. (1963) "A Rapidly Convergent Descent Method for Minimization", Volume 6. The Computer Journal 2 163-168

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]''' Nocedal, Jorge and Wright, Stephen J. (2006) "Numerical optimization". Springer, 2nd Edition

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]''' M. Bazaraa and H. Sherali and C. M. Shetty. (2006) "Unconstrained Optimization". Nonlinear Programming. John Wiley & Sons, Ltd 343-467

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' Neculai Andrei. (2020) "Nonlinear Conjugate Gradient Methods for Unconstrained Optimization", Volume. Springer International Publishing, 1st Edition

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' Martínez, José Mario. (1992) "Fixed-Point Quasi-Newton Methods", Volume 29. SIAM Journal on Numerical Analysis 5 1413-1434

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]''' J. E. Dennis and Homer F. Walker. (1981) "Convergence Theorems for Least-Change Secant Update Methods", Volume 18. Society for Industrial and Applied Mathematics. SIAM Journal on Numerical Analysis 6 949&#8211;987

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]''' Xie, Yuchen and Byrd, Richard H. and Nocedal, Jorge. (2020) "Analysis of the BFGS Method with Errors", Volume 30. SIAM Journal on Optimization 1 182-209

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]''' William E. Boyce and Richard C. DiPrima. (2000) "Elementary differential equations and boundary value problems.". Elementary differential equations and boundary value problems. John Wiley, 7th Edition

<div id="cite-11"></div>
'''[11]''' Fletcher, R. (2000) "Newton-Like Methods". Practical Methods of Optimization. John Wiley & Sons, Ltd 44-79

<!-- comment-->