ALGORITMO PARA LA SOLUCIÓN NUMÉRICA DE SISTEMAS DE ECUACIONES NO LINEALES MEDIANTE UNA ESTRATEGIA DE OPTIMIZACIÓN GLOBAL BASADA EN ANÁLISIS DE INTERVALOS

AN ALGORITHM FOR NUMERICAL SOLUTION OF NONLINEAR EQUATIONS SYSTEMS USING A STRATEGY OF GLOBAL OPTIMIZATION BASED ON INTERVAL ANALYSIS

ALGORÍTMO PARA A SOLUÇÃO NUMÉRICA DE SISTEMAS DE EQUAÇÕES NÃO LINEARES MEDIANTE UMA ESTRATÉGIA DE OTIMIZAÇÃO GLOBAL BASEADA EM ANÁLISE DE INTERVALOS

Luis Antonio Gómez^*, Edilberto José Reyes^**, Carlos Rodrigo Correa^***

^*Ingeniero Electrónico, Universidad Industrial de Santander. Profesor, Escuela de Ingenierías Eléctrica, Electrónica y de Telecomunicaciones, Universidad Industrial de Santander. Bucaramanga, Colombia. luisgomezardila@gmail.com.
^**Matemático y Magíster en Matemáticas, Universidad Nacional de Colombia. Profesor Asociado, Escuela de Matemáticas, Universidad Industrial de Santander. Bucaramanga, Colombia. ereyes@uis.edu.co.
^***Ingeniero Químico, Universidad Nacional de Colombia. Ph.D, Polymer Science and Engineering, Lehigh University. Profesor Titular, Escuela de Ingenierías Eléctrica, Electrónica y de Telecomunicaciones, Universidad Industrial de Santander. Bucaramanga, Colombia. crcorrea@uis.edu.co.

Artículo recibido 15-II-2012. Aprobado 13-VIII-2012
Discusión abierta hasta junio de 2013

RESUMEN

En este artículo se presenta un algoritmo para la solución numérica de sistemas de ecuaciones no lineales. Para este propósito el sistema de ecuaciones se convierte en una función de valor real para luego ser minimizada, en el sentido global, en un dominio inicial dado (una caja en ⁿ) usando análisis de intervalos. El algoritmo diseñado tiene la capacidad de determinar la existencia o no de soluciones al sistema de ecuaciones en una caja dada. Las soluciones del sistema de ecuaciones, si existen dentro de la caja dada, son expresadas mediante encerramientos por subcajas cuyo tamaño es menor que la exactitud establecida. No hay restricción acerca de la relación entre el número de ecuaciones y el número de incógnitas del sistema. Se realiza además un análisis de la convergencia del algoritmo y se muestran los resultados de su aplicación para algunos problemas de prueba.

PALABRAS CLAVE: sistemas de ecuaciones no lineales; optimización global; minimización; análisis de intervalos; caja en ⁿ.

ABSTRACT

In this paper an algorithm for the numerical solution of nonlinear equations systems is presented. For this purpose the system of equations becomes a function of real value which will be minimized, in the global sense, in a given initial domain (a box in ⁿ) using analysis of intervals. The designed algorithm has the ability to determine the existence or not of solutions to the system of equations in a given box. The solutions of the system of equations, if they exist inside the given box, are expressed by means of enclosures by sub-boxes whose size is smaller than the established accuracy. There is not restriction about the relationship between the number of equations and the number of unknowns of the system. It is also carried out an analysis of the convergence of the algorithm, and the results of their application are shown for some test problems.

KEY WORDS: systems of nonlinear equations; global optimization; minimization; interval analysis; box in ⁿ.

RESUMO

Em este artigo apresenta-se um algoritmo para a solução numérica de sistemas de equações não lineares. Para este propósito o sistema de equações converte-se em uma função de valor real para logo ser minimizada, no sentido global, em um domínio inicial dado (uma caixa em ⁿ) usando análise de intervalos. O algoritmo desenhado tem a capacidade de determinar a existência ou não de soluções ao sistema de equações em uma caixa dada. As soluções do sistema de equações, se existirem dentro da caixa dada, são expressas mediante fechamentos por subcaixas cujo tamanho é menor que a exatidão estabelecida. Não há restrição a respeito da relação entre o número de equações e o número de incógnitas do sistema. Realiza-se ademais uma análise da convergência do algoritmo e mostram-se os resultados de sua aplicação para alguns problemas de prova.

PALAVRAS-CÓDIGO: sistemas de equações não lineares; otimização global; minimização; análise de intervalos; caixa em ⁿ.

1. INTRODUCCIÓN

En el análisis y síntesis de modelos matemáticos de procesos, por lo general, se requiere la resolución de sistemas de ecuaciones no lineales para la determinación de estados de las variables del sistema que permitan que el proceso tenga un comportamiento deseado. Se requiere, además, que la técnica usada en el proceso de resolución de sistemas de ecuaciones no lineales suministre una alta exactitud en los resultados obtenidos, teniendo en cuenta que esto también depende de la exactitud del sistema de cómputo usado. Lo anterior debido a que la no linealidad en los sistemas prácticos y, por tanto, en sus modelos, hace posible que el sistema se comporte de una manera indeseada al no poder controlar la exactitud de las variables que determinan el desempeño del proceso.

Un método para afrontar esta situación es transformar el problema, de uno en el cual se tienen que hallar todas las soluciones a un sistema de ecuaciones no lineales en un dominio dado, en otro equivalente en el cual se tienen que hallar todos los puntos de óptimo global (en este caso mínimo global) de una función objetivo sobre el dominio dado (Burden y Faires, 2002). Una vez resuelto el problema de optimización global asociado al sistema de ecuaciones no lineales es posible determinar si este último tiene o no solución; si es así, la solución al sistema de ecuaciones será la misma que la solución al problema de optimización global asociado.

Debido a que, en la práctica, la búsqueda por métodos numéricos de soluciones a sistemas de ecuaciones no lineales (así como la resolución de problemas de optimización) está restringida a una región acotada del espacio determinado por las variables presentes en el sistema (Hansen, 1979), y debido a que las expresiones que determinan el sistema de ecuaciones por resolver son, por lo general, funciones continuas (o, al menos, acotadas), es posible emplear métodos basados en análisis de intervalos para la resolución de esta clase de problemas.

El análisis de intervalos, una teoría originalmente propuesta por Ramon E. Moore como una solución al problema de la acumulación de errores por redondeo en computadores digitales, es una herramienta matemática de fácil comprensión la cual, además de tener la capacidad de suministrar una gran exactitud en los resultados obtenidos al trabajar con sus métodos, tiene la ventaja de ya estar implementada en algunos sistemas de cómputo simb ólico como MATLAB (Hansen, 2004). No obstante, un problema siempre presente en cálculos numéricos extensos es la propagación de errores debido a la representación en máquina de los números reales. El redondeo hacia fuera es una herramienta que permite controlar este problema. INTLAB, una toolbox para MATLAB desarrollada por Siegfried M. Rump, tiene la capacidad de implementar aritmética de intervalos usando redondeo hacia fuera e implementar funciones elementales cuyos argumentos son intervalos (Moore, Kearfott y Cloud, 2009).

2. EL PROBLEMA

Sea el conjunto de los números reales y sea un subconjunto no vacío de ⁿ. Considere el siguiente sistema de ecuaciones:

donde, para i=1,..., m, f_i es una función cuyo dominio contiene a , y con recorrido en los números reales.

El problema entonces se puede establecer como sigue: diseñar e implementar un algoritmo que permita hallar todas las soluciones del sistema (*) en el conjunto . Es decir, hallar todas las a tales que f_i (a)=0 para i=1,..., m.

3. FUNDAMENTOS MATEMÁTICOS

Se presentan algunos teoremas básicos de la teoría de conjuntos y del análisis matemático clásico cuyas demostraciones pueden ser consultadas en Bartle (1976) y Moschovakis (2006). Un estudio detallado de los conjuntos totalmente ordenados y completamente ordenados puede consultarse en Moschovakis (2006) y Ó Searcóid (2007).

3.1 Imagen directa y minimización

Definición 1. Sean X, Y conjuntos, A X y f: X Y una función. Se define la imagen directa de A bajo f como el conjunto f [A]:={f(x): x A}. Note que f [A] Y para todo A X.

Teorema 1. Sean X, Y conjuntos y f:XY una función. Supóngase que {X_α}_αεJ es una colección de subconjuntos de X. Entonces,

Definición 2. Sean X, Y conjuntos, a X y f:XY una función. Si Y es un conjunto totalmente ordenado, entonces se dice que f obtiene su mínimo sobre X en a si y solo si f(a) ≤ f(x) para todo x X.

3.2 Sistemas de ecuaciones y optimización

Como se presenta a continuación, es posible demostrar que la solución de un sistema de ecuaciones se puede convertir en un problema de optimización. Considere el problema planteado en el numeral 2.

Sea f: la función definida por:

Note que f está bien definida, y además que la imagen de f sobre consiste solo en números reales no negativos.

Luego, dado que f(x) ≥ 0 para todo x y que es completamente ordenado (Bartle, 1976), existe el ínfimo de f[] sobre y además este ínfimo es no negativo. Por tanto, si existe el mínimo de f sobre , este mínimo debe ser no negativo. Más aun, si el sistema (*) tiene solución en la caja dada , el mínimo de f sobre existe y es cero; esto se muestra en el siguiente resultado:

Teorema 2. Suponga que el sistema (*) tiene solución en , y sea a . Entonces:

a es una solución para el sistema (*) si y solo si a minimiza f.

Demostración. Si satisface el sistema (*) entonces f_i(a)=0 para cada i=1,..., m. Luego, f(a)=0 y, como f(x) ≥ 0 para todo x , entonces a es un punto de mínimo para f.

Ahora, si minimiza a f pero no satisface el sistema (*) entonces f(a) debe ser positivo, ya que f(x)≥ 0 para todo x . Como el sistema tiene solución en , existe x* tal que f(x*)=0 y x* ≠ a. Luego, f(x*) < f(a) lo cual contradice que a sea un punto de mínimo para f .

Note la importancia de la condición general sobre la consistencia del sistema (*) en , ya que dado un sistema de ecuaciones siempre es posible construir f, y si a la minimiza esto no implica en general que el sistema tenga solución.

Por tanto, podemos convertir el problema de hallar soluciones a un sistema de ecuaciones (no lineales) en un conjunto dado , en un problema de optimización (para este caso minimización), para una función f (construida como se muestra más arriba) en el conjunto . Por tanto un algoritmo básico, fundamentado en el anterior teorema, para la búsqueda de las soluciones del sistema (*) en el conjunto , es el siguiente:

Entrada: El sistema de ecuaciones dado (*) y el conjunto .
Paso 1: Construir f.
Paso 2: Minimizar f sobre .
Paso 3: Sea a un punto de mínimo para f. Si f(a)= 0 entonces a satisface (*); en caso contrario, el sistema (*) no tiene solución en .

Es necesario, al aplicar este algoritmo, asegurar la existencia del mínimo de f requerida en el paso 2; para este fin, más adelante, se darán condiciones suficientes sobre f, y por tanto sobre las funciones f_i, para que este mínimo se obtenga.

3.3 Sucesiones

Definición 3. Sea X un conjunto no vacío. Una sucesión en X es una función cuyo dominio es el conjunto de los números naturales y con recorrido en el conjunto X.

Por lo general, si : X es una sucesión, entonces se acostumbra identificar a con la lista de sus imágenes. Por tanto, la sucesión se representa con el símbolo (_n)_nεN.

Definición 4. Una sucesión (X_i)_iεN de subconjuntos en ⁿ es llamada sucesión anidada si solo si X_i+1

Definición 5. Sea (_n)_nεN una sucesión de números reales. Se dice que (_n)_nεN es una sucesión convergente si existe s tal que, dado cualquier δ > 0, exista n* tal que |_n - s|<δ para todo n ≥ n*. El número real s se denomina límite de la sucesión.

Sea (_n)_nεN una sucesión de números reales convergente. Si s es el límite de esta sucesión, entonces se acostumbra a usar la notación (_n)_nεNs para expresar que la sucesión converge a dicho número.

Uno de los principales resultados acerca de sucesiones de números reales es que si (_n)_nεN es una sucesión convergente, entonces el límite de la sucesión es único.

3.4 Intervalos y cajas en ⁿ

Definición 6. Sean a, b números reales tales que a ≤ b. El intervalo acotado cerrado en con extremos a y b, el cual se representa como [a, b], se define como sigue:

Es inmediato que todo intervalo acotado cerrado en es un subconjunto no vacío, ya que al menos a y b pertenecen al intervalo. La colección de todos los subconjuntos de de esta forma se denota con I(); esto es:

I():={X : X es un intervalo acotado cerrado no vacío}.

Luego, si X I() existen a, b con a ≤ b, tales que X = [a, b]. En general, por brevedad en la notación, el extremo inferior a de X se denota como y el extremo superior b como ; así, si X I(), entonces X = [, ]. Si X I() es tal que = entonces se dice que el intervalo X es un intervalo degenerado o un intervalo punto (Munack, 1992). Luego, si X es un intervalo degenerado, existe x tal que X =[x, x]. Por tanto, se puede identificar el conjunto de los números reales con el subconjunto de I() consistente en todos los intervalos degenerados; esto es x [x, x] para todo x (Moore, Kearfott y Cloud, 2009). En general, si D es no vacío, entonces I(D) se define como la colección de todos los subconjuntos de D, los cuales son intervalos acotados cerrados no vacíos.

Definición 7. Sea (X_n)_nεN una sucesión en I() y sea a . Se dice que la sucesión (X_n)_nεN converge a a si y solo si las correspondientes sucesiones de extremos inferiores y extremos superiores de los intervalos convergen a a. Esto es,

(X_n)_nεN a si y solo si (_n)_nεN a y (_n)_nεN a.

Definición 8. Sean X, Y I(). Se dice que X ≤ Y si y solo si X ≤ Y.

Definición 9. Sea X I(). El punto medio y el ancho de X, que se denotan respectivamente por m(X) y ω(X), se definen por:

Definición 10. Si un subconjunto de ⁿ es un paralelepípedo rectangular cerrado acotado con lados paralelos a los ejes coordenados, decimos que es una caja en ⁿ.

Note que ⁿ es una caja si solo si = ⁿ_i=1, donde X_i I() para i=1,..., n. Es decir, una caja de ⁿ es un producto cartesiano de n elementos de I(). Por tanto una caja en ⁿ puede ser representada por una n-upla ordenada de intervalos en I(); específicamente, la caja = ⁿ_i=1 X_i puede ser representada por (X₁,..., Xn). El ancho de la caja se define como ω():=max{ω(X_i):i=1,..., n}, y su punto medio está definido como m():=(m(X₁),..., m(X_n)).

Se denota por I(ⁿ) al conjunto de todas las cajas en ⁿ y, en general, si D ⁿ es un subconjunto no vacío, I(D) denota el conjunto de todas las cajas contenidas en D (Ratschek, 1985).

Teorema 3. Sea (X_i)_iεN una sucesión anidada de cajas en ⁿ, entonces existe a ⁿ tal que a _i para todo i . Incluso más, si el ancho de la sucesión de cajas anidadas converge a cero, esto es (ω(X_i)_iεN 0, entonces tal elemento es único.

3.5 Conjuntos compactos y funciones continuas

Definición 11. Un conjunto X ⁿ se llama compacto si y solo si es un conjunto acotado cerrado. En particular, toda caja en ⁿ es un conjunto compacto.

Teorema 4. Sea {X_α}_αεJ una colección de subconjuntos compactos de ⁿ, entonces:

1. La intersección _αεJ X_α es un subconjunto compacto de ⁿ.
2. Si el conjunto de J índices es finito, entonces la unión _αεJ X_α es un subconjunto compacto de ⁿ.

Teorema 5. Sea (X_i)_iεN una sucesión anidada de subconjuntos compactos no vacíos de ⁿ. Entonces la intersección _iεN X_i es un subconjunto compacto no vacío de ⁿ.

Definición 12. Sean X ⁿ, a X y f: X una función. Se dice que f es continua en a si y solo si para cualquier ε>0 existe δ>0 tal que para todo x X, ||x - a|| < δ implica |f(x)- f(a)|<ε.

Teorema 6. Sea X un conjunto compacto no vacío y sea f: X una función continua sobre X. Entonces, existen a, b X tales que f(a)≤ f(x)≤ f(b) para todo x X. Más aun, la imagen de f sobre X es un intervalo acotado cerrado en ; esto es f[X] I().

Luego, para la función f definida en (**), si es compacto y f es continua sobre , entonces existe el mínimo de f sobre . Además para que f sea continua sobre es suficiente que cada f_i sea también continua sobre .

Por tanto, para el desarrollo y subsecuente análisis de la convergencia del algoritmo, supondremos que el conjunto es una caja en ⁿ y que cada f_i es una función continua sobre . Con esto se asegura la existencia del mínimo para f requerido en el paso 2 del algoritmo general.

3.6 Aritmética y análisis de intervalos

Se define a continuación la aritmética de subconjuntos de en términos de la aritmética de números reales.

Definición 13. Sean A, B subconjuntos no vacíos, y supóngase que representa cualquiera de las operaciones aritméticas usuales de suma, resta, multiplicación o división (+, -, ', /). Entonces:

En el caso de la división, para que la operaci ón esté bien definida, se debe cumplir que 0 B. La suma y la multiplicación así definidas cumplen las propiedades asociativa y conmutativa, y además existen identidades para estas dos operaciones. Los siguientes dos teoremas se deducen en forma inmediata de la definición de la aritmética de subconjuntos de .

Teorema 7 (Principio de inclusión). Sean A, B, C, D , subconjuntos no vacíos, tales que C A y D B. Entonces C D A B; donde es cualquiera de la operaciones de suma, resta, multiplicación o división (en el caso de la división 0 D y 0 B).

Demostración. Sea u C D, entonces existen a C y b D tales que u=a b. Como C A y D B, entonces a A y b B; luego u A B .

En lo que sigue si A, B son subconjuntos no vacíos de , entonces el producto A B será representado por AB.

Teorema 8 (Ley subdistributiva). Sean A, B, C subconjuntos no vacíos, entonces A(B+C)AB+AC.

Demostración. Sea u A(B+C), entonces existen a A b B y c C, y tales que u = a(b+c). Luego, dado que en la aritmética de números reales se cumple la ley distributiva del producto respecto a la suma, se tiene u = ab+ac. Por tanto, u AB+AC .

Note que en general no se cumple la igualdad A(B+C)=AB+AC para subconjuntos de . Por ejemplo, al establecer A=[1,4], B=[-1,2] y C=[1,2], y se tiene que A(B+C)=[1,4][0,4]=[0,16] pero AB+AC=[-4,8]+[1,8]=[-3,16]. Luego, se verifica para este caso que A(B+C) AB+AC pero A(B+C) ≠ AB+AC.

Si las operaciones antes definidas se restringen a I(), entonces estas operaciones son cerradas en el sentido de que si A, B I() entonces A B I(). Específicamente, si A=[a,b] y B=[c,d] son intervalos acotados cerrados, entonces

En el caso de la división se debe asegurar que 0 B. Los conjuntos S1 y S2 se definen como sigue:

Definición 14. Sean X una caja en ⁿ y f: X una función continua sobre X. Se define f: I(X)I() por (Y):=f[Y] para todo Y I(X).

Según el teorema 2, dado que f es continua, se tiene que (Y) I() para todo Y I(X); por lo tanto, esta función está bien definida. Además note que la función calcula para cada elemento en I(X) su imagen directa bajo la función f. Luego, dado que X es una caja en ⁿ se tiene que X I(X) y, por ello,(X) contiene el valor mínimo de la función f sobre X. Más aun f(x) ([x,x]) para todo x X.

Definición 15. Sean X una caja en ⁿ y f: X una función continua sobre X. Se dice que una función F:I(X)I() es una función de inclusión para f si y solo si para todo Y I(X) se tiene que (Y) F(Y).

Definición 16. Sean X una caja en ⁿ, f:X una función continua sobre X y F:I(X)I() una función de inclusión para f. F se llama isótona si y solo si para Y, Z I(X),

Y Z implica que F(Y) F(Z).

Si X es una caja en (un intervalo acotado cerrado) y g: X es una función continua para la cual el cálculo de su imagen sobre subcajas contenidas en X no representa dificultad alguna (como Sin, Exp, Log, etc.), se puede definir una función de inclusión G para g a través de ; es decir, G(Y):= (Y) para todo Y I(X) (Ratschek y Voller, 1991). Más aun, G definida de esta forma es isótona. Por lo general, las funciones de este tipo están predeclaradas en los sistemas de cómputo simbólico, y nos referiremos a ellas con este término. Las funciones de inclusión para funciones predeclaradas, definidas como se muestra en este párrafo, se denominan funciones de inclusión predeclaradas.

Definición 17 (Extensión natural de intervalo). Sea X una caja en ⁿ, y sea f: X una función continua sobre X tal que la expresión que la determina no contiene conectivos lógicos. Supóngase que F:I(X) I() es una función de inclusión para f. Si F se obtiene al reemplazar, en la expresión que determina la función f, cada presencia de la variable x (la cual toma valores en X) por la variable Y (la cual toma valores en I(X)) y cada presencia de una función predeclarada g en f por su función de inclusión predeclarada G, y si las operaciones aritméticas usuales en f son reemplazadas por las correspondientes operaciones de la aritmética de intervalos, entonces F se llama una extensión natural de intervalo de f.

Las funciones de inclusión tipo extensión natural de intervalo, además de ser isótonas, son de sumo interés, ya que f([x,x]) = F([x,x]) para todo x X; esto es f(x) F([x,x]) para todo x X. Incluso más, se tiene el siguiente resultado:

Teorema 9. Sea X una caja en ⁿ, y sea f: X una función continua sobre X. Supóngase que F:I(X)I() es una extensión de intervalo natural para f. Entonces, para todo Y I(X), ω(F(Y))0 si ω(Y)0.

Algunos paquetes de cómputo simbólico que implementan aritmética de intervalos, como INTLAB, automáticamente crean extensiones naturales de intervalo para funciones continuas en varias variables, lo cual facilita la construcción de funciones de inclusión isótonas para una gran variedad de funciones.

4. ALGORITMO

Dada una caja en ⁿ y el sistema de ecuaciones no lineales (*), se construye para cada f_i su respectiva función de inclusión F_i (una extensión natural de intervalo). La función de inclusión para la función f definida por (**) queda determinada, para todo Y I(), por

Note que de esta forma, F(Y) ≥ 0 para todo Y I().

El siguiente algoritmo permite desarrollar los pasos 2 y 3 del algoritmo básico propuesto en el numeral 3.2.

Entradas: Número n de dimensiones para trabajar. Caja inicial . Exactitud deseada en la determinación de las soluciones (δ). Margen de error (ε) para el valor mínimo.

1. Establecer Y=.
2. Calcular F(Y). Establecer [u, v]:=F(Y).
3. Si u>0 o v<0, entonces ir al paso 18.
4. Si ω([u, v])<ε y ω(Y)<δ, entonces la caja inicial es un encerramiento de una posible solución para el sistema (*). Fin.
5. Hacer (Y₁,[u₁, v₁], w₁, r₁)=(Y, [u, v], ω(Y), ω([u, v])). Se inicializa la lista L={(Y₁, [u₁, v₁], w₁, r₁)}.
6. Sea (Y, [u, v], w, r) el primer elemento de la lista L; donde Y=ⁿ _i=1 I_i, y cada I_i I().
7. Sea k:=min{i: w=ω(I_i)}.
8. Se biseca Y en la dirección k y se obtienen las subcajas V₁ y V₂.
9. Para i {1,2}, se calcula [u_i, v_i ]=F(V_i), w_i=ω(V_i) y r_i=ω([u_i, v_i]).
10. Si para algún i, u_i>0 o v_i<0, entonces se elimina V_i. En caso contrario, se agrega (V_i, [u_i, v_i], w_i, r_i) al final de la lista L, teniendo en cuenta que los wi sean decrecientes.
11. Se elimina (Y, [u, v], w, r) de la lista L.
12. Si L=Ø (el conjunto vacío), entonces ir al paso 18.
13. Si L=Ø se renumeran los elementos, teniendo en cuenta el orden en que se agregaron (esto asegura que los wi sean decrecientes).
14. Si w₁ ≥ δ, entonces ir al paso 6.
15. Si para todo j, r_j < ε entonces cada subcaja Y_j en la lista L es un encerramiento de una posible solución para el sistema (*). Fin.
16. Si para algún j, r_j ≥ ε entonces se renumera la lista L en orden decreciente de los r_i.
17. Se procesa la lista como sigue:
    1. Sea (Y, [u, v], w, r) el primer elemento de la lista L; donde Y=ⁿ _i=1 I_i y cada I_i I().
    2. Sea k:=min{i: w=ω(I_i)}.
    3. Se biseca Y en la dirección k y se obtienen las subcajas V₁ y V₂.
    4. Para i {1,2}, se calcula [u_i, v_i]=F(V_i), w_i=ω(V_i) y r_i=ω([u_i, v_i]).
    5. Si para algún i, u_i>0 o v_i<0, entonces se elimina V_i. En caso contrario, se agrega (V_i, [u_i, v_i], w_i, r_i) a la lista L.
    6. Se elimina (Y, [u, v], w, r) de la lista L.
    7. Si L=Ø (el conjunto vacío), entonces ir al paso 18.
    8. Si L≠Ø se renumeran los elementos, teniendo en cuenta que los ri queden en orden decreciente. Ir al paso 15.
18. El sistema (*) no tiene solución en . Fin.
19. Cada subcaja Y_i de la lista L es un encerramiento de una posible solución para el sistema (*). Fin.

5. ANÁLISIS DEL ALGORITMO

5.1 Funcionamiento del algoritmo

El funcionamiento del algoritmo se basa, en lo fundamental, en suprimir de modo iterativo subcajas donde es seguro que no se encuentran soluciones al sistema de ecuaciones (*), debido a que para todo elemento en esas subcajas se tiene que la función objetivo f definida en (**) tiene imagen estrictamente positiva.

En cada iteración del algoritmo se evalúa la función de inclusión F en la primera subcaja Y de mayor anchura presente en la lista L. Si la imagen de F sobre esta subcaja consta solo de valores positivos, la subcaja es eliminada de la lista; de lo contrario, se procede a bisecar Y a lo largo del lado de mayor longitud de la subcaja, para crear así las nuevas subcajas V₁ y V₂; se elimina Y de la lista. Si F(V_i) consta solo de valores positivos para algún i {1, 2}, se elimina V_i; de lo contrario, se agrega a la lista. Si en alguna iteración la lista L está vacía, entonces el sistema no tiene solución en la caja dada; de lo contrario, se repite el proceso.

Si en alguna iteración todas las subcajas Y presentes en la lista L satisfacen ω(Y)<δ, entonces se procede a determinar ω(F(Y)) para cada una de ellas. Si se satisface ω(F(Y))<ε para todas las subcajas de la lista, entonces el algoritmo termina y cada una de estas subcajas es un encerramiento de una posible solución para el sistema de ecuaciones. En caso contrario, si para alguna subcaja Y se tiene que ω(F(Y))≥ ε, entonces esta subcaja se procesa de forma análoga a como se explicó en el párrafo anterior. El proceso continúa hasta que la lista quede vacía, en cuyo caso el sistema no tendría solución, o hasta que todo elemento de la lista satisfaga ω(F(Y))<ε, en cuyo caso cada subcaja de la lista sería un posible encerramiento de una solución para el sistema de ecuaciones.

Por tanto, si después de aplicar el algoritmo a un sistema de ecuaciones, la lista L es no vacía, entonces cada subcaja de la lista es un encerramiento para una posible solución del sistema, cuyo ancho es menor que la exactitud δ establecida. Además, para todo elemento x de cada subcaja, se cumple |f(x)|< ε y, por tanto, |f_i(x)|< para cada i=1,..., m.

Luego, si el sistema de ecuaciones (*) tiene solución en , es posible que, además de las subcajas que encierran dichas soluciones, existan otras subcajas en la lista L que no contienen soluciones al sistema (*), pero sobre las cuales la función f determinada por (*) satisfaga |f(x)|<ε para todo elemento x de estas subcajas; esto se debe a la propiedad de continuidad que adquiere la función f al ser construida a partir de las funciones f_i, que son continuas. Estos puntos son llamados seudosoluciones generadas por el algoritmo para el sistema (*). Por lo tanto, el conjunto obtenido al aplicar este algoritmo es un superconjunto del conjunto de soluciones del sistema de ecuaciones (*) en .

Note además que no se impone restricción acerca de la relación entre el número de ecuaciones y la cantidad de incógnitas del sistema, lo que permite considerar problemas más generales. Tampoco se impone restricción alguna acerca del posible número de soluciones del sistema de ecuaciones (*), ya que el algoritmo puede determinar todas las subcajas necesarias para encerrar dichas soluciones incluso si la cantidad es infinita, esto debido a que todo subconjunto acotado de ⁿ puede ser encerrado por una unión de cajas.

5.2 Convergencia del algoritmo

Considere de nuevo el problema planteado en el numeral 2. Sea * el conjunto solución al sistema (*) en , y sea L_n la colección de subcajas generadas en la n-ésima iteración del algoritmo. Note que L_n consta a lo sumo de n elementos.

Supóngase, para el análisis de la convergencia del algoritmo, que se establece como cero la exactitud deseada en la determinación de las soluciones, esto es δ=0, y que durante el proceso no se satisfacen los criterios de terminación del algoritmo. Luego, L_n≠Ø para todo n y, además, esto implica que 0 F(). Más aun, para todo n y para todo Y L_n se tiene que 0 F(Y) (p1), de lo contrario si para algún n>1 existiera un Y en L_n tal que F(Y)>0, entonces esta subcaja habría sido eliminada en la iteración n-1.

El proceso de bisección, tal como se define en el algoritmo, implica que la sucesión de los máximos de las anchuras de las subcajas presentes en la listas converge a cero; esto es,

Además, como f es continua y F es isótona, se tiene,

Ahora, para cada n se define _n:= _YεLnY. Es inmediato que _n≠Ø para todo n , debido a que la lista L_n es no vacía. Además, como cada elemento de lista L_n es un subconjunto compacto, entonces cada _n también lo es. Además, debido al proceso de bisección, la sucesión (_n)_nεN es una sucesión anidada de subconjuntos compactos de ⁿ.

Se define := _{nεN n}. Note que es un subconjunto no vacío y compacto en ⁿ.

Proposición 1. Si L_n≠Ø para todo n , entonces * es no vacío.

Prueba. Supóngase que * es vacío. Como es compacto no vacío en ⁿ y f es continua sobre , existe a tal que f(a) ≤ f(x) para todo x X. Además se tiene que f(a)>0. Dado que se cumple (p3), existe n* tal que para todo n ≥ n*

Sea Y L_n*. Entonces se cumple que F(Y)<f(a)/10 y, dado que f[Y] F(Y), se tiene que f(x)<f(a)/10 para todo x Y, lo cual contradice que a sea un punto de mínimo para f en .

El número 10, usado en la prueba de la proposici ón anterior, se eligió teniendo en cuenta el redondeo hacia fuera.

Proposición 2. Si L_n≠Ø para todo n , entonces * .

Prueba. Sea x *. Como f(x)=0 se tiene que cualquier subcaja que contenga a x no será eliminada en el proceso. Luego, para todo n existe Y en L_n tal que x Y, y por tanto, x _n para todo n . Así, por consiguiente x .

Luego, las proposiciones 1 y 2 demuestran que el conjunto solución que se obtiene al aplicar el algoritmo converge a un superconjunto del conjunto solución del sistema de ecuaciones (*), como se había mencionado en la sección 5.1.

5.3 Consideraciones acerca de la velocidad del algoritmo

La velocidad del algoritmo está determinada por la anchura de la caja inicial ω(), el número de ecuaciones m, el número de dimensiones del dominio n (esto es, el número de incógnitas), la exactitud deseada en las soluciones δ, el margen de error establecido para el proceso de optimización ε y, especialmente, el número de soluciones al sistema dentro del dominio dado, entre otras cosas. Además, en la aritmética real puede existir más de una forma de representar la expresión que determina una función dada. Por ejemplo, las expresiones x+xx y x(1+x) determinan el mismo valor para todo x y, por tanto, representan una misma función de valor real; sin embargo, la expresión X+XX y la expresión X(1+X) determinan intervalos diferentes para algunos elementos X en I(); por ejemplo, tome X=[-1,0] (Ratschek y Voller, 1991). Más aun, X(1+X) X+XX para todo X I(). Luego, la velocidad de convergencia del algoritmo también depende de la forma en que se escriben las expresiones para cada una de las funciones que determinan el sistema de ecuaciones. En consecuencia, para aumentar la velocidad de convergencia del algoritmo es necesario reescribir la expresión que determina cada función de tal forma que, en lo posible, cada variable solo ocurra una vez.

6. IMPLEMENTACIÓN DEL ALGORITMO, PRUEBAS Y ANÁLISIS DE RESULTADOS

La implementación del algoritmo fue realizada en MATLAB, usando la toolbox INTLAB. Consta de dos programas tipo script donde el segundo se utiliza para definir las funciones que determinan el sistema.

Como primer ejemplo considere el siguiente sistema de dos ecuaciones con dos incógnitas:

Al aplicar el programa a este sistema en la caja determinada por [-1, 1]×[-1, 1] el resultado obtenido fue: "El sistema de ecuaciones no tiene solución en la caja dada". Por tanto, en esta caja no existe solución al sistema dado.

No obstante, en la caja [-1, 5]×[-1, 5] el resultado obtenido fue: "Cada una de las siguientes cajas son un encerramiento de una posible solución para el sistema de ecuaciones:

La exactitud en las soluciones es menor que: ancho = 7,152557373046875e-007 (ver tabla 1). Para cada x en las subcajas y cada función que determina el sistema de ecuaciones, abs(fi(x)) es menor que: rta = 2,582377455379701e-006".

Ello muestra que en esta nueva caja sí hay solución. Para el cómputo se usó δ=1×10^-6 y ε=1×10^-6. Luego una solución, a seis cifras decimales, para el sistema es x=0,726246, y=2,067306, y se observa que cada caja satisface el requisito establecido para la exactitud, y además para cada punto en las subcajas presentadas se tiene que, al evaluar las expresiones que determinan el sistema, la distancia a cero es menor de 2,6×10^-6, lo cual es muy inferior a =0,001, que es una cota superior para el margen de error esperado, según se explicó en la sección 5.1.

Lo anterior también demuestra la capacidad del algoritmo para determinar si dentro de una caja dada existen soluciones al sistema de ecuaciones, y en caso afirmativo las computa.

A continuación presentamos una serie de pruebas que tienen como dominio cajas cuyos lados son el intervalo [-20, 20]. Se establece δ=1×10^-8 y ε=1×10^-6 (ver tabla 2).

Considere el siguiente sistema de cuatro ecuaciones con tres incógnitas:

Al aplicar el programa a este sistema se obtuvo el siguiente resultado: "Cada una de las siguientes cajas son un encerramiento de una posible solución para el sistema de ecuaciones:

La exactitud en las soluciones es menor que ancho = 9,313225746154785e-009. Para cada x en las subcajas, y cada función que determina el sistema de ecuaciones, abs(fi(x)) es menor que: rta = 8,139288356488952e-008".

Según los resultados obtenidos, cada subcaja presentada por el programa satisface la condición de exactitud impuesta, y además para cada punto en las subcajas presentadas se tiene que, al evaluar las expresiones que determinan el sistema, la distancia a cero es menor de 8,14×10^-8. Al analizar los resultados se puede observar que una posible solución para este sistema es x=1, y=1, z=1 y de hecho lo es. Las demás subcajas que no contienen esta solución muestran que hay la posibilidad de que existan soluciones muy cercanas a la ya encontrada, no obstante, también es posible que su presencia se deba a la propiedad de continuidad de la función objetivo, como se explicó.

Considere ahora el siguiente sistema de ecuaciones:

Al aplicar el programa desarrollado se obtuvo: "Cada una de las siguientes cajas son un encerramiento de una posible solución para el sistema de ecuaciones:

La exactitud en las soluciones es menor que ancho = 9,313225746154785e-009. Para cada x en las subcajas y cada función que determina el sistema de ecuaciones, abs(fi(x)) es menor que: rta = 9,243792935043797e-008". Se observa que x=1,5352800, y=10,9461823 y x=-2,729338, y=11,3599756 son dos posibles soluciones al sistema (ver tabla 3).

Cada subcaja tiene un ancho menor que el establecido, y la distancia a cero de cada expresi ón que determina el sistema para puntos en las subcajas obtenidas es menor de 9,25×10^-8.

En adelante, para cada sistema de prueba, solo se presentarán las posibles soluciones, la exactitud con que se hallaron δ_s y la distancia a cero de las imágenes de las funciones que determinan el sistema sobre las cajas producidas ε_s.

Para el sistema

al aplicar el programa al sistema de ecuaciones y analizar sus resultados se obtuvieron las siguientes posibles soluciones: x=0,5 y=0,86602540; x=0.5, y=-0.86602540. Con δ_s menor que el establecido y εs menor de 5,72×10^-8.

Considere ahora el siguiente sistema de 5 ecuaciones con 5 incógnitas,

Se obtuvo que la solución, satisfaciendo la exactitud requerida y con un ε_s menor de 1,6×10^-6, es:

x₁=-2,56756527
x₂=2,6133588
x₃=11,1487223
x₄=0,2874417
x₅=-4,9128008

7.OBSERVACÍONES Y CONCLUSIONES

El algoritmo desarrollado, basado en la aritm ética de intervalos, por no imponer condiciones sobre la relación entre la cantidad de ecuaciones y la cantidad de incógnitas, brinda un método general para la resolución numérica de sistemas de ecuaciones no lineales, el cual proporciona una alta exactitud en las soluciones halladas, no obstante un eventual elevado tiempo de cómputo, según los criterios del problema en particular.

Algunos métodos para acelerar la búsqueda de soluciones a sistemas de ecuaciones en un dominio y que fueron probados fueron: dividir la caja inicial en subcajas para aplicar el algoritmo a cada una de ellas por separado, aplicar el algoritmo con una exactitud baja, por ejemplo , y luego aplicar otros métodos que, aunque no suministren gran exactitud al aplicarlos solos, sí proporcionan más velocidad de cómputo.

Es importante la selección de la caja inicial y de su tamaño, tanto para acelerar el cómputo como para tener seguridad de encerrar soluciones del sistema si estas existiesen.

Entre las líneas de investigación para trabajos futuros están la modificación e implementación del algoritmo para la búsqueda de soluciones en sistemas con variables complejas, y también la búsqueda de soluciones restringidas a los enteros para su aplicación a la solución de sistemas de ecuaciones diofánticos.

REFERENCIAS

Bartle, Robert G. The elements of real analysis. 2^nd ed. New York: John Wiley and Sons, 1976. 496 p.         [ Links ]

Burden, Richard L. y Faires, J. Douglas. Análisis numérico. 7^a ed. Math Learning, Thomson, 2002. 628 p.         [ Links ]

Hansen, E. R. (1979). "Global optimization using interval analysis: The one-dimensional case". Journal of Optimization Theory and Applications, vol. 29, No. 3 (November), pp. 331-344.         [ Links ]

Hansen, E. R. Global optimization using interval analysis. New York: Marcel Dekker and Sun Microsystems, 2004.         [ Links ]

Moore, Ramon E.; Kearfott, R. Baker and Cloud, Michael J. Introduction to interval analysis. Philadelphia: SIAM, 2009. 223 p.         [ Links ]

Moschovakis, Yiannis. Notes on set theory. 2^nd ed. New York: Springer, 2006.         [ Links ]

Munack, H. (1992). "On global optimization using interval arithmetic". Computing, vol. 48, No. 3-4, pp. 319-336.         [ Links ]

Ó Searcóid, Mícheál. Metric spaces. Springer Undergraduate Mathematical Series. New York: Springer, 2007. 304 p.         [ Links ]

Ratschek, H. (1985). "Inclusion functions and global optimization". Mathematical Programming, vol. 33, No. 3 (December), pp. 300-317.         [ Links ]

Ratschek, H. and Voller, R. L. (1991). "What can interval analysis do for global optimization?". Journal of Global Optimization, vol. 1, pp. 111-130.         [ Links ]