Ecuaciones No Lineales
La exposición en esta sección sigue de cerca el capítulo 3 del libro de Miranda y Fackler.
Los problemas que envuelven sistemas de ecuaciones no lineales son muy comunes en economía. Éstos se presentan en dos formatos:
- Búsqueda de raíces:
$$f(x) = 0$$
con \(f:R^n \rightarrow R^n\) y \(x\in R^n\)
- Punto fijo:
$$x = g(x)$$
con \(g:R^n \rightarrow R^n\) y \(x\in R^n\)
Ambas formas son equivalentes:
- Problema de raíces expresado como un problema de punto fijo:
$$g(x) = x - f(x)$$
- Problema de punto fijo expresado como un problema de raíces:
$$f(x) = x-g(x)$$
En muchas aplicaciones económicas, no es posible encontrar una solución analítica a un problema no lineal y por tanto la solución tiene que hallarse de manera numérica.
Problemas:
- Los métodos numéricos con iterativos y son sensibles a la condición inicial.
- Los sistemas no lineales pueden tener más de una solución.
Vamos a analizar métodos que no usan derivadas (método de bisección y método de iteración de la función) y métodos de usan derivadas (Método de Newton y Método de quasi-Newton).
Existen métodos más sofisticados basados en los métodos mencionados y sólo realizaremos una descripción escueta.
Métodos Numéricos y Ecuaciones No Lineales: Ideas Básicas
Método de Bisección
Es el método más simple y robusto para hallar la raíz de una función no lineal definida en \(R\) en el intervalo \([a,b]\).
Buscamos resolver:
$$f(x)=0$$
con \(f:R \rightarrow R\).
Es un método iterativo que evalúa el signo en los puntos extremos y en el punto medio del intervalo y en función de ello reduce el intervalo de búsqueda. Repite este proceso hasta que el intervalos resultante sea menor que un nivel de tolerancia.
Requerimientos: La solución debe estar contenida en \([a,b]\) y la función debe ser continua en el intervalo.
La siguiente función implementa el método de Bisección:
function root = bisect(f,a,b)
% bisection implementa el Método de Biseccion
% Uso: root = bisect(f,a,b)
if f(a)*f(b)>0
error('La función debe ser de signo distinto en los extremos');
else
tol = 1.5e-8;
s = sign(f(a));
x = (a+b)/2;
d = (b-a)/2;
while d>tol
d = d/2;
if s == sign(f(x))
x = x+d;
else
x = x-d;
end
end
root = x;
end
end
Para probar el método creamos un m-file con la siguiente función:
function fval = ff(x)
fval = x^3-2;
end
Buscamos la raíz en el intervalos \([1,2]\).
x = bisect(@ff,1,2);
disp(x);
1.2599
Iteración de la Función
El método de iteración de la función es un técnica iterativa usada para encontrar un punto fijo \(x=g(x)\) con \(g:R^n \rightarrow R^n\). Este método también puede ser aplicado al problema de búsqueda de raíces usando la transformación apropiada (ver más arriba).
La iteración se inicia con una conjetura inicia \(x^(0)\) para el punto fijo y se la actualiza usando la siguiente regla:
$$x^{(k+1)} \leftarrow g(x^{(k)})$$
Este proceso continúa iterando hasta que \(g(x^{(k)})\) se "suficientemente" cercano a \(x^{(k+1)}\).
La convergencia está garantizada si:
- \(g\) es diferenciable.
- la conjetura inicial \(x^{(0)}\) está suficientemente cerca del punto fijo \(x^*\).
- \(||g'(x^*)||<1\).
La siguiente figura ilustra la iteración de una función en \(R\):
La siguiente función implementa el método de iteración de la función:
function [x, gval] = fixpoint(g,x0)
% fixpoint implementa el Método de Iteración de la Función
% Uso: [x, gval] = fixpoint(f,x0)
tol = 1.5e-8;
maxiter = 1000;
x = x0;
for i=1:maxiter
gval = g(x);
if norm(gval-x)<tol
return;
end
x = gval;
end
if i == maxiter;
disp('No se obtuvo convergencia');
end
end
Para probar el método creamos un m-file con la siguiente función:
function gval = gg(x)
gval = x^0.5;
end
Buscamos el punto fijo partiendo de una conjetura de \(x=4\).
[xstar, gx] = fixpoint(@gg,4);
disp(xstar);
disp(gx);
1.0000
1.0000
El Método de Newton
En la práctica muchos problemas se resuelven usando el método de Newton o variaciones de él. De hecho las funciones de Matlab para realizar ésta tarea usan justamente variaciones del método de Newton.
El método de Newton se basa en el principio de linelización sucesiva. Así, problemas no lineales son reemplazados por una secuencia de problemas lineales cuyas soluciones convergen a la solución del problema original.
Este método puede ser aplicado tanto a problemas de búsqueda de raíz como problemas de punto fijo. En el último caso debe aplicarse la transformación apropiada (ver más arriba).
El método es iterativo y empieza con una conjetura inicial \(x^{(0)}\) para la raíz del problema. La actualización desde \(x^{(k)}\) a \(x^{(k+1)}\) ocurre resolviendo un problema de búsqueda de raíces lineal. Aplicando la aproximación de Taylor de primer orden:
$$f(x) \approx f(x^{(k)}) + f'(x^{(k)})(x-x^{(k)}) = 0$$
Por tanto:
$$ x^{(k+1)} \leftarrow x^{(k)} - [f'(x^{(k)})]^{-1} f(x^{(k)})$$
La iteración concluye cuando \(x^{(k)}\) se "suficientemente" cercano a \(x^{(k+1)}\).
La convergencia está garantizada si: - \(f\) es continuamente derivable. - La conjetura inicial \(x^{(0)}\) está suficientemente cerca de la raíz \(x^*\). - \(f'\) es invertible.
La siguiente figura ilustra la iteración de una función en \(R\):
La siguiente función implementa el método de Newton:
function x = newton(f,x0)
% newton implementa el Método de Newton
% Uso: x = newton(f,x0)
tol = 1.5e-8;
maxiter = 1000;
x = x0;
for i=1:maxiter
[fval,fjac] = f(x);
x = x - fjac\fval;
if norm(fval) < tol
break;
end
end
if i == maxiter;
disp('No se obtuvo convergencia');
end
end
Para probar el método usemos la función ya creada \(f(x) = x^3 - 2\). Buscamos la raíz partiendo de la conjetura \(x=1\).
function [fval, fjac] = ff2(x)
fval = x^3-2;
fjac = 3*x^2;
end
x = newton(@ff2,1);
disp(x);
1.2599
El Método de Quasi-Newton
El problema del método de Newton es que requiere la matriz Jacobiana para realizar la linealización. En muchos problemas es difícil encontrar dicha matriz.
El Método de Quasi-Newton sigue la misma lógica que el método de Newton pero reemplaza la matriz Jacobiana por una aproximación fácil de computar.
El costo viene asociado a que es más lento en converger y requiere de una aproximación inicial a la matriz Jacobiana.
Método de la Secante
El Método de la Secante es el método de Quasi-Newton univariado más usado. La idea es reemplazar la derivada con una aproximación construida con los valores de la iteración anterior. La regla de actualización sería entonces:
$$ x^{(k+1)} \leftarrow x^{(k)} - \frac{x^{(k)}-x^{(k-1)}}{f(x^{(k)})-f(x^{(k-1)})} f(x^{(k)})$$
La siguiente función implementa el método de la secante:
function x = secante(f,x0)
% secante implementa el Método de la Secante
% Uso: x = secante(f,x0)
tol = 1.5e-8;
maxiter = 1000;
x1 = 0.5*x0;
for i=1:maxiter
fval0 = f(x0);
fval1 = f(x1);
x = x1 - ((x1-x0)/(fval1-fval0))*fval1;
if abs(fval1) < tol
break;
end
x0 = x1;
x1 = x;
end
if i == maxiter;
disp('No se obtuvo convergencia');
end
end
Para probar el método usemos la función:
function fval = ff(x)
fval = x^3-2;
end
Buscamos la raíz conjeturando \(x=1\).
x = secante(@ff,1);
disp(x);
1.2599
Método de Broyden
El Método de Broyden es el método de Quasi-Newton multivariado más usado. El método parte con una conjetura inicial para la raíz \(x^{(0)}\) y para la matriz Jacobiana \(A^{(0)}\). Las ecuaciones de actualización son para la raíz:
$$f(x) \approx f(x^{(k)}) + A^{(k)}(x-x^{(k)}) = 0$$
$$ x^{(k+1)} \leftarrow x^{(k)} - [A^{(k)}]^{-1} f(x^{(k)})$$
y para la matriz jacobiana:
$$f(x^{(k+1)}) - f(x^{(k)}) = A^{(k)}(x^{(k+1)} - x^{(k)})$$
$$A^{(k+1)} \leftarrow A^{(k)} + \left[f(x^{(k+1)}) - f(x^{(k)}) - A^{(k)}d^{(k)}\right] \frac{d^{(k)}}{d^{(k)} d^{(k)}}$$
con \(d^{(k)} = x^{(k+1)}-x^{(k)}\).
El toolbox CompEcon de Miranda y Fackler (2002) tiene implementado éste y los otros métodos descritos antes.
Las Funciones fzero y fsolve
Matlab cuenta dos funciones diseñadas para lidiar con problemas no lineales: (1) fzero y (2) fsolve.
fzero
fzero es un función para buscar la raíz de una función definida en los reales. Usa una combinación de métodos entre los cuales se encuentra los de bisección y el de la secante.
La sintaxis de esta función es:
[x,fval,exitflag] = fzero(fun,x0,opt)
donde fun es una función definida en \(R\), x0 es la conjetura inicial (el intervalo de búsqueda si se provee un vector \(2 \times 1\)) y opt define las opciones del método.
Entre los resultados tenemos: x la raíz, fval el valor de la función evaluada en la raíz, y exitflagun indicador con el resultado de método ($x>0$ indica que se encontró la solución).
Las opciones se manejan usando la función optimset de Matlab. Sólo para tomar un ejemplo:
opt = optimset('Display','iter','TolX',num,'MaxIter',num)
No usar opt sólo hace que la función use las opciones por defecto.
[x,fval,exitflag] = fzero(@ff,1)
x =
1.2599
fval =
0
exitflag =
1
Ahora usemos las opciones de la función.
%opts = optimset('Display','iter');
%x = fzero(@ff,1,opts);
fsolve
fsolve es un función para resolver numéricamente sistemas de ecuaciones no lineales. Esta función usa variaciones (más robustas y eficientes) del método de Newton.
La sintaxis de esta función es:
[x, fval, exitflag] = fsolve(fun,x0,options)
donde fun es un vector de n función definidas en \(R^n\), \(x0 \in R^n\) es la conjetura inicial y opt define las opciones del método (al igual que en el caso anterior).
Resolvamos el siguiente sistema de ecuaciones no lineales usando fsolve.
$$e^{-e^{-(x+y)}} = y(1+x^2)$$
$$x \cos y + y \sin x = \frac{1}{2}$$
Primero requerimos generar un m-file con las ecuaciones:
function f = sistemanl(x)
f = zeros(length(x),1);
f(1) = exp(-exp(-x(1)+x(2))) - x(2)*(1+x(1)^2);
f(2) = x(1)*cos(x(2)) + x(2)*sin(x(1)) - 0.5;
end
La conjetura inicial será \(x=y=0\).
x0 = [0,0];
[x, f, ef] = fsolve(@sistemanl,x0)
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.
x =
0.3931 0.3366
f =
1.0e-10 *
-0.2027
-0.0095
ef =
1
sistemanl(x)
ans =
1.0e-10 *
-0.2027
-0.0095
Ahora usemos las opciones de la función.
% opts = optimset('Display','iter');
% x0 = [0,0];
% [x, f, flag]= fsolve(@sistemanl,x0,opts)