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 BDF (por sus siglas en inglés Backward Differentiation Formulas) [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 ${\textstyle {\boldsymbol {f}}:\mathbb {R} ^{n}\to \mathbb {R} }$ una función suave. Consideremos el siguiente problema

$\min {\boldsymbol {f}}({\boldsymbol {x}});$ ${\hbox{ s. a}}{\boldsymbol {x}}\in \mathbb {R} ^{n}.$

(1)

El primer método cuasi-Newton (denotado por DFP) fue creado por Davidon, Fletcher y Powell [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 ${\textstyle {\boldsymbol {x}}^{k}}$. Tenemos la siguiente notación: ${\textstyle {\boldsymbol {f}}^{k}={\boldsymbol {f}}({\boldsymbol {x}}^{k})}$, ${\textstyle \nabla {\boldsymbol {f}}^{k}=\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k})}$ y ${\textstyle \nabla ^{2}{\boldsymbol {f}}^{k}=\nabla ^{2}{\boldsymbol {f}}({\boldsymbol {x}}^{k})}$ donde ${\textstyle \nabla ^{2}{\boldsymbol {f}}}$ denota el hessiano de la función ${\textstyle {\boldsymbol {f}}}$, por lo que la aproximación está dada por

${\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}},$

(2)

donde ${\textstyle {\boldsymbol {B_{k}}}}$ es una matriz simétrica definida positiva de tamaño ${\textstyle n\times n}$, que se actualizará en cada iteración. La matriz ${\textstyle {\boldsymbol {B_{K}}}}$ no es el hessiano de la función ${\textstyle {\boldsymbol {f}}}$, pero es una aproximación de ella. Notemos que ${\textstyle {\boldsymbol {m}}_{k}({\boldsymbol {0}})={\boldsymbol {f}}^{k}}$ y dado que ${\textstyle {\boldsymbol {m}}_{k}}$ es una función convexa, alcanza el mínimo en

${\boldsymbol {p}}^{k}=-{\boldsymbol {B}}_{k}^{-1}\nabla {\boldsymbol {f}}^{k},$

(3)

el cual será usado en la nueva iteración

${\boldsymbol {x}}^{k+1}={\boldsymbol {x}}^{k}+\alpha _{k}{\boldsymbol {p}}^{k},$

(4)

donde la longitud de paso ${\textstyle \alpha _{k}}$ se escoge de tal forma que cumpla las condiciones de Wolfe.

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

${\boldsymbol {B}}_{k+1}({\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{k})=\nabla {\boldsymbol {f}}^{k+1}-\nabla {\boldsymbol {f}}^{k}.$

(5)

Para simplificar la ecuación anterior se definen los vectores

${\boldsymbol {s}}^{k}={\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{k}=\alpha _{k}{\boldsymbol {p}}^{k};$

(6.a)

${\boldsymbol {y}}^{k}=\nabla {\boldsymbol {f}}^{k+1}-\nabla {\boldsymbol {f}}^{k}.$

(6.b)

Así, la ecuación (5) se reescribe como

${\boldsymbol {B}}_{k+1}{\boldsymbol {s}}^{k}={\boldsymbol {y}}^{k},$

(7)

la cual es conocida como la ecuación secante.

Del hecho de que la matriz ${\textstyle {\boldsymbol {B}}_{k+1}}$ es definida positiva y por la ecuación secante, se tiene que

$({\boldsymbol {s}}^{k})^{t}{\boldsymbol {y}}^{k}>0.$

(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 ${\textstyle {\boldsymbol {x}}^{k+1}}$ y ${\textstyle {\boldsymbol {x}}^{k}}$, pero no siempre se cumplirá cuando la función no es convexa. En este último caso, se tiene que forzar a que se cumpla (8) imponiendo las condiciones de Wolfe o las condiciones fuertes de Wolfe en ${\textstyle \alpha _{k}}$ [4].

Cuando la condición de curvatura se cumple, la ecuación secante (7) tiene infinitas soluciones, ya que hay ${\textstyle {\frac {n(n+1)}{2}}}$ grados de libertad en una matriz simétrica y la ecuación secante representa ${\textstyle n}$ condiciones y hay ${\textstyle n}$ 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 ${\boldsymbol {B}}_{k+1}$, 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 ${\textstyle {\boldsymbol {B}}_{k+1}}$ es pidiendo que, entre todas las matrices simétricas que cumplen la ecuación secante, ${\textstyle {\boldsymbol {B}}_{k+1}}$ debe ser la más cercana a la matriz actual ${\textstyle {\boldsymbol {B}}_{k}}$, en otras palabras, tenemos que resolver el siguiente subproblema

$\min _{\boldsymbol {B}}\left\Vert {\boldsymbol {B}}-{\boldsymbol {B}}_{k}\right\Vert ;$

(9.a)

${\hbox{s. a }}{\boldsymbol {B}}={\boldsymbol {B}}^{t};$ ${\boldsymbol {Bs}}^{k}={\boldsymbol {y}}^{k},$

(9.b)

donde ${\textstyle {\boldsymbol {s}}^{k}}$ y ${\textstyle {\boldsymbol {y}}^{k}}$ satisfacen la condición de curvatura y ${\textstyle {\boldsymbol {B}}_{k}}$ es una matriz simétrica definida positiva. Muchas normas pueden ser usadas en (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

$\left\Vert {\boldsymbol {A}}\right\Vert _{\boldsymbol {W}}^{2}=tr\left({\boldsymbol {WA}}^{t}{\boldsymbol {WA}}\right).$

(10)

El peso ${\textstyle {\boldsymbol {W}}}$ puede ser escogido como cualquier matriz que satisfaga la relación ${\textstyle {\boldsymbol {Wy}}^{k}={\boldsymbol {s}}^{k}}$. Se asume ${\textstyle {\boldsymbol {W}}={\overline {\boldsymbol {G}}}_{k}^{-1}}$, donde ${\textstyle {\overline {\boldsymbol {G}}}_{k}}$ es el hessiano promedio, dado por

${\overline {\boldsymbol {G}}}_{k}=\int _{0}^{1}\nabla ^{2}{\boldsymbol {f}}({\boldsymbol {x}}^{k}+\tau \alpha _{k}{\boldsymbol {p}}^{k})d\tau {.}$

(11)

Fácilmente puede verificarse que ${\textstyle {\overline {\boldsymbol {G}}}_{k}^{-1}{\boldsymbol {y}}^{k}={\boldsymbol {s}}^{k}}$.

Con la norma y matriz de peso descritas arriba, la única solución de (9) es

${\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},$

(12)

donde ${\textstyle \gamma _{k}={\dfrac {1}{({\boldsymbol {y}}^{k})^{t}{\boldsymbol {s}}^{k}}}.}$

Esta fórmula es llamada la fórmula de actualización DFP, ya que fue propuesta por Davidon en 1959 [2], y posteriormente estudiada, implementada y popularizada por Fletcher y Powell [3].

Sea ${\textstyle {\boldsymbol {H}}_{k}={\boldsymbol {B}}_{k}^{-1}}$. Usando la fórmula de Sherman-Morrison-Woodbury [4] y la ecuación (12) se sigue que

${\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}}}.$

(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 (13) son matrices de rango 1, así ${\textstyle {\boldsymbol {H}}_{k}}$ 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.

Punto inicial ${\textstyle {\boldsymbol {x}}^{0}}$, tolerancia ${\textstyle \varepsilon }$, aproximación del inverso del hessiano ${\textstyle {\boldsymbol {H}}_{0}}$. ${\textstyle {\boldsymbol {x}}^{*}}$ mínimo de ${\textstyle {\boldsymbol {f}}}$. ${\textstyle k\leftarrow 0.}$ $\left\Vert \nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k})\right\Vert >\varepsilon$Calcular la dirección de búsqueda ${\boldsymbol {p}}^{k}=-{\boldsymbol {H}}_{k}\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k}).$ (14) Calcular ${\textstyle {\boldsymbol {x}}^{k+1}={\boldsymbol {x}}^{k}+\alpha _{k}{\boldsymbol {p}}^{k}}$, donde ${\textstyle \alpha _{k}}$ es calculado mediante un procedimiento de búsqueda en la línea que satisfaga las condiciones de Wolfe, ${\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},$ (15.a) $\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},$ (15.b) con ${\textstyle 0<c_{1}<c_{2}<1.}$ Definir ${\textstyle {\boldsymbol {s}}^{k}={\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{k}}$ y ${\textstyle {\boldsymbol {y}}^{k}=\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k+1})-\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k})}$. Calcular ${\textstyle {\boldsymbol {H}}_{k+1}}$ por medio de (13). ${\textstyle k\leftarrow k+1.}$

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 (12). Las condiciones impuestas en las aproximaciones del hessiano ${\textstyle {\boldsymbol {B}}_{k}}$, pueden ser impuestas en sus inversas ${\textstyle {\boldsymbol {H}}_{k}}$. Por lo tanto, la actualización ${\textstyle {\boldsymbol {H}}_{k+1}}$ debe ser simétrica, definida positiva, y satisfacer la ecuación secante (7), ahora escrita como ${\textstyle {\boldsymbol {H}}_{k+1}{\boldsymbol {y}}^{k}={\boldsymbol {s}}^{k}.}$

La condición de la cercanía a ${\textstyle {\boldsymbol {H}}_{k}}$ es ahora especificada de forma análoga a (9)

$\min _{\boldsymbol {H}}\left\Vert {\boldsymbol {H}}-{\boldsymbol {H}}_{k}\right\Vert ;$

(16.a)

${\hbox{s. a }}{\boldsymbol {H}}={\boldsymbol {H}}^{t};$ ${\boldsymbol {Hy}}^{k}={\boldsymbol {s}}^{k}.$

(16.b)

La norma considerada es de nuevo la norma pesada de Frobenius descrita anteriormente, donde la matriz de peso ${\textstyle {\boldsymbol {W}}}$ es ahora cualquier matriz que satisfaga ${\textstyle {\boldsymbol {Ws}}^{k}={\boldsymbol {y}}^{k}}$. (Por concreción, se asume de nuevo que ${\textstyle {\boldsymbol {W}}}$ está dado por el inverso del hessiano promedio ${\textstyle {\overline {\boldsymbol {G}}}_{k}}$ definida en (11)). La única solución ${\textstyle {\boldsymbol {H}}_{k+1}}$ de (16) está dada por

${\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},$

(17)

donde

$\rho _{k}={\frac {1}{({\boldsymbol {y}}^{k})^{t}{\boldsymbol {s}}^{k}}},$

(18)

o bien, desarrollando las multiplicaciones

${\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).$

(19)

La demostración del siguiente teorema puede encontrarse en [5].

Teorema 1: Si la fórmula de actualización ${\textstyle BFGS}$ (19), es ahora escrita como ${\textstyle {\boldsymbol {H}}_{k+1}={\boldsymbol {H}}_{k}+{\boldsymbol {E}}}$, entonces ${\textstyle {\boldsymbol {E}}}$ resuelve el problema

$\min _{\boldsymbol {E}}\left\Vert {\boldsymbol {E}}\right\Vert _{\boldsymbol {W}};$ (20.a) $s.a{\boldsymbol {E}}={\boldsymbol {E}}^{t},$ (20.b) ${\boldsymbol {Ey}}^{k}={\boldsymbol {\eta }},$ (20.c)

donde ${\textstyle {\boldsymbol {W}}}$ satisface ${\textstyle {\boldsymbol {Ws}}^{k}={\boldsymbol {y}}^{k}}$  y  ${\textstyle {\boldsymbol {\eta }}={\boldsymbol {s}}^{k}-{\boldsymbol {H}}_{k}{\boldsymbol {y}}^{k}}$.

De acuerdo con Nocedal [4], si ${\textstyle {\boldsymbol {H}}_{k}}$ es definida positiva, entonces ${\textstyle {\boldsymbol {H}}_{k+1}}$ es definida positiva.

Punto inicial ${\textstyle {\boldsymbol {x}}_{0}}$, tolerancia ${\textstyle \varepsilon }$, aproximación del inverso del hessiano ${\textstyle {\boldsymbol {H}}_{0}}$. ${\textstyle {\boldsymbol {x}}^{*}}$ mínimo de ${\textstyle {\boldsymbol {f}}}$. ${\textstyle k\leftarrow 0.}$ $\left\Vert \nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k})\right\Vert >\varepsilon$Calcular la dirección de búsqueda ${\boldsymbol {p}}^{k}=-{\boldsymbol {H}}_{k}\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k}).$ (21) Calcular ${\textstyle {\boldsymbol {x}}^{k+1}={\boldsymbol {x}}^{k}+\alpha _{k}{\boldsymbol {p}}^{k}}$, donde ${\textstyle \alpha _{k}}$ es calculado mediante un procedimiento de búsqueda en la línea que satisfaga las condiciones de Wolfe (1). Definir ${\textstyle {\boldsymbol {s}}^{k}={\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{k}}$ y ${\textstyle {\boldsymbol {y}}^{k}=\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k+1})-\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k})}$. Calcular ${\textstyle {\boldsymbol {H}}_{k+1}}$ por medio de (19). ${\textstyle k\leftarrow k+1.}$

Algoritmo. 2 Algoritmo BFGS

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 [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 ${\textstyle {\boldsymbol {\bar {f}}}:\mathbb {R} ^{2}\to \mathbb {R} }$ definida por

${\boldsymbol {\bar {f}}}(x,y)={\boldsymbol {f}}(x,y)+{\boldsymbol {e}}(x,y),$

(22)

donde ${\textstyle {\boldsymbol {f}}(x,y)=0.00001(x^{4}+y^{4})}$ y ${\textstyle {\boldsymbol {e}}(x,y)=0.1{\boldsymbol {f}}(x,y)sin(x^{2}+y^{2})}$, el punto inicial ${\textstyle {\boldsymbol {x}}^{0}=(3,1)^{t}}$, una tolerancia de ${\textstyle 10^{-8}}$ y la matriz identidad como matriz inicial. En este ejemplo, ${\textstyle {\boldsymbol {e}}(x,y)}$ 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

Figura 1: Resultados de los métodos DFP y BFGS para el punto inicial ${\boldsymbol {x}}^{0}=(3,1)^{t}$.

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 ${\textstyle {\boldsymbol {B}}_{k+1}}$ 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 ${\textstyle {\boldsymbol {B}}_{k+1}}$. 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 (9) se resuelve

$\min _{\boldsymbol {B}}cond({\boldsymbol {B}});$

(23.a)

${\hbox{s. a }}{\boldsymbol {B}}={\boldsymbol {B}}^{t};$ ${\boldsymbol {{B}s}}^{k}={\boldsymbol {y}}^{k},$

(23.b)

donde ${\textstyle cond({\boldsymbol {B}})}$ está dado por ${\textstyle {\frac {\lambda _{n}}{\lambda _{1}}}}$, con ${\textstyle 0<\lambda _{1}<\cdots {<\lambda }_{n}}$ valores propios de la matriz ${\textstyle {\boldsymbol {B}}}$ (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 ${\textstyle {\boldsymbol {f}}}$ es una función de ${\textstyle \mathbb {R} ^{2}}$ a ${\textstyle \mathbb {R} }$. Consideramos ${\textstyle {\boldsymbol {y}}^{k}=(y_{1}^{k},y_{2}^{k})}$, ${\textstyle {\boldsymbol {s}}^{k}=(s_{1}^{k},s_{2}^{k})}$ y

${\boldsymbol {B}}={\begin{bmatrix}b_{11}&b_{12}\\b_{12}&b_{22}\end{bmatrix}}.$

(24)

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

$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}})}}}},$

(25)

donde ${\textstyle tr({\boldsymbol {B}})}$ y ${\textstyle det({\boldsymbol {B}})}$ denotan la traza y el determinante de ${\textstyle {\boldsymbol {B}}}$, respectivamente. Ahora bien, dado que la matriz tiene que ser simétrica y cumplir la ecuación secante, tenemos que

$b_{11}s_{1}^{k}+b_{12}s_{2}^{k}=y_{1}^{k},$

$b_{12}s_{1}^{k}+b_{22}s_{2}^{k}=y_{2}^{k},$

de donde se obtiene que

$b_{11}={\frac {y_{1}^{k}-b_{12}s_{2}^{k}}{s_{1}^{k}}},$

$b_{22}={\frac {y_{2}^{k}-b_{12}s_{1}^{k}}{s_{2}^{k}}},$

$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),$

$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}.$

Sustituyendo las ecuaciones anteriores en (25), se obtiene una función de ${\textstyle \mathbb {R} }$ a ${\textstyle \mathbb {R} }$, la cual se minimiza calculando la primera derivada (respecto a ${\textstyle b_{12}}$) e igualándola a 0. Con lo cual se obtiene que

$b_{12}={\frac {y_{1}^{k}y_{2}^{k}\left\Vert {\boldsymbol {s}}^{k}\right\Vert ^{2}-s_{1}^{k}s_{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}}},$

que al sustituirlo en las expresiones de ${\textstyle b_{11}}$ y ${\textstyle b_{22}}$ se llega a

$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}}}.$

Con lo cual se obtienen las entradas de la matriz a la cual llamaremos ${\textstyle {\boldsymbol {B}}_{k+1}^{MCN}}$, ya que es la que tiene el menor número de condición (minor condition number), la cual resuelve el subproblema (23). Además, sustituyendo el valor de ${\textstyle b_{12}}$ en la expresión del determinante, se obtiene que ${\textstyle {\boldsymbol {B}}_{k+1}^{MCN}}$ es definida positiva. Haciendo los cálculos y simplificaciones necesarias se sigue que

${\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}}},$

(29)

donde ${\textstyle ({\boldsymbol {z}}^{k})^{t}=(({\boldsymbol {s}}^{k})^{\perp })^{t}=(s_{2}^{k},-s_{1}^{k})}$. 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 ${\textstyle {\boldsymbol {H}}_{k}={\boldsymbol {B}}_{k}^{-1}}$, teniendo así una expresión del inverso de (29)

${\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}}},$

(30)

donde ${\textstyle ({\boldsymbol {w}}^{k})^{t}=(({\boldsymbol {y}}^{k})^{\perp })^{t}=(y_{2}^{k},-y_{1}^{k})}$.

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 ${\textstyle {\boldsymbol {B}}_{k+1}^{MCN}}$ sea la matriz que resuelva el problema (9) con la misma norma ${\textstyle \left\Vert {\boldsymbol {A}}\right\Vert _{\boldsymbol {W}}=tr({\boldsymbol {WA}}^{t}{\boldsymbol {WA}})}$, pero con una matriz de peso adecuada. Así, el problema ahora es encontrar ${\textstyle {\boldsymbol {W}}}$ de tal forma que ${\textstyle {\boldsymbol {Wy}}^{k}={\boldsymbol {s}}^{k}}$ y

${\boldsymbol {B}}_{k+1}^{MCN}=\arg \min _{\boldsymbol {B}}\left\Vert {\boldsymbol {B}}-{\boldsymbol {B}}_{k}\right\Vert _{\boldsymbol {W}};$ ${\hbox{ s. a }}{\boldsymbol {B}}={\boldsymbol {B}}^{t};$ ${\boldsymbol {{B}s}}^{k}={\boldsymbol {y}}^{k}.$

Notando que ${\textstyle {\boldsymbol {B}}}$ y ${\textstyle {\boldsymbol {B}}_{k}}$ son matrices simétricas se tiene que

$\left\Vert {\boldsymbol {B}}-{\boldsymbol {B}}_{k}\right\Vert _{\boldsymbol {W}}=tr\left(\left({\boldsymbol {W}}({\boldsymbol {B}}-{\boldsymbol {B}}_{k})\right)^{2}\right).$

Además, de la restricción ${\textstyle {\boldsymbol {Wy}}^{k}={\boldsymbol {s}}^{k}}$ y haciendo un análisis análogo al realizado para ${\textstyle {\boldsymbol {B}}_{k+1}^{MCN}}$ se obtiene que

$w_{11}={\frac {s_{1}^{k}-w_{12}y_{2}^{k}}{y_{1}^{k}}},$ (31) $w_{22}={\frac {s_{2}^{k}-w_{12}y_{1}^{k}}{y_{2}^{k}}}.$ (32)

Si se denota por ${\textstyle c_{1}=b_{11}-b_{11}^{k}}$, ${\textstyle c_{2}=b_{12}-b_{12}^{k}}$ y ${\textstyle c_{3}=b_{22}-b_{22}^{k}}$, se tiene que

${\boldsymbol {g}}(w_{12})=\left\Vert {\boldsymbol {B}}-{\boldsymbol {B}}_{k}\right\Vert _{\boldsymbol {W}}$ $=2\left(c_{1}w_{12}+{\frac {c_{2}(s_{2}^{k}-w_{12}y_{1}^{k})}{y_{2}^{k}}}\right)\left(c_{3}w_{12}+{\frac {c_{2}(s_{1}^{k}-w_{12}y_{2}^{k})}{y_{1}^{k}}}\right)$ $+\left(c_{2}w_{12}+{\frac {c_{1}(s_{1}^{k}-w_{12}y_{2}^{k})}{y_{1}^{k}}}\right)^{2}+\left(c_{2}w_{12}+{\frac {c_{3}(s_{2}^{k}-w_{12}y_{1}^{k})}{y_{2}^{k}}}\right)^{2}.$

Obteniendo los puntos estacionarios de la expresión anterior (con respecto a ${\textstyle b_{ij}}$), obtenemos que

$w_{12}={\frac {b_{22}^{k}(s_{2}^{k})^{2}(y_{1}^{k})^{2}-2b_{12}^{k}(s_{2}^{k})^{2}(y_{1}^{k})(y_{2}^{k})+b_{11}^{k}(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}}}$

(33)

Sustituyendo (33) en (31) y (32), se halla a la matriz de peso ${\textstyle {\boldsymbol {W}}}$.

4 Algoritmo

Definimos ahora, un algoritmo que encuentre la solución de (1)

Punto inicial ${\textstyle {\boldsymbol {x}}^{0}}$, tolerancia ${\textstyle \varepsilon }$, aproximación del inverso del hessiano ${\textstyle {\boldsymbol {H}}_{0}}$. ${\textstyle {\boldsymbol {x}}^{*}}$ mínimo de ${\textstyle {\boldsymbol {f}}}$. ${\textstyle k\leftarrow 0.}$ $\left\Vert \nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k})\right\Vert >\varepsilon$Calcular la dirección de búsqueda ${\boldsymbol {p}}^{k}=-{\boldsymbol {H}}_{k}\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k}).$ (34) Calcular ${\textstyle {\boldsymbol {x}}^{k+1}={\boldsymbol {x}}^{k}+\alpha _{k}{\boldsymbol {p}}^{k}}$, donde ${\textstyle \alpha _{k}}$ es calculado mediante un procedimiento de búsqueda en la línea que satisfaga las condiciones de Wolfe (1). Definir ${\textstyle {\boldsymbol {s}}^{k}={\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{k}}$ y ${\textstyle {\boldsymbol {y}}^{k}=\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k+1})-\nabla {\boldsymbol {f}}({\boldsymbol {x}}^{k})}$. Calcular ${\textstyle {\boldsymbol {H}}_{k+1}}$ por medio de (30). ${\textstyle k\leftarrow k+1.}$

Algoritmo. 3 Algoritmo MCN

5 Convergencia

Para analizar la convergencia del método MCN, usaremos los resultados dados en [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 [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 ${\textstyle \nabla {\boldsymbol {f}}({\boldsymbol {x}})={\boldsymbol {0}}}$.

Sean ${\textstyle X}$ un espacio lineal (vectorial) de dimensión finita, ${\textstyle D\subset X}$, ${\textstyle \Omega \subset \mathbb {R} ^{n}}$, ${\textstyle {\boldsymbol {\phi }}:\Omega \times D\to \mathbb {R} ^{n}}$. Se consideran los métodos iterativos definidos por

${\boldsymbol {x}}^{k+1}={\boldsymbol {\phi }}({\boldsymbol {x}}^{k},E_{k}).$

(35)

Definición 1: Decimos que ${\textstyle {\boldsymbol {x}}\in \Omega }$ es un punto fijo de ${\textstyle {\boldsymbol {\phi }}}$ si y sólo si,

${\boldsymbol {x}}={\boldsymbol {\phi }}({\boldsymbol {x}},E),\quad {\hbox{para todo }}E\in D.$

(36)

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

${\boldsymbol {\phi }}({\boldsymbol {x}},{\boldsymbol {B}})={\boldsymbol {x}}-{\boldsymbol {B}}^{-1}{\boldsymbol {F}}({\boldsymbol {x}}).$

(37)

Claramente, el conjunto de puntos ${\textstyle {\boldsymbol {x}}\in \mathbb {R} ^{n}}$ tales que ${\textstyle {\boldsymbol {x}}={\boldsymbol {\phi }}({\boldsymbol {x}},{\boldsymbol {B}})}$ para todo ${\textstyle {\boldsymbol {B}}\in D}$, es el conjunto de soluciones de ${\textstyle {\boldsymbol {F}}({\boldsymbol {x}})={\boldsymbol {0}}}$. Los métodos definidos por ${\textstyle {\boldsymbol {x}}^{k+1}={\boldsymbol {\phi }}({\boldsymbol {x}}^{k},{\boldsymbol {B}}_{k})}$ 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 (35) esté bien definida y converja a un punto fijo de ${\textstyle {\boldsymbol {\phi }}}$ con una convergencia lineal. Se consideran dos suposiciones, A1 y A2 que se establecen a continuación.

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

1. Existe ${\textstyle {\boldsymbol {x}}^{*}\in \Omega }$, ${\textstyle L,p>0}$ tales que

   $\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},$

   (38)

para todo ${\textstyle {\boldsymbol {x}}\in \Omega }$, ${\textstyle E\in D}$. Esto implica que para todo ${\textstyle {\boldsymbol {x}},{\boldsymbol {z}}\in \Omega }$, ${\textstyle E\in D}$,

$\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},$

(39)

donde

$\sigma ({\boldsymbol {x}},{\boldsymbol {z}})=\max \left\{\left\vert {\boldsymbol {x}}-{\boldsymbol {x}}^{*}\right\vert ,\left\vert {\boldsymbol {z}}-{\boldsymbol {x}}^{*}\right\vert \right\}.$

(40)

2. Para todo ${\textstyle E\in D}$,

   ${\boldsymbol {x}}^{*}={\boldsymbol {\phi }}({\boldsymbol {x}}^{*},E).$

   (41)

Suposición A2. Existe ${\textstyle E_{*}\in D}$ tal que

$\left\vert {\boldsymbol {\phi }}'({\boldsymbol {x}}^{*},E_{*})\right\vert =r^{*}<1.$

(42)

Lema 1: Suponga que  ${\textstyle {\boldsymbol {\phi }},{\boldsymbol {x}}^{*},E_{*}}$ satisfacen las Suposiciones A1 y A2. Sea ${\textstyle r\in (r^{*},1)}$. Entonces, existen ${\textstyle {\overline {\varepsilon }}={\overline {\varepsilon }}(r),{\overline {\delta }}={\overline {\delta }}(r)}$ tales que

$\left\vert {\boldsymbol {\phi }}({\boldsymbol {x}},E)-{\boldsymbol {x}}^{*}\right\vert \leq \left\vert {\boldsymbol {x}}-{\boldsymbol {x}}^{*}\right\vert ,$

(43)

siempre que ${\textstyle \left\vert {\boldsymbol {x}}-{\boldsymbol {x}}^{*}\right\vert \leq {\overline {\varepsilon }}}$, ${\textstyle \left\Vert E-E_{*}\right\Vert \leq {\overline {\delta }}}$.

Teorema 2: Suponga que ${\textstyle {\boldsymbol {\phi }},{\boldsymbol {x}}^{*},E_{*}}$ satisfacen las Suposiciones A1 y A2. Sea ${\textstyle r\in (r^{*},1)}$. Entonces, existen ${\textstyle \varepsilon =\varepsilon (r),\delta =\delta (r)}$ tales que, si ${\textstyle \left\vert {\boldsymbol {x}}^{0}-{\boldsymbol {x}}^{*}\right\vert \leq \varepsilon }$ y ${\textstyle \left\Vert E_{k}-E_{*}\right\Vert \leq \delta }$ para todo ${\textstyle k=0,1,\ldots }$, entonces la sucesión generada por (35) está bien definida, converge a ${\textstyle {\boldsymbol {x}}^{*}}$ y satisface

$\left\vert {\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{*}\right\vert \leq \left\vert {\boldsymbol {x}}^{k}-{\boldsymbol {x}}^{*}\right\vert ,$

(44)

para todo ${\textstyle k=0,1,2,\ldots }$.

Para el caso de los métodos cuasi-Newton, a fin de verificar la condición (42), se debe calcular ${\textstyle {\boldsymbol {\phi }}'({\boldsymbol {x}},E)}$ que en este caso está dada por:

${\boldsymbol {\phi }}'({\boldsymbol {x}},{\boldsymbol {B}})={\boldsymbol {I}}-{\boldsymbol {B}}^{-1}{\boldsymbol {F}}'({\boldsymbol {x}}),$

de donde se tiene que cumplir

$\left\vert {\boldsymbol {I}}-{\boldsymbol {B}}^{-1}{\boldsymbol {F}}'({\boldsymbol {x}})\right\vert =r_{*}<1.$

Si se toma ${\textstyle {\boldsymbol {B}}_{*}={\boldsymbol {F}}'({\boldsymbol {x}}^{*})}$, se tiene ${\textstyle r_{*}=0}$.

Teorema 3: Suponga que  ${\textstyle {\boldsymbol {\phi }},{\boldsymbol {x}}^{*},E_{*}}$ y la sucesión ${\textstyle \{{\boldsymbol {x}}^{k}\}}$ satisfacen las hipótesis del Teorema 2. Suponga además que ${\textstyle {\boldsymbol {x}}^{k}\neq {\boldsymbol {x}}^{*}}$ para todo ${\textstyle k=0,1,2,\ldots }$ y que

$\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.$

(45)

Entonces

${\overline {\lim }}{\frac {\left\vert {\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{*}\right\vert }{\left\vert {\boldsymbol {x}}^{k}-{\boldsymbol {x}}^{*}\right\vert }}\leq r_{*}.$

(46)

Corolario 1: Bajo las hipótesis del Teorema 3, si ${\textstyle \displaystyle \lim _{k\to \infty }\left\Vert E^{k}-E^{*}\right\Vert =0}$, entonces

${\overline {\lim }}{\frac {\left\vert {\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{*}\right\vert }{\left\vert {\boldsymbol {x}}^{k}-{\boldsymbol {x}}^{*}\right\vert }}\leq r_{*}.$

Teorema 4: Suponga que ${\textstyle {\boldsymbol {\phi }},{\boldsymbol {x}}^{*},E_{*}}$ satisfacen las Suposiciones A1 y A2, y que ${\textstyle r^{*}=0}$. Suponga que una sucesión es generada por (35) y que, para todo ${\textstyle k=0,1,2,\ldots }$

$\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},$

(47)

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

$\left\vert {\boldsymbol {x}}^{k+1}-{\boldsymbol {x}}^{*}\right\vert \leq (L+M)\left\vert {\boldsymbol {x}}^{k}-{\boldsymbol {x}}^{*}\right\vert ^{p+1}.$

(48)

Habíamos visto que en el caso de los métodos cuasi-Newton, ${\textstyle {\boldsymbol {\phi }}'({\boldsymbol {x}},{\boldsymbol {B}})={\boldsymbol {I}}-{\boldsymbol {B}}^{-1}{\boldsymbol {F}}'({\boldsymbol {x}})}$. Así,

${\boldsymbol {\phi }}'({\boldsymbol {x}}^{*},{\boldsymbol {B}}_{k})={\boldsymbol {I}}-{\boldsymbol {B}}_{k}^{-1}{\boldsymbol {F}}'({\boldsymbol {x}}^{*})$

y

$\left[{\boldsymbol {I}}-{\boldsymbol {\phi }}'({\boldsymbol {x}}^{*},{\boldsymbol {B}}_{k})\right]^{-1}=\left[{\boldsymbol {F}}'({\boldsymbol {x}}^{*})\right]^{-1}{\boldsymbol {B}}_{k},$ $\left[{\boldsymbol {I}}-{\boldsymbol {\phi }}'({\boldsymbol {x}}^{*},{\boldsymbol {B}}_{*})\right]^{-1}=\left[{\boldsymbol {F}}'({\boldsymbol {x}}^{*})\right]^{-1}{\boldsymbol {B}}_{*}.$

Por lo tanto, la condición (45) se convierte en

$\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.$

Esto es equivalente a

$\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.$

(49)

Esta es la condición de Dennis-Walker de convergencia de los métodos cuasi-Newton a una tasa ideal (ver [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 (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

${\boldsymbol {f}}:\mathbb {R} ^{2}\to \mathbb {R} ,$ (50) $(x_{1},x_{2})\longmapsto x_{1}^{4}+(x_{1}+x_{2})^{2}+\left(e^{x_{2}}-1\right)^{2},$ (51)

con el punto inicial ${\textstyle {\boldsymbol {x}}^{0}=(-10,17)}$, una tolerancia de ${\textstyle 10^{-8}}$ y la matriz identidad como matriz inicial. Se obtuvieron los resultados mostrados en la Figura 2 (que corresponde a una ventana emergente del lenguaje de programación):

Figura 2: Resultados de los métodos.

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 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.

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.

${\textstyle {\boldsymbol {x}}^{0}}$	BFGS [tiempo]	DFP [tiempo]	MCN [tiempo]
(-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]
(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

${\boldsymbol {f}}:\mathbb {R} ^{2}\to \mathbb {R} ,$ (52) $(x_{1},x_{2})\longmapsto 0.00001(x_{1}^{4}+x_{2}^{4}),$ (53)

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

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.

Figura 5: Gráfica de la función ${\boldsymbol {f}}(x,y)=0.00001(x^{4}+y^{4})$.

Tabla. 2 Número de iteraciones y tiempo, medido en segundos, de los tres métodos usando distintos puntos iniciales.

${\textstyle {\boldsymbol {x}}^{0}}$	BFGS [tiempo]	DFP [tiempo]	MCN [tiempo]
(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]
(3,1)	24 [3.69]	977 [137.59]	15 [2.94]

6.2 Caso con error

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

${\boldsymbol {\tilde {f}}}(x,y)={\boldsymbol {f}}(x,y)+{\boldsymbol {e}}(x,y),$

(54)

donde ${\textstyle {\boldsymbol {e}}(x,y)}$ está definido por

${\boldsymbol {e}}(x,y)=0.1{\boldsymbol {f}}(x,y)\sin {(x+y)}.$

(55)

Se hizo uso del mismo punto inicial ${\textstyle {\boldsymbol {x}}^{0}=(2,1)^{t}}$, tolerancia ${\textstyle 10^{-8}}$ y matriz inicial. Los resultados obtenidos fueron los siguientes

Figura 6: Resultados de los métodos.

Figura 7: Gráfica de la función ${\boldsymbol {\tilde {f}}}(x,y)={\boldsymbol {f}}(x,y)+{\boldsymbol {e}}(x,y)$.

Tabla. 3 Número de iteraciones y tiempo, medido en segundos, de los tres métodos usando distintos puntos iniciales.

${\textstyle {\boldsymbol {x}}^{0}}$	BFGS [tiempo]	DFP [tiempo]	MCN [tiempo]
(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]
(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 (54), donde ${\textstyle {\boldsymbol {f}}(x,y)}$ está definido por (52) y ahora

${\boldsymbol {e}}(x,y)=0.1{\boldsymbol {f}}(x,y)\sin {(x^{2}+y^{2})}.$

(56)

Se consideró el punto inicial ${\textstyle {\boldsymbol {x}}^{0}=(2,1)^{t}}$, una tolerancia de ${\textstyle 10^{-8}}$ y la matriz identidad como matriz inicial. En la Figura 8 se muestran los resultados obtenidos.

Figura 8: Resultados de los métodos considerando error.

Figura 9: Gráfica de la función considerando errores.

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.

${\textstyle {\boldsymbol {x}}^{0}}$	BFGS [tiempo (seg)]	DFP [tiempo (seg)]	MCN [tiempo (seg)]
(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]
(-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 [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:

${\frac {d{\boldsymbol {y}}_{1}}{dt}}={\boldsymbol {f}}_{1}(t,{\boldsymbol {y}}_{1},{\boldsymbol {y}}_{2}),$

(57)

${\frac {d{\boldsymbol {y}}_{2}}{dt}}={\boldsymbol {f}}_{2}(t,{\boldsymbol {y}}_{1},{\boldsymbol {y}}_{2}),$

(58)

con condiciones iniciales

${\boldsymbol {y}}_{1}(t_{0})=y_{1}^{0},\quad {\boldsymbol {y}}_{2}(t_{0})=y_{2}^{0}.$

(59)

Las funciones ${\textstyle {\boldsymbol {f}}_{1}}$ y ${\textstyle {\boldsymbol {f}}_{2}}$ son suficientemente suaves. Supongamos que las funciones ${\textstyle {\boldsymbol {f}}_{1}}$ y ${\textstyle {\boldsymbol {f}}_{2}}$ cumplen las condiciones del Teorema de Existencia y Unicidad dado en [10], por lo cual el problema de valor inicial (57)-(59) tiene una única solución en algún intervalo del eje ${\textstyle t}$ que contiene a ${\textstyle t_{0}}$. Deseamos determinar valores aproximados ${\textstyle {\boldsymbol {y}}_{1}^{1},{\boldsymbol {y}}_{1}^{2},\ldots ,{\boldsymbol {y}}_{1}^{n},\ldots }$ y ${\textstyle {\boldsymbol {y}}_{2}^{1},{\boldsymbol {y}}_{2}^{2},\ldots ,{\boldsymbol {y}}_{2}^{n},\ldots }$ de las soluciones ${\textstyle {\boldsymbol {y}}_{1}={\boldsymbol {\phi }}(t)}$, ${\textstyle {\boldsymbol {y}}_{2}={\boldsymbol {\Phi }}(t)}$ en los puntos ${\textstyle t_{n}=t_{0}+nh}$, con ${\textstyle n=1,2,\ldots }$.

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

${\boldsymbol {y}}'={\boldsymbol {f}}(t,{\boldsymbol {y}}),\quad {\boldsymbol {y}}(t_{0})={\boldsymbol {y}}_{0},$

(60)

donde ${\textstyle {\boldsymbol {y}}}$ es un vector con componentes ${\textstyle {\boldsymbol {y}}_{1}}$ y ${\textstyle {\boldsymbol {y}}_{2}}$, ${\textstyle {\boldsymbol {f}}}$ es el vector de funciones con componentes ${\textstyle {\boldsymbol {f}}_{1}}$ y ${\textstyle {\boldsymbol {f}}_{2}}$ y ${\textstyle {\boldsymbol {y}}_{0}}$ es el vector con componentes ${\textstyle y_{1}^{0},y_{2}^{0}}$.

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

${\boldsymbol {y}}_{n+1}={\boldsymbol {y}}_{n}+h{\boldsymbol {f}}_{n}.$

(61)

Las condiciones iniciales son usadas para determinar ${\textstyle {\boldsymbol {f}}_{0}}$, el cual es el vector tangente a la gráfica de la solución ${\textstyle {\boldsymbol {y}}={\boldsymbol {\Phi }}(t)}$ en el plano ${\textstyle y_{1}y_{2}}$. Nos movemos en dirección de este vector tangente por un paso de tiempo ${\textstyle h}$ con el fin de encontrar el siguiente punto ${\textstyle {\boldsymbol {y}}_{1}}$. Después, calculamos un nuevo vector tangente ${\textstyle {\boldsymbol {f}}_{1}}$, nos movemos a lo largo de éste por un paso de tiempo ${\textstyle h}$ para encontrar ${\textstyle {\boldsymbol {y}}_{2}}$, y así sucesivamente.

Una variación del método de Euler puede ser obtenido de la siguiente manera. Dado que ${\textstyle {\boldsymbol {y}}={\boldsymbol {\Phi }}(t)}$ es una solución del problema de valor inicial (60), obtenemos

${\boldsymbol {\Phi }}'(t)={\boldsymbol {f}}(t,{\boldsymbol {y}}(t)).$

(62)

Aproximando la derivada en (62) por el cociente de diferencias hacia atrás ${\textstyle {\dfrac {{\boldsymbol {\phi }}(t_{n})-{\boldsymbol {\phi }}(t_{n-1})}{h}},}$ obteniendo

${\boldsymbol {\phi }}(t_{n})-{\boldsymbol {\phi }}(t_{n-1})\approx h{\boldsymbol {f}}(t_{n},{\boldsymbol {y}}_{n})),$

o bien

${\boldsymbol {y}}_{n}={\boldsymbol {y}}_{n-1}+h{\boldsymbol {f}}(t_{n},{\boldsymbol {y}}_{n})).$

(63)

Valores iniciales ${\textstyle t_{0}}$ y ${\textstyle {\boldsymbol {y}}_{0}}$, longitud de paso ${\textstyle h}$ y número de pasos ${\textstyle m}$. Aproximación ${\textstyle {\boldsymbol {y}}}$. $n$ desde 1 hasta $m$${\textstyle t_{n}=t_{n-1}+h}$. Resolver (63).

Algoritmo. 4 El Método de Euler Implícito

Aunque la implementación del Algoritmo 4 parece sencilla, tenemos el problema de resolver el sistema de ecuaciones (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

${\boldsymbol {g}}({\boldsymbol {a}})={\boldsymbol {y}}_{n-1}+h{\boldsymbol {f}}(t_{n},{\boldsymbol {a}})-{\boldsymbol {a}},$

y trabajando con el problema de optimización

$\min _{{\boldsymbol {a}}\in \mathbb {R} ^{2}}{\frac {\left\Vert g({\boldsymbol {a}})\right\Vert ^{2}}{2}}.$

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ó ${\textstyle {\boldsymbol {y}}_{n-1}}$ como punto inicial, la matriz identidad como matriz inicial y una tolerancia de ${\textstyle 10^{-8}}$.

Consideremos el problema dado en [1], el cual involucra un sistema lineal de coeficientes constantes no homogéneo de dimensión 2.

${\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}}.$

(64)

El problema tienen la siguiente solución exacta

${\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}}.$

(65)

Se tomaron un total de 1 000 puntos de la partición del intervalo temporal ${\textstyle [0,10]}$. En la Figura 10a, se muestran la solución exacta del sistema y la obtenida usando el método BFGS en el MEI. En la Figura 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.

(a) Solución con el método BFGS.	(b) Solución con el método MCN.
Figura 10: Solución numérica encontrada usando los dos métodos.

(a) Error absoluto del método BFGS.	(b) Error absoluto del método MCN.
(c) Error relativo del método BFGS.	(d) Error relativo del método MCN.
Figura 11: Error de los métodos en cada componente de la solución.

Ahora, consideremos el problema no lineal dado en [1].

${\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,$ ${\boldsymbol {y}}_{2}'={\dfrac {1}{{\boldsymbol {y}}_{2}}}-e^{t^{2}}-2te^{-t^{2}},\quad {\boldsymbol {y}}_{2}(1.5)=e^{-9/4}.$

(66)

El problema tienen la siguiente solución exacta

${\boldsymbol {y}}_{1}(t)={\dfrac {1}{t}},$ ${\boldsymbol {y}}_{2}(t)=e^{-t^{2}}.$

(67)

Se tomaron un total de 2 000 puntos de la partición del intervalo temporal ${\textstyle [1.5,1.7]}$. En la Figura 12a, se muestran la solución exacta del sistema y la obtenida usando el método BFGS en el MEI. En la Figura 12b, se muestran la solución exacta del sistema y la obtenida usando el MCN.

(a) Solución con el método BFGS.	(b) Solución con el método MCN.
Figura 12: Solución numérica encontrada usando los dos métodos.

(a) Error absoluto del método BFGS.	(b) Error absoluto del método MCN.
(c) Error relativo del método BFGS.	(d) Error relativo del método MCN.
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

${\boldsymbol {y}}_{1}'={\dfrac {4{\boldsymbol {y}}_{1}}{t}},\quad {\boldsymbol {y}}_{1}(1)=10^{-4},$ ${\boldsymbol {y}}_{2}'={\dfrac {4{\boldsymbol {y}}_{2}}{t}},\quad {\boldsymbol {y}}_{2}(1)=10^{-4},$

(68)

con ${\textstyle a=10^{-4}}$. El problema tienen la siguiente solución exacta

${\boldsymbol {y}}_{1}(t)=10^{-4}t^{4},$ ${\boldsymbol {y}}_{2}(t)=10^{-4}t^{4}.$

(69)

Se tomaron un total de 1 000 puntos de la partición del intervalo temporal ${\textstyle [1,2]}$. En la Figura 14a, se muestran la solución exacta del sistema y la obtenida usando el método BFGS en el MEI. En la Figura 14b, se muestran la solución exacta del sistema y la obtenida usando el MCN.

(a) Solución con el método BFGS.	(b) Solución con el método MCN.
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 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.

(a) Error absoluto del método BFGS.	(b) Error absoluto del método MCN.
(c) Error relativo del método BFGS.	(d) Error relativo del método MCN.
Figura 15: Error de los métodos en cada componente de la solución.

Ahora, consideremos el ejemplo dado en [10] dado por un modelo de Lotka-Volterra

${\boldsymbol {y}}_{1}'={\boldsymbol {y}}_{1}-0.5{\boldsymbol {y}}_{1}{\boldsymbol {y}}_{2},\quad {\boldsymbol {y}}_{1}(0)=2.5,$ ${\boldsymbol {y}}_{2}'=-0.75{\boldsymbol {y}}_{2}+0.25{\boldsymbol {y}}_{1}{\boldsymbol {y}}_{2},\quad {\boldsymbol {y}}_{2}(0)=1,$

(70)

donde ${\textstyle {\boldsymbol {y}}_{1}}$ y ${\textstyle {\boldsymbol {y}}_{2}}$ denotan la población de presa y depredador, respectivamente. Se tomaron un total de 7 000 puntos de la partición del intervalo temporal ${\textstyle [0,30]}$. En la Figura 16a, se muestra la solución obtenida usando el método BFGS en el MEI. En la Figura 16b, se muestra la solución obtenida usando el MCN.

(a) Solución con el método BFGS.	(b) Solución con el método MCN.
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

[1] Lambert, J. D. (1991) "Numerical Methods for Ordinary Differential Systems: The Initial Value Problem". John Wiley & Sons, Inc.

[2] Davidon, W C. (1959) "VARIABLE METRIC METHOD FOR MINIMIZATION", Volume. Technical Report ANL–5990 (revised), Argonne National Laboratory

[3] Fletcher, R. and Powell, M. J. D. (1963) "A Rapidly Convergent Descent Method for Minimization", Volume 6. The Computer Journal 2 163-168

[4] Nocedal, Jorge and Wright, Stephen J. (2006) "Numerical optimization". Springer, 2nd Edition

[5] M. Bazaraa and H. Sherali and C. M. Shetty. (2006) "Unconstrained Optimization". Nonlinear Programming. John Wiley & Sons, Ltd 343-467

[6] Neculai Andrei. (2020) "Nonlinear Conjugate Gradient Methods for Unconstrained Optimization", Volume. Springer International Publishing, 1st Edition

[7] Martínez, José Mario. (1992) "Fixed-Point Quasi-Newton Methods", Volume 29. SIAM Journal on Numerical Analysis 5 1413-1434

[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–987

[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

[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

[11] Fletcher, R. (2000) "Newton-Like Methods". Practical Methods of Optimization. John Wiley & Sons, Ltd 44-79