En Matlab®, tenemos disponible como rutina de propósito general la rutina fzero. Esta rutina es híbrida, ya que combina varias técnicas en forma heurística para lograr convergencia. La rutina decide si realiza una bisección, el método de la secante o una interpolación cuadrática inversa, dependiendo de las condiciones locales. Se recomienda practicar su uso empleando la rutina fzerogui desarrollada por Moler [1]. Esta función está limitada solamente al caso de funciones continuas y no es capaz de hallar raíces donde la función no cruza el eje cambiando de signo, como es el caso de x2 en x = 0.
Para el caso particular de los polinomios, Matlab® tiene la función roots que calcula todas las raíces. En caso de que se quiera hacer cero una ecuación cuando hay datos reales, con errores de medición, se puede plantear el problema original 2.1 como la minimización de (f)2 y aplicar entonces la función de minimización escalar fminbnd.
2.4 Sistemas de ecuaciones: el método de Newton y sus variantes
La naturaleza vectorial de un sistema de ecuaciones, en contraste con el carácter escalar de las ecuaciones en una incógnita vistas anteriormente, no modifica mucho la manera en que se resuelve iterativamente el problema. No obstante, veremos algunos detalles que se deben considerar en el momento de implementar los distintos algoritmos.
Ya hemos visto que el método de Newton corresponde a una interpolación lineal local (porque la serie de Taylor def(x) era local). De modo análogo procedemos para los sistemas de ecuaciones: aproximamos
Donde el gradiente, o jacobiano J, de f viene dado por
Por lo tanto, para obtener el nuevo valor x(k+1) tenemos que resolver el sistema lineal
Para la corrección d(k) = x(k+1) -x(k). Este es el método de Newton para n variables. Nótese que es necesario que la inversa del jacobiano exista (no debe haber singularidad). La siguiente tabla muestra el algoritmo estándar del método en n variables.
TABLA 2.2. Algoritmo típico del método de Newton
2.4.1 Variaciones del método de Newton
En la ecuación 2.15, J(k) tiene que ser evaluado en cada iteración a partir de expresiones para las derivadas parciales. Para sistemas con muchas ecuaciones, esto no es práctico. Existen varias alternativas para reducir la cantidad de trabajo o simplificar la técnica de solución.
a) Podemos remplazar J(k) por una matriz A constante, con lo que se obtiene el método de la cuerda paralela:
b) Si escogemos A = I, tenemos una iteración de punto fijo en cada variable por separado, con pérdida de rapidez de convergencia. Una elección más apropiada es A = J(0), con lo cual la iteración queda:
Esta se utiliza frecuentemente en resolución de ecuaciones diferenciales ordinarias por métodos implícitos (capítulo 3). Generalmente, el jacobiano J se revalúa cada m iteraciones, para mejorar la rapidez de convergencia.
c) Cuando las derivadas parciales no se entregan explícitamente, el jacobiano se puede determinar por diferencias finitas. Esto conduce al método de Newton discreto:
Comentarios:
De las aproximaciones recién comentadas se desprenden dos desventajas:
i) Hay que evaluar el jacobiano (a veces por diferencias finitas). Esto implica n2 operaciones por paso (se puede reducir si se actualiza el jacobiano cada m iteraciones).
ii) La solución para d(k) implica n3/3 operaciones por iteración, aunque si el jacobiano se fija, esto significa n2 operaciones/iteración.
En general, las funciones f son caras de evaluar en problemas de ingeniería de procesos, por lo que los métodos de Newton son costosos en términos de tiempo computacional. Para disminuir el número de operaciones/iteración, se recurre a los métodos de cuasi-Newton, que usan una aproximación al jacobiano, los que son actualizados en cada paso a fin de reducir el número de operaciones. Para más detalles acerca de estos métodos, se recomienda consultar el artículo de Sargent [2] y el de Dennis y Mori [3].
El siguiente ejemplo muestra un código Matlab® para resolver un sistema de n ecuaciones usando el método de Newton discreto 2.18. Se calcula cada columna ‘j’ del jacobiano introduciendo una perturbación a la variable ‘j’ mediante diferencias finitas.
Ejemplo 2.6. Aplicación del método de Newton discreto en N variables
Se desea resolver el sistema de ecuaciones no lineales:
Se usa el método de Newton discreto. El método programado debe controlar que el número de iteraciones no exceda un valor máximo, y como criterio de error utilizar la norma del vector f: err = ||f(x(k))||= norm(f) que debe ser menor a 1·10-5. ¿Cuántas raíces distintas encuentra usted?
La siguiente tabla muestra los resultados de aplicar el método de Newton discreto al sistema de ecuaciones anterior, con distintos vectores iniciales x(1). En este caso encontramos tres raíces diferentes. Esto se debe a que si eliminamos x2 y x3 en las ecuaciones del sistema, hallamos una ecuación cúbica para x1, la cual posee tres raíces reales distintas. La función que implementa el método de Newton y el sistema de ecuaciones se entregan a continuación de la tabla; hay que hacer notar que el vector f(x) debe ser un vector columna.
TABLA 2.3. Solución del ejemplo 2.6
A continuación el listado del macro de Matlab® que resuelve el ejemplo 2.6.
function [raiz,f,err,nit]=newton_dis_N(fun,x1,tol,maxit)
% Esta rutina resuelve un sistema de ecuaciones no lineales usando el método clásico de
% Newton-Raphson con jacobiano discreto. Variables de entrada:
% fun: función externa que calcula el vector f(x)
% x0: vector inicial estimado para la solución
% tol: tolerancia en error relativo de x
% maxit: máximo número de iteraciones aceptable
% variables de salida:
% raíz: vector solución estimada
% f: valor de la función en la solución estimada
% err: estimación del error alcanzado
% nit: número de iteraciones realizadas por el método
if nargin<4,
maxit=100;
end