Academia.edu no longer supports Internet Explorer.
To browse Academia.edu and the wider internet faster and more securely, please take a few seconds to upgrade your browser.
M etodos Num ericos con Matlab
Alberto Morales VazquezRelated papers
Life Cycle and Performances Progressing Of Artillery Barrels
2011
The present researches of advancing the gun barrels and especially for artillery barrels, are dealing with materials other than steel or on combination b etween materials which are advanced and new innovative process for increasing of the performances and life cycle of all gun calibers. Considerable efforts we re made for developing new technique, in addition to the in vestigation of new materials. In this research, 3-techniques are describes for in creasing the mechanical properties of thick-walled cylinders which are subjected to high internal load s due to explosives. In autofrettage process, an el astic- plastic state produces at the cylinder wall by load ing of the barrel tube with a controlled inner pres sure. In elastic process, the tube is unloading, this case d oes not recover the initial dimensions, therefore; the residual stress and strain stares can be determined. The practical data were obtained from a standard te nsion test specimen, collected from the D-30 gun barrel ...
From the Private to the Public Space: Domestic and Urban Architecture of Islamic Granada
A Companion to Islamic Granada, 2022
M¶etodos Num¶ericos con Matlab
Mishelle Segu¶³
1
1
mishelle@gmx.de
Esta es una primera versi¶on, por lo que se agradece cualquier comentario o sugerencia.
1
Introducci¶on
El prop¶osito de estas notas es dar una introducci¶on a los sistemas de software interactivo
Matlab (Matrix laboratory) y Octave, los cuales son equivalentes, con la ventaja de que
Octave es un software libre2 . Esto con el objetivo de que posteriormente los alumnos puedan
crear sus propios programas de programaci¶on din¶amica, los cuales pueden incluir gr¶a¯cas,
regresiones, o cualquier computo num¶erico que se desee.
Dado que Matlab est¶a dise~
nado para resolver problemas num¶ericamente, es decir, que
tienen un aritm¶etica de precisi¶on ¯nita, nos va a producir soluciones aproximadas y no
exactas. Sin embargo, esto no resta credibilidad a este programa, ya que es uno de los
programas m¶as ultizados entre los economistas que realizan calibraciones, y es muy f¶acil de
usar.
2
Si tienen Windows, lo pueden bajar de http : ==sourceforge:net=project=showfiles:php?group id =
2888 y se llama \Octave 2.1.73 for Windows Installer".
2
Temario
1. Introducci¶on a Matlab: Sint¶axis, c¶alculos sencillos y gr¶a¯cas b¶asicas
2. Integraci¶on y diferenciaci¶on num¶erica
3. Resolver un sistema de ecuaciones no lineales con m¶etodos de gradiente (Newton,
secante)
4. Resolver una ecuaci¶on en diferencias con m¶etodos num¶ericos (Gauss-Seidel).
5. Programaci¶on din¶amica num¶erica:
² Iteraci¶on de la funci¶on de valor
² Programaci¶on din¶amica determin¶istica
² Programaci¶on din¶amica estoc¶astica
6. Otros m¶etodos de la aproximaci¶on: polinomios, interpolaci¶on
7. Aplicaciones econ¶omicas
3
1
Introducci¶
on a Matlab
1.1
Sint¶
axis y c¶
alculos sencillos
Antes que nada se~
nalemos que cuando una quiere iniciar Matlab, ¶este puede tardar unos
segundos en abrir (de 10 a 15 seg.).
Algunas cosas que uno debe de saber de Matlab son:
² Las min¶
usculas y las may¶
usculas no son equivalentes
² Un punto y coma al ¯nal del comando har¶a que no se vea en la pantalla el resultado
² Matlab utiliza par¶entesis () y corchetes [], los cuales no son intercambiables
² Las teclas con las °echas hacia arriba y hacia abajo ayudan para desplazarse entre los
comandos previos que se realicen
² Para seleccionar ayuda, hay que teclear help seguido del comando del cual se requiere
la ayuda
² Matlab se puede cerrar tecleando quit o exit
Explicar las barras de herramienta
Notemos que en la ventana de comandos, cuando asignamos un nombre a una expresi¶on,
es s¶olo un nombre y NO representa una variable matem¶atica, como lo hace Maple.
Veamos ahora unos ejercicios de s¶³mbolos y puntuaci¶on3 para realizarlos en la computadora:
Ver ejercicios1.1.pdf
3+7
a = sqrt(3)
b = sin(pi)
c = [1 : 10]
d=1:2:8
e = [1 4];
e0
e: ¤ e
f = [1 4; 2 5]
f ^2
3
http://www.indiana.edu/~statmath/math/matlab/gettingstarted/
4
Recordemos ahora algunas propiedades b¶asicas de las matrices. Digamos que tenemos el
vector ¯la g = [1 2 3] y el vector columna h = [4; 5; 6], entonces podemos como los vectores
son conformables, podemos multiplicar g ¤ h; o tambi¶en podemos computar el producto
interno con la funci¶on dot(g; h):
Hacerlo y comprobar que da 32
Recordar que el producto a ¤ a no est¶a de¯nido, dado que las dimensiones no son
compatibles, con lo cual si lo tecleamos, Matlab nos arroja un error.
Matlab tiene muchas funciones matem¶aticas para las matrices, por ejemplo exp(g);
log(h), sqrt(g); abs(g); sin(g);etc.
Por default Matlab nos muestra el resultado con 4 d¶³gitos decimales, aunque siempre
los guarda completos en la memoria y computa el equivalente a 14 d¶³gitos decimales. Si
deseamos ver dichos d¶³gitos, escribimos format long, por ejemplo sqrt(g). Si deseamos
regresar al formato anterior, s¶olo escribimos format.
Notemos tambi¶en que Matlab nos muestra los n¶
umeros muy grandes y muy peque~
nos
en notaci¶on de exponencial, con una factor de 10 denotado por e.Por ejemplo, 2^(¡24).
Para redondear n¶
umeros, Matlab tiene los comandos round; f ix; ceil y f loor; por
ejemplo round(g):
Matlab tambi¶en puede realizar varias funciones de an¶alisis de datos, por ejemplo sumar,
sacar la media,mediana, desviaci¶on est¶andar o ver la diagonal: sum(g); mean(g); median(g);
std(g); diag(f ):
Veamos a continuaci¶on algunas funciones de algebra lineal que Matlab realiza. Por
ejemplo, recordando que h es un vector columna, tal vez quisieramos resolver el sistema
lineal B ¤ x = h, en donde B = [¡3 0 1; 2 5 ¡7; ¡1 4 8]; a trav¶es del operador n, de
la forma x = Bnh: Recordemos que la divisi¶on inclinada a la derecha a=b, signi¯ca ab , pero
inclinada a la izquierda anb, signi¯ca ab :
Ver que el resltado es -1.3717, 1.3874, -0.1152.
Uno podr¶³a comprobar el resultado a trav¶es de la norma euclideana del residuo: norm(B¤
x ¡ h); que es practicamente cero. Los valores propios (eigenvalues) se pueden encontrar usando eig, por ejemplo j = eig(B):
Si queremos saber el tama~
no de la matriz, simplemente escribimos size(B): Si la queremos transponer ponemos B 0 :
Como ya vimos podemos generar una secuencia se n¶
umeros a trav¶es de los dos puntos,
como en el caso c; el cual tambi¶en lo podemos observar sin los corchetes. Notamos que en
5
este caso va del 1 al 10 en valores de uno, pero >c¶omo hacemos para que incrementen on
otras unidades? Esto es muy sencillo, s¶olo se incorporan los incrementos especi¯candolo
como m : s : n; en donde m es el valor con el que se desea inicializar, s es el paso con el cual
se desea incrementar, y n es el valor ¯nal. Por ejemplo, k = 1 : 3 : 10 ¶o l = 1 : 0:2 : 10:
En Matlab tambi¶en podemos hacer matrices m¶as grandes usando otras que ya ten¶³amos,
por ejemplo:
C = [B; g]:
Si deseamos s¶olo un elemento de la matriz, lo se~
nalamos como C(i; j); en donde i es
la ¯la y j es la columna, por ejemplo C(2; 3); que nos de ¡7: Pero >qu¶e tal si queremos
una submatriz? Entonces s¶olo seleccionamos C(i1 : i2; j1 : j2);por ejemplo, C(2 : 3; 1 : 2):
En Matlab tambi¶en podemos escoger s¶olo las columnas o las ¯las, por ejemplo C(:; j) es la
columna j de C, y C(i; :) es la ¯la i:
Matlab tambi¶en tiene por default algunas matrices, como la identidad y la de ceros
y unos, a trav¶es de eye, zeros y ones.
Por ejemplo, I3 = eye(3; 3), Y = zeros(3; 5);
Z = ones(2): Si lo que deseamos es una secuencia de s¶olo un n¶
umero repetido varias veces,
por ejemplo si el n¶
umero cero lo queremos repetido 3 veces, escribimos x(3) = 0:
Posteriormente veremos que es u
¶ til generar valores aleatorios a trav¶es de los comandos
rand y randn, los cuales vienen de la distribuci¶on uniforme [0,1] y de la normal (0,1),
respectivamente. Por ejemplo, F = rand(3) y G = randn(1; 5):
Matlab puede computar el in¯nito, a trav¶es del comando inf, sin embargo recordemos
que no podemos dividir inf/inf, porque matem¶aticamente no est¶a de¯nido y en la pantalla
de Matlab aparecer¶³a NaN(Not a Number).
Por u
¶ ltimo veamos que Matlbab, como muchos lenguajes de programaci¶on, trabaja con
bucles (loops). Un ejemplo es:
1+
1
1+
1
1
1+ 1+1
Dicha soluci¶on se puede obtener generando el bucle
m = 2;
for n = 1 : 3; m = 1 +
1
;
m
end
m
6
Si deseamos ver todas nuestras variables escribimos who y se desplegaran todos los
nombres de nuestras variables. Aunque si se desea mayor informaci¶on acerca de las variables,
escribimos whos.
1.2
Gr¶
a¯cas
Las gr¶a¯cas m¶as sencillas de hacer son las gr¶a¯cas de puntos en el plano cartesiano. Por
ejemplo,
x = [0; 3; 6:1; 7];
y = [3; 8; 10; 9:2];
plot(x; y)
N¶otese que por default, Matlab junta los puntos con l¶³neas rectas, pero para tener s¶olo los
puntos ponemos
plot(x; y; 0o0)
Tambi¶en podemos gra¯car una funci¶on con los puntos de x; por ejemplo,
plot(x; sin(x))
Podemos tambi¶en dar no s¶olo puntos, sino intervalos, por ejemplo,
x = (0 : :1 : 2 ¤ pi);
y = sin(x);
plot(x; y)
7
¶o
x = (¡5 : :1 : 5);
y = x:=(1 + x:^2);
plot(x; y)
Para ver las diferentes opciones que tenemos para gr¶a¯car, nos podemos ir a la ayuda,
donde ah¶³ nos se~
nalan c¶omo poner diferentes colores, marcadores y tipos de l¶³neas. Veamos
algunos ejemplo:
plot(x; y; 0g ¡ 0)
plot(x; y; 0p0;0 MarkerSize 0 ; 10)
Podemos a~
nadirle a las gr¶a¯cas t¶³tulo y nombre a los ejes, por ejemplo
plot(x; y); :::
title(0Gr¶af ica 10); :::
xlabel(0eje x0)
El comando hold on lo usamos para poner en una misma gr¶a¯ca, otra gr¶a¯ca, por
ejemplo usando el ejemplo anterior, posteriormente escribimos:
hold on
v=
8
w=
plot(v; w)
Subplot hace que aparezcan varias gr¶a¯cas en una misma ventana.
Realizemos ahora algunos ejercicios de gr¶a¯cas4
x = ¡10 : 0:5 : 10;
y = x:^2;
plot(x; y)
t = 0 : 0:1 : 2 ¤ pi;
x = cos(t);
y = sin(t);
plot(x; y)
t = 0 : pi=5 : 2 ¤ pi;
u = cos(t);
v = sin(t);
f igure
plot(u; v)
plot(x; y; 0r ¡ 0; u; v; 0b¤ : 0)
figure
subplot(1; 2; 1)
plot(x; y)
title(0F ino0)
subplot(1; 2; 2)
plot(u; v)
title(0T osco0)
Ver ejercicios1.2.pdf
4
idem
9
2
Integraci¶
on y diferenciaci¶
on num¶
erica
En este cap¶³tulo usaremos m¶etodos num¶ericos para resolver problemas de c¶alculo y ecuaciones
diferenciales. Veremos algunas t¶ecnicas num¶ericas para encontrar derivadas e integrales
de¯nidas.
Por ejemplo, c¶omo calculamos el ¶area de un triangulo usando Matlab. Lo que podemos
empezar a hacer es un archivo .m, en el cual de¯namos nuestro objetivo:
function A=Atri(b,h)
A=b*h/2;
Entonces en la ventana de Matlab s¶olo escribimos por ejemplo Atri(2,5) para saber
el ¶area de un tri¶angulo de base 2 y altura 5.
2.1
Integraci¶
on num¶
erica
El problema general en la integraci¶on num¶erica, tambi¶en llamado cuadratura, es el de comR
putar D f (x)dx en donde f : Rn ! R es una funci¶on integrable sobre el dominio D ½ Rn .
Todas las f¶ormulas de integraci¶on num¶erica usan un n¶
umero ¯nito de evaluaciones del inteR
grando, f , y usan una suma ponderada de estos valores para aproximar D f (x)dx . Existen
varios m¶etodos, aunque a veces los m¶as sencillos no son tan e¯cientes y los m¶as complejos
tardan mucho y necesitan varios requisitos. Nosotros no veremos los m¶etodos muy complejos.
Las f¶ormulas de cuadratura de las cotas de Newton eval¶
uan a f en un n¶
umero ¯nito
de puntos, usan esta informaci¶on para construir una aproximaci¶on polinomial por pedazos
R
de f , e integran esta funci¶on f para aproximar D f(x)dx.
10
f(x)
Q
U
V
S
R
P
a
b
x
Gr¶a¯ca 2.1. Reglas de las cotas de Newton
En la gr¶a¯ca 2.1 la funci¶on f pasa por los puntos P; Q y R. La integral
Rb
a
f (x)dx es
el area bajo la la funci¶on f y el eje de las abscisas. Podemos observar tres aproximaciones
directamente. La caja aUQVb aproxima a f con una funci¶on constante de f en Q, que es
el punto medio de [a; b]: El trapezoide aP Rb aproxima a f con una l¶³nea recta a trav¶es de
los puntos P y R: Por u
¶ltimo, el ¶area bajo la curva punteada P QSR aproxima a f con
una par¶abola a trav¶es de P; Q y R: Estas aproximaciones est¶an basadas en uno, dos y tres
evaluaciones de f, respectivamente, las cuales son usadas para computar polinomios de orden
uno, dos y tres. Este enfoque nos lleva a las f¶ormulas cuadr¶aticas de las cotas de Newton.
Antes de ver estas f¶ormulas de las cotas de Newton y las propiedades de los errores,
recordemos que es el teorema de Taylor.
Teorema de Taylor
Es el teorema m¶as utilizado en el an¶alisis num¶erico y veamoslo en sus versiones en R
y en Rn :
Teorema 2.1. (Teorema de Taylor en R) Si f 2 C n+1 [a; b] y x; x0 2 [a; b]; entonces
f(x) = f (x0 ) + (x ¡ x0 )f 0 (x0 ) +
+::: +
(x ¡ x0 )2 00
f (x0 )
2!
(x ¡ x0 )n (n)
f (x0 ) + Rn+1 (x)
n!
11
en donde
1
Rn+1 (x) =
n!
Z
x
x0
(x ¡ t)n f (n+1) (t)dt =
(x ¡ x0 )n+1 (n+1)
f
(»)
(n + 1)!
xi
para » 2 [x; x0 ]
El teorema de Taylor dice, esencialmente, que uno puede usar informaci¶on de las
derivadas en un s¶olo punto para construir una aproximaci¶on polinomail de una funci¶on
en un punto.
Teorema 2.2. (Teorema de Taylor en Rn ) Sea f : Rn ! R y es C k+1 ; entonces para
x0 2: Rn ,
n
n
n
X
@f 0
1 X X @ 2f
0
f(x) = f (x ) +
(x )(xi ¡ xi ) +
(x0 )(xi ¡ x0i )(xj ¡ x0j )
@x
2!
@x
@x
i
i
j
i=1
i=1 j=1
0
n
n
@kf
1 X X
:::
(x0 )(xi1 ¡ x0i1 ):::(xik ¡ x0ik ) + ©(jjx ¡ x0 jjk+1 )
+::: +
k! i =1 i =1 @xi1 :::@xik
1
2.1.1
k
Regla del punto medio
Esta es la f¶ormula de cuadratura m¶as sencilla, y viene dada por
Z
a
b
f(x)dx = (b ¡ a)f (
a+b
(b ¡ a)3 00
)+
f (»)
2
24
para alguna » 2 [a; b]: Como haremos tambi¶en m¶as adelante, los primeros t¶erminos compren-
den la regla de integraci¶on y el u
¶ltimo t¶ermino es el error de dicha regla. Esta regla es un
claro ejemplo de una regla abierta, que es una regla que no utiliza los puntos extremos.
R¼
Por ejemplo, si queremos integrar la funci¶on S = 0 sin(x)
dx con la regla del punto
x
medio, vamos a obtener
Z
¼
0
sin( ¼2 )
sin(x)
1
dx ¼ ¼ ¼ = ¼ ¼ = 2
x
2
2
Vemos en la siguiente gr¶a¯ca el verdadero valor del ¶area con aquel encontrado con esta
regla.
12
Gr¶a¯ca 2.2. Area dada por la integral S y la aproximaci¶on usando la regla del punto medio
Sin embargo, esta regla es poco precisa para aproximar el valor deseado, por lo que
es m¶as conveniente romper el intervalo [a; b] en peque~
nos intervalos, aproximar la integral
en cada uno de estos peque~
nos intervalos y sumar estas aproximaciones. El resultado es
una regla compuesta. Sea n ¸ 1 el n¶
umero de intervalos y de¯namos a h como h =
(b¡a)
n
y
xj = a + (j ¡ 21 )h; para j = 1; 2; :::; n: Entonces la regla compuesta del punto medio es:
Z
b
f (x)dx = h
a
n
X
f (xj ) +
j=1
h2 (b ¡ a) 00
f (»)
24
para alguna » 2 [a; b]: Notemos que el error es proporcional a h2 ; es decir, que si doblamos
el n¶
umero de intervalos, esto reduce a la mitad el tama~
no del paso h y por lo tanto reduce
el error un 75 porciento. Es por esto que esta regla compuesta del punto medio converge
cuadr¶aticamente para f 2 C 2 :
Dado que esta regla es muy sencilla, pocas veces se llega a utilizar en la pr¶actica en
econom¶³a, por lo que demos ¶enfasis a las siguientes dos reglas.
2.1.2
Regla trapezoidal
Se basa en la aproximaci¶on lineal de f usando s¶olo los valores de f en los puntos extremos
de [a; b]: La regla trapezoidal consiste en:
Z
a
b
f (x)dx =
(b ¡ a)
(b ¡ a)3 00
[f (a) + f (b)] ¡
f (»)
2
12
13
para alguna » 2 [a; b]: Esta regla es el ejemplo m¶as simple para la regla cerrada; la cual es
una regla que usa los puntos extremos.
2
Por ejemplo, f(x) = e¡x ; con a = 0 y b = 2: Usando la regla trapezoidal tenemos que
Z
2
0
2
2
2
e¡x dx ¼ 1[e¡0 + e¡2 ] ¼ 1: 0183
La funci¶on y la aproximaci¶on lineal se pueden ver en la siguiente gr¶a¯ca:
2
Gr¶a¯ca 2.3. f (x) = e¡x y una l¶inea recta obtenida por la regla trapezoidal
En t¶erminos de Matlab, lo escribiriamos en un archivo .m de la siguiente manera:
function q=trapecio(f,a,b)
ya=feval(f,a);
yb=feval(f,b);
q=(b-a)*(ya+yb)/2;
De¯niendo cualquier funci¶on en otro archivo .m, por ejemplo:
function y=ejemplo(x)
y=exp(-x^2);
Luego en la ventana de Matlab s¶olo escribimos trapecio('ejemplo',0,2), si es que
queremos integrar de 0 a 2.
14
Esta regla tambi¶en tiene el problema de que trata de aproximar con una recta, lo cual
suena un poco il¶ogico para aproximar funciones en C n ; en donde n ¸ 2: Por lo tanto, aqui
tambi¶en querremos utilizar la regla compuesta trapezoidal. Sea h =
(b¡a)
n
, xj = a + jh; para
j = 1; 2; :::; n; y denotemos fi a f (xj ); entonces dicha regla es:
Z
b
f (x)dx =
a
h2 (b ¡ a) 00
h
[f0 + 2f1 + ::: + 2fn¡1 + fn ] ¡
f (»)
2
12
para alguna » 2 [a; b]:
Por ejemplo, f (x) = x1 ; con a = 1, b = 2 y n = 2: Usando la regla compuesta trapezoidal
tenemos que
Z
2
1
1
1 1
1
1
17
1
dx ¼ [f (1) + 2f (1:5) + f(2)] = [ + 2
+ ]=
¼ 0:7083
x
4
4 1
1:5 2
24
La funci¶on y las aproximaci¶on lineales se pueden ver en la siguiente gr¶a¯ca:
Gr¶a¯ca 2.4. y =
1
x
y la aproximaci¶on trapezoidal en [1,1.5] y [1.5,2]
En t¶erminos de Matlab, lo escribiriamos en un archivo .m de la siguiente manera:
function I=Trap(f,a,b,n)
h=(b-a)/n
S=feval(f,a);
for i=1:n-1
x(i)=a+h*i;
15
S=S+2*feval(f,x(i));
end
S=S+feval(f,b);
I=h*S/2;
De¯niendo cualquier funci¶on en otro archivo .m, por ejemplo:
function y=ejemplo2(x)
y=1/x;
Luego en la ventana de Matlab s¶olo escribimos Trap('ejemplo2',1,2,2), si es que
queremos integrar de 1 a 2.
Checar 330.0285 M672A
Regla de Simpson
2.1.3
Una aproximaci¶on por pedazos lineal de f en la regla compuesta trapezoidal puede ser
innecesaria si f es suave. Otra manera es utilizar una aproximaci¶on por pedazos cuadr¶atica
de f , la cual utiliza los valores de f en a; b y en el punto medio 21 (a +b): La regla de Simpson
en el intervalo [a; b] es:
Z
b
f (x)dx =
a
(b ¡ a)
a+b
(b ¡ a)5 (4)
[f (a) + 4f(
) + f (b)] ¡
f (»)
6
2
2880
para alguna » 2 [a; b]:
Por ejemplo, f (x) =
Z
1
0
1
;
1+x2
con a = 0 y b = 1: Usando la regla de Simpson tenemos que
1
1
1
1
1
1
1
47
dx ¼ [f (0) + 4f( ) + f (1)] = [
+4
]=
¼ 0:7833
2 +
2
2
2
1
1+x
6
2
6 1+0
1+1
60
1+ 2
El valor exacto de la integral es arctan(1) =
¼
4
¼ 0:7853
La funci¶on y el polinomio cuadr¶atico que pasa por los puntos (0,1),(.5,.8) y (1,.5) se
pueden ver en la siguiente gr¶a¯ca:
16
Gr¶a¯ca 2.5. f (x) =
1
1+x2
y la aproximaci¶on cuadr¶atica usando la regla de Simpson
No es de sorprenderse que el valor aproximado de la integral con la regla de Simpson
sea muy bueno, porque las dos funciones son similares.
En t¶erminos de Matlab, lo escribiriamos en un archivo .m de la siguiente manera:
function s=Simpson(f,a,b)
ya=feval(f,a);
yb=feval(f,b);
c=(a+b)/2;
yc=feval(f,c);
q=(b-a)*(ya+4*yc+yb)/6;
De¯niendo cualquier funci¶on en otro archivo .m, por ejemplo:
function y=ejemplo(x)
y=exp(-x^2);
Luego en la ventana de Matlab s¶olo escribimos Simpson('ejemplo',0,2), si es que
queremos integrar de 0 a 2.
Construyamos ahora la correspondiente regla compuesta de (n + 1)-puntos sobre [a; b]:
Sea n ¸ 2 un n¶
umero par de intervalos, entonces h =
17
(b¡a)
n
, xj = a + jh; para j = 0; 1; :::; n;
y la regla compuesta de Simpson es:
Sn (f) =
h
h4 (b ¡ a) (4)
[f0 + 4f1 + 2f2 + 4f3 + ::: + 4fn¡1 + fn ] ¡
f (»)
3
180
para alguna » 2 [a; b]: Esta regla lo que hace b¶asicamente es que toma tres particiones
consecutivas de xj ; usa la funci¶on cuadr¶atica para aproximar a f y la integra para aproximar
la integral sobre el intervalo.
2
Por ejemplo, f(x) = e¡x ; con a = 0, b = 2 y n = 4: Usando la regla compuesta de
Simpson tenemos que
Z
0
2
2
1
3
1
[f(0) + 4f ( ) + 2f (1) + 4f ( ) + f (2)]
6
2
2
1 ¡02
1 2
3 2
2
2
(e
+ 4e¡( 2 ) + 2e¡1 + 4e¡( 2 ) + e¡2 ) = 0:8818
=
6
e¡x dx ¼
La funci¶on y las aproximaciones cuadr¶aticas se pueden ver en la siguiente gr¶a¯ca:
2
Gr¶a¯ca 2.6. y = e¡x y su aproximaci¶on usando la regla de Simpson con n = 4
En t¶erminos de Matlab, lo escribiriamos en un archivo .m de la siguiente manera:
function J=Simp(f,a,b,n)
h=(b-a)/n;
S=feval(f,a);
for i=1:2:n-1
x(i)=a+h*i;
18
S=S+4*feval(f,x(i));
end
for i=2:2:n-2
x(i)=a+h*i;
S=S+2*feval(f,x(i));
end
S=S+feval(f,b);
J=h*S/3
De¯niendo cualquier funci¶on en otro archivo .m, por ejemplo:
function y=ejemplo(x)
y=exp(-x^2);
Luego en la ventana de Matlab s¶olo escribimos Trap('ejemplo2',0,2,4), si es que
queremos integrar de 0 a 2.
Notemos que usando localmente una aproximaci¶on cuadr¶atica de f tenemos un error de
orden h4 ; mientras que usando localmente una aproximaci¶on lineal de la regla trapezoidal tenemos un error de orden h2 : Por lo que con cualquier aproximaci¶on lineal de funciones suaves,
aproximaciones de altos grados dar¶an como resultado errores asimpt¶oticamente menores.
Veamos un ejemplo num¶erico en la tabla 2.1, en donde vemos que en algunos casos se
necesitan varios puntos para acercarnos a los valores verdaderos. Vemos que ambas reglas
tienen di¯cultad en la u
¶ltima columna, debido al kink en x = ¡0:05: Si ¶este fuera en x = 0;
entonces cualquier regla trapezoidal con un n¶
umero de puntos impares podr¶³a computar la
integral exactamente.
Ejercicio 2.2. Demostrar que esto se cumple.
Tabla 2.1. Algunas integrales sencillas
Regla
N¶
umero de puntos
R1
0
1
x 4 dx
19
R1
0
ex dx
R1
¡1
max[0; x + 0:05]dx
Trapezoidal
4
0.7212
1.7342
0.6056
7
0.7664
1.7223
0.5583
10
0.7797
1.7200
0.5562
13
0.7858
1.7193
0.5542
3
0.6496
1.4662
0.4037
7
0.7816
1.7183
0.5426
11
0.7524
1.6231
0.4844
15
0.7922
1.7183
0.5528
0.8000
1.7183
0.5513
Simpson
Verdadera
En general, cuando queramos una buena aproximaci¶on de la integral, va a ser necesario
utilizar muchos puntos para alcanzar la convergencia con nuestra funci¶on deseada.
Si no queremos de¯nir una regla en particular, podemos aproximar la integral en Matlab/Octave con quad y quad8, y no tenemos que de¯nir cuantos puntos queremos aproximar.
R 2¼
Por ejemplo, queremos integrar la funci¶on I = 0 e¡x sin(10x)dx, entonces escribimos
>>f=inline('exp(-x).*sin(10*x)');
>>quad(f,0,2*pi)
ans = 0.098825
Esto lo podemos comprobar resolviendo la integral por partes que nos dar¶³a
I=¡
1 ¡x
e [sin(10x) + 10 cos(10x)]j2¼
0 ¼ 0:0988
101
Si deseamos resolver integrales con m¶etodos de ordenes altas, entonces en vez de usar
el comando quad('f',xmin,xmax) usamos quad8(¢). El primer comando usa la regla de
Simpson, mientras que quad8 usa cotas de Newton m¶as elevadas de ordenes. Tambi¶en
notemos que tanto Matlab como Octave ya tienen incluida el comando del trapecio a trav¶es
de trapz(x,y), en donde x y y son vectores.
Sin embargo, si queremos resolver la integral exactamente, escribimos en Matlab el
comando syms x t para decirle a Matlab que queremos integrar. Si queremos una integral inde¯nida, pues s¶olo escribimo la integral, por ejemplo int(-2*x/(1+x^2)^2) returns
1/(1+x^2). Pero si queremos que nos integre en un intervalo [a; b] escribimos por ejemplo,
int(x*log(1+x),0,1) returns 1/4.
int
Realizemos a continuaci¶on algunos ejemplos usando:
20
a) La regla compuesta trapezoidal con 2 subintervalos
b) La regla compuesta trapezoidal con 10 subintervalos
c) La regla de Simpson compuesta con 2 subintervalos
d) La regla de Simpson compuesta con 10 subintervalos
1.
2.
R4
2x dx
0
R 1 1+x
dx
0 1+x3
1. octave:1> Trap('ejem1',0,4,2)
h=2
ans = 25
octave:2> Trap('ejem1',0,4,10)
h = 0.40000
ans = 21.779
octave:3> Simp('ejem1',0,4,2)
h=2
J = 22
ans = 22
octave:4> Simp('ejem1',0,4,10)
h = 0.40000
J = 21.641
ans = 21.641
octave:5> Trap('ejem2',0,1,2)
h = 0.50000
ans = 1.1667
octave:6> Trap('ejem2',0,1,10)
h = 0.10000
ans = 1.2075
21
octave:7> Simp('ejem2',0,1,2)
h = 0.50000
J = 1.2222
ans = 1.2222
octave:8> Simp('ejem2',0,1,10)
h = 0.10000
J = 1.2092
ans = 1.2092
2.2
Diferenciaci¶
on num¶
erica
La diferenciaci¶on num¶erica se usa frecuentemente en problemas num¶ericos. Por ejemplo, en
optimizaci¶on y en problemas de ecuaciones no lineales vamos a usar aproximaciones ¯nitas de
gradientes, Hesianos y Jacobianos. De hecho, esto por lo general se computa num¶ericamente,
porque las soluciones anal¶³ticas son muy dif¶³ciles, adem¶as de que demoran mucho tiempo.
Las derivadas num¶ericas tambi¶en son muy importantes en algunos m¶etodos de ecuaciones
en diferencia. Lo que haremos ser¶a encontrar estimaciones para la derivada o pendiente de
una funci¶on usando los valores de la funci¶on s¶olo en un conjunto de puntos discretos.
La derivada est¶a de¯nida como
f (x + ") ¡ f(x)
"!0
"
f 0 (x) = lim
lo que sugiere la f¶ormula
: f (x + h) ¡ f (x)
f 0 (x) =
h
en donde h queremos que sea peque~
na.
Veamos ahora c¶omo diferenciar exactamente en Matlab.
Usando tambi¶en antes el comando syms x t escribimos an¶alogamente como lo hicimos
con la integral, pero usamos ahora el comando diff, por ejemplo diff(sin(x^2))returns
2*cos(x^2)*x.
22
2.2.1
Primeras derivadas
Diferencias adelantadas, atrasadas y centrales Las f¶ormulas m¶as sencillas se basan
enutilizar una l¶³nea recta para interpolar, dados los datos, es decir, que usan dos puntos
dados para estimar la derivada. Por lo que vamos a suponer que xi¡1 y xi+1 existen.
Sea f (xi¡1 ) = yi¡1 ; f(xi ) = yi y f(xi+1 ) = yi+1 y vamos a suponer que la distancia
entre las x's es constante, es decir, h = xi+1 ¡ xi = xi ¡ xi¡1 :
Entonces veamos nuestras tres f¶ormulas:
Diferencia adelanta
f 0 (xi ) ¼
yi+1 ¡ yi
f (xi+1 ) ¡ f (xi )
=
xi+1 ¡ xi
xi+1 ¡ xi
Diferencia atrasada
f 0 (xi ) ¼
f (xi ) ¡ f (xi¡1 )
yi ¡ yi¡1
=
xi ¡ xi¡1
xi ¡ xi¡1
Un enfoque m¶as balanceado da una aproximaci¶on de la derivada en xi usando los valores
de f (xi¡1 ) y f (xi+1 ): Tomando el promedio de las dos diferencias arriba mencionadas, nos
da la f¶ormula de la diferencia central.
Diferencia central
f 0 (xi ) ¼
yi+1 ¡ yi¡1
f(xi+1 ) ¡ f (xi¡1 )
=
xi+1 ¡ xi¡1
xi+1 ¡ xi¡1
En la gr¶a¯ca 2.7 podemos observar estas tres f¶ormulas de diferencia con una funci¶on
aleatoria f 0 (x)
23
10
8
6
4
2
0
1
2
Diferencia central
Diferencia atrasada
3
Diferencia adelantada
Gr¶a¯ca 2.7. Aproximaciones de f¶ormulas de diferencias
Para ilustrar estas tres f¶ormulas, consideremos los puntos (1,2), (2,4), (3,8), (4,16) y
(5,32), en donde (xi ; yi ); con i = 0; 1; 2; 3; 4: Entonces la diferencia adelantada para f 0 (x2 ) =
f 0 (3) con h = 1 nos dar¶³a
f 0 (x2 ) ¼
f (x3 ) ¡ f(x2 )
y3 ¡ y2
=
= 16 ¡ 8 = 8
x3 ¡ x2
1
Analogamente para la diferencia adelantada nos dar¶³a f 0 (x2 ) ¼ 4, y para la diferencia
central nos da f 0 (x2 ) ¼ 6:
Sin embargo, podemos utilizar estas f¶ormulas, pero con h's m¶as altas, por ejemplo con
h = 2 para la diferencia central nos dar¶³a f 0 (x2 ) ¼ 7:5: Los datos que acabamos de utilizar
son de la funci¶on y = f (x) = 2x , por lo que podemos comparar las estimaciones de las
derivadas con los verdaderos valores. En este caso f 0 (x) = 2x (log 2); as¶³ que f 0 (3) ¼ 5:544:
F¶
ormulas generales de tres puntos A veces interpolar una funci¶on a trav¶es de
un polinomio es mejor que a trav¶es de una recta, adem¶as de que puede utilizar m¶as puntos.
Por lo que vamos a ampliar las dos primeras f¶ormulas antes mencionadas.
Diferencia adelanta con tres puntos
f 0 (xi ) ¼
¡f (xi+2 ) + 4f (xi+1 ) ¡ 3f (xi )
¡yi+2 + 4yi+1 ¡ 3yi
=
xi+2 ¡ xi
xi+2 ¡ xi
24
Diferencia atrasada con tres puntos
f 0 (xi ) ¼
3f (xi ) ¡ 4f(xi¡1 ) + 2f (xi¡2 )
3yi ¡ 4yi¡1 + 2yi¡2
=
xi ¡ xi¡2
xi ¡ xi¡2
Aplicadas al ejemplo anterior vemos que con la diferencia adelantada de tres puntos
0
f (x2 ) ¼ 4; y la diferencia atrasada de tres puntos nos da f 0 (x2 ) ¼ 5; por lo que vemos que
entre m¶as puntos utilicemos, m¶as nos acercamos al verdadero valor.
Empero, estas aproximaciones requieren que la distancia entre las x's sea constante,
pero no siempre este es el caso, por lo que el polinomio interpolador de Lagrange s¶olo necesita
que x1 < x2 < x3 :
Recordemos la forma general del polinomio que pasa por n puntos (x1 ; y1 ); :::; (xn ; yn )
por lo que tiene n t¶erminos que corresponden a cada un de los puntos:
p(x) = L1 y1 + L2 y2 + ::: + Ln yn
en donde
Lk (x) =
(x ¡ x1 ):::(x ¡ xk¡1 )(x ¡ xk+1 ):::(x ¡ xn )
(xk ¡ x1 ):::(xk ¡ xk¡1 )(xk ¡ xk+1 ):::(xk ¡ xn )
El numerador es el producto
Nk (x) = (x ¡ x1 ):::(x ¡ xk¡1 )(x ¡ xk+1 ):::(x ¡ xn )
por lo que podemos escribir
p(x) = c1 N1 + c2 N2 + ::: + cn Nn
en donde ck son los coe¯cientes del polinomio interpolador de Lagrange
ck =
yk
(xk ¡ x1 ):::(xk ¡ xk¡1 )(xk ¡ xk+1 ):::(xk ¡ xn )
En nuestro an¶alisis con tres puntos el polinomio interpolador de Lagrange ser¶³a:
L(x) = L1 (x)y1 + L2 (x)y2 + L3 (x)y3
La aproximaci¶on a la primera derivada de f viene de f 0 (x) ¼ L0 (x); que se puede
25
escribir como
L0 (x) = L01 (x)y1 + L02 (x)y2 + L03 (x)y3
en donde
2x ¡ x2 ¡ x3
(x1 ¡ x2 )(x1 ¡ x3 )
2x ¡ x1 ¡ x3
L02 (x) =
(x2 ¡ x1 )(x2 ¡ x3 )
2x ¡ x1 ¡ x2
L03 (x) =
(x3 ¡ x1 )(x3 ¡ x2 )
L01 (x) =
Por lo tanto,
f 0 (x) ¼
2x ¡ x2 ¡ x3
2x ¡ x1 ¡ x3
2x ¡ x1 ¡ x2
y1 +
y2 +
y3
(x1 ¡ x2 )(x1 ¡ x3 )
(x2 ¡ x1 )(x2 ¡ x3 )
(x3 ¡ x1 )(x3 ¡ x2 )
Por ejemplo, tenemos los puntos (-2,4), (0,2) y (2,8), entonces en un archivo .m escribimos los siguinetes comandos que nos ayudan a determinar los coe¯cientes del polinomio
interpolador de Lagrange y el resultado de dicho polinomio.
function p=Lagrange coef(x,y)
% Calcula los coeficientes de las funciones de Lagrange
n=length(x);
for k=1:n
d(k)=1;
for i=1:n
if i~=k
d(k)=d(k)*(x(k)-x(i));
end
c(k)=y(k)/d(k)
end
end
function p=Lagrange eval(t,x,c)
% Evalua el polinomio interpolador de Lagrange en x=t
m=length(x);
l=length(t);
26
for i=1:l
p(i)=0;
for j=1:m
N(j)=1;
for k=1:m
if (j~=k)
N(j)=N(j)*(t(i)-x(k));
end
end
p(i)=p(i)+N(j)*c(j);
end
end
En el cual en este caso nos va a dar p(x) = x2 + x + 2 , con lo cual evaluado en x = 2;
nos da p(2) = 8:
>> x=[-2 0 2]
x=
-2 0 2
>> y=[4 2 8]
y=
428
>> Lagrange coef(x,y)
c=
4
c=
-2
c=
0.5000
c=
0.5000 1.0000
c=
0.5000 1.0000
c=
0.5000 -0.5000
27
c=
0.5000 -0.5000 2.0000
c=
0.5000 -0.5000 1.0000
c=
0.5000 -0.5000 1.0000
>> Lagrange eval(2,x,[.5 -.5 1])
ans =
8
2.2.2
Segundas derivadas
Las f¶ormulas de derivadas m¶as altas se pueden encontrar diferenciando varias veces el polinomio interpolador o con las expansiones de Taylor. Por ejemplo, dados tres puntos con
distancias iguales entre ellos, la f¶ormula de la segunda derivada es
f 00 (xi ) ¼
1
[f(xi+1 ) ¡ 2f(xi ) + f (xi¡1 )]
h2
as¶³ que usando los puntos (2,4), (3,8), (4,16) con h = 1 evaluado en x2 = 3, tenemos
f 00 (3) ¼ [f (4) ¡ 2f (3) + f (2)] = 4
Con la ayuda del teorema de Taylor, podemos escribir las siguientes f¶ormulas de diferencia central:
000
f (xi ) ¼
0000
f (xi ) ¼
1
[f(xi+2 ) ¡ 2f(xi+1 ) + 2f (xi¡1 ) ¡ f (xi¡2 )]
2h3
1
[f (xi+2 ) ¡ 4f (xi+1 ) + 6f (xi ) ¡ 4f (xi¡1 ) ¡ f (xi¡2 )]
h4
28
3
Resolver un sistema de ecuaciones no lineales con
m¶
etodos de gradiente
Algunos conceptos de equilibrio est¶an expresados como sistemas de ecuaciones no lineales.
Estos problemas generalmente toman dos formas: ceros y puntos ¯jos. Si f : Rn ! Rn ,
entonces un cero de f es cualquier x tal que f (x) = 0, y un punto ¯jo de f es cualquier x
tal que f (x) = x: Estos son b¶asicamente el mismo poroblema, ya que x es un punto ¯jo de
f (x) si y s¶olo si (sii) es un cero de f(x) ¡ x: En este cap¶³tulo veremos m¶etodos num¶ericos
para resolver ecuaciones no lineales.
El concepto de equilibrio general de Arrow-Debreu consiste esencialmente en encontrar
un vector de precios tal que el exceso de demanda es cero, y este en efecto es un problema de
un sistema de ecuaciones no lineales. Otro ejemplo es el de equilibrios de Nash con estrategias
continuas y las trayectorias de transici¶on de los sistemas deterministicos din¶amicos. As¶³ que
como pueden ver en econom¶³a existen varios problemas de estos y nosotros lo que queremos
es tratar de resolverlos con m¶etodos num¶ericos.
Primero veremos m¶etodos para resolver problemas en una dimensi¶on y despu¶es veremos
para resolver en dimensi¶on ¯nita. Como suele suceder con estos m¶etodos, es que no existe
el m¶etodo perfecto, sini que cada uno tiene sus ventajas y sus desventajas, adem¶as de que
unos son mejores para determinados modelos que otros. Para entender mejor la aplicaci¶on de
estos m¶etodos en econom¶³a, veremos algunos ejemplos interesantes, aunque todav¶³a sencillos.
3.1
Problemas unidimensionales
Sea f : R ! R y queremos resolver el problema f(x) = 0: Este caso es muy importante, ya
que es la base para muchos m¶etodos multidimensionales. A los problemas para encontrar
cero de la funci¶on no lineal, tambi¶en se le llama ra¶³ces de una ecuaci¶on no lineal.
3.1.1
Bisecci¶
on
El m¶etodo de la bisecci¶on es una t¶ecnica sistem¶atica de b¶
usqueda para encontrar el cero
de una funci¶on continua. Este m¶etodo se basa en encontrar un intervalo en el cual un cero
se sabe que existe, dividiendo el intervalo en dos subintervalos iguales y determinar cual
subintervalo contiene al cero.
Sea f continua con f (a) < 0 < f (b) para alguna a; b; a < b, como se puede ver en la
gr¶a¯ca 3.1.
29
f(x)
a
c
b
Gr¶a¯ca 3.1. M¶etodo de la bisecci¶on
Dadas estas condiciones, el teorema del valor medio nos dice que 9f(») = 0; » 2 (a; b),
entonces el m¶etodo de la bisecci¶on usa este resultado varias veces para computar un cero.
Consideremos c =
1
(a
2
+ b); es decir el punto intermedio de [a; b]. Si f (c) = 0; pues ya
tenemos nuestro resultado. Si f (c) > 0, como en la gr¶a¯ca, entonces existe un cero en
(a; c) y el m¶etodo de la bisecci¶on contin¶
ua con (a; c): An¶alogo para el caso en que f(c) < 0:
Entonces lo que hace el m¶etodo de la bisecci¶on es que va a trabajar cada vez en intervalos
m¶as peque~
nos, aunque cabe se~
nalar que puede haber varios ceros, y el objetivo este m¶etodo
es encontrar un cero. Sin embrago, podemos usar este m¶etodo varias veces para encontrar
m¶as ceros.
Veamos a continuaci¶on el algor¶³tmo de este m¶etodo:
Objetivo: Encontrar f (x) = 0; f : R ! R
Paso 1: Encontrar xI < xD tal que f (xI )f (xD ) < 0 y escoger la tolerancia "; ± > 0
Paso 2: Computar el punto medio xM =
xI +xD
2
M
I
Paso 3: Re¯nar los l¶³mites, es decir, si f (x )f (x ) < 0; entonces xD = xM y no cambia
xI , en otro caso xI = xM y no cambiar xD
Paso 4: Veri¯car si xD ¡ xI · "(j1 + jxI j + jxD j) ¶o jf (xM )j · ±, entonces detener y
reportar la soluci¶on en xM , en otro caso, volver al paso 1.
Mientras que el m¶etodo de la bisecci¶on es simple, nos dice los componentes importantes
para resolver cualquier ecuaci¶on no lineal, es decir, hacemos un guess inicial, computamos
30
las iteraciones y veri¯camos que la u
¶ ltima iteraci¶on sea una soluci¶on aceptable. Este m¶etodo
tiene las ventajas de que es sencillo y siempre que existe una soluci¶on la vamos a encontrar,
sin embargo, el proceso puede ser lento.
Veamos un ejemplo del m¶etodo de la bisecci¶on. Queremos encontrar la aproximaci¶on
p
num¶erica de 3; entonces aproximamos el cero de y = f (x) = x2 ¡ 3: Dado que f (1) = ¡2
y f (2) = 1; entonces tomamos como puntos iniciales a = 1 y b = 2 (y(a) = ¡2 y y(b) = 1),
por lo que nuestro punto medio es x1 = 1:5 y evaluandolo en la funci¶on nos da y1 = f (1:5) =
¡0:75: Dado que y(a) y y1 tienen el mismo signo, y y(b) y y1 tienen signos opuestos, sabemos
que existe un cero en el intervalo [x1 ; b]: Entonces ahora tomamos a = 1:5 y continuando de
esta manera vamos a obtener el siguiente cuadro para los primeros cinco pasos.
Cuadro 3.1.C¶alculos de
Paso
a
b
p
3 usando el m¶etodo de la bisecci¶on
xi
y(a)
y(b)
yi
1
1
2
1.5
-2
1
-0.75
2
1.5
2
1.75
-0.75
1
0.0625
3
1.5
1.75
1.625
-0.75
0.0625
-0.3594
4
1.625
1.75
1.6875 -0.3594 0.0625
-0.1523
5
1.6875
1.75
1.7188 -0.1523 0.0625
-0.0459
Tenemos una estimaci¶on del cero en cada iteraci¶on, y sabemos que despu¶es de k iteraciones el error va a ser a lo mucho
b1 ¡a1
;
2k
en donde a1 y b1 es el intervalo original. Como
podemos observar, puede haber ocaciones en que estemos muy cerca del cero, pero que al
siguiente paso nos alejemos.
Veamos ahora el m¶etodo de la bisecci¶on en Matlab:
function [x,y] = bisect(fun,a,b,tol,max)
% Resuelve f(x)=0 usando el m¶
etodo de la bisecci¶
on
% fun
% [a,b]
Nombre de la funci¶
on entre ap¶
ostrofes
Intervalo que contiene al cero
% tol
Tolerancia permitida al computar el cero
% maxit
N¶
umero m¶
aximo de iteraciones permitido
% x
Vector de aproximaciones de cero
% y
Vector de valores de la funci¶
on fun(x)
31
a(1)=a;
b(1)=b;
ya(1)=feval(fun,a(1));
yb(1)=feval(fun,b(1));
if ya(1)*yb(1)>0
error( ' La funci¶
on tiene el mismo signo en sus puntos extremos')
end
for i = 1:max
x(i) = (a(i) + b(i))/2;
y( i ) = feval (fun , x(i));
if ((x(i)-a(i)) < tol )
disp( ' El m¶
etodo de la bisecci¶
on ha convergido');
break;
end
if y(i) == 0
disp('se encontr¶
o el cero exacto');
break;
elseif y(i)*ya(i)<0
a(i+1)=a(i);
ya(i+1)=ya(i);
b(i+1)=x(i);
yb(i+1)=y(i);
else
a(i+1)=x(i);
ya(i+1)=y(i);
b(i+1)=b(i);
yb(i+1)=yb(i);
end;
iter = i;
end
if (iter >= max)
disp('no se encontr¶
o el cero a la tolerancia deseada');
end
n= length(x);
32
k = 1:n;
out = [k' a(1:n)' b(1:n)' x' y'];
disp(' paso a b x y')
disp(out)
Para el caso que realizamos manualmente veamos la funci¶on:
function y=fun1(x)
y=x^2-3;
Con lo cual en la pantalla de Matlab s¶olo escribimos Bisect('fun1',1,2,0.1,5).
>> Bisect('fun1',1,2,0.1,5)
El m¶etodo de la bisecci¶on ha convergido
paso a b x y
1.0000 1.0000 2.0000 1.5000 -0.7500
2.0000 1.5000 2.0000 1.7500 0.0625
3.0000 1.5000 1.7500 1.6250 -0.3594
4.0000 1.6250 1.7500 1.6875 -0.1523
ans =
1.5000 1.7500 1.6250 1.6875
El m¶etodo de la bisecci¶on es lento, pero seguro; y su tasa de convergencia lineal es de
1
:
2
La ventaja de este m¶etodo es que puede funcionar para problemas que causen di¯cultad
con otros m¶etodos.
3.1.2
M¶
etodo de Newton
Otra desventaja del m¶etodo de la bisecci¶on es que s¶olo asume que f es continua. El m¶etodo
de Newton usa las propiedades de suavidad de f para formular un m¶etodo que es r¶apido
cuando funciona, pero que puede no siempre convergir. El m¶etodo de Newton aproxima a
f (x) con una sucesi¶on de funciones lineales y aproxima un la nulidad de f con los ceros de
las aproximaciones lineales. Si f es suave, entonces estas aproximaciones son crecientemente
precisas, y las iteraciones sucesivas converger¶an r¶apido al cero de f. B¶asicamente el m¶etodo
de Newton reduce un problema no lineal a una sucesi¶on de problemas lineales, en donde los
ceros son f¶aciles de computar.
33
Formalmente, el m¶etodo de Newton es una simple interaci¶on. Suponer que nuestro
primer guess es xk . Entonces en xk construimos la aproximaci¶on lineal de f en xk ; dando
como resultado la funci¶on
g(x) ´ f 0 (xk )(x ¡ xk ) + f(xk )
Las funciones f (x) y g(x) son tangentes en xk y generalmente cerca en una vecindad de xk :
Entonces en vez de resolver la nulidad de la funci¶on f , resolvemos la nulidad de la funci¶on
g, esperando que las dos tengan nulidades similares. Nuestro nuevo guess, xk+1 ; ser¶a el cero
de g(x), implicando que
xk+1 = xk ¡
f (xk )
f 0 (xk )
(3.1)
La siguiente gr¶a¯ca muestra los pasos de x1 a x5 ; en donde se ve que el nuevo guess no
alcanza la nulidad de la funci¶on, pero esperamos que la secuencia xk converja a la nulidad
de f . El teorema 3.1 provee las condiciones su¯cientes para la convergencia.
Gr¶a¯ca 3.2. M¶etodo de Newton
00
Teorema 3.1. Sea f 2 C 2 y f (x¤ ) = 0: Si x0 est¶a su¯cientemente cerca a x¤ ; f 0 (x¤ ) 6= 0
¤
(x )
¤
y j ff 0 (x
¤ ) j < 1, entonces la secuencia de Newton xk de¯nida en (3.1) converge a x ; y es
cuadr¶aticamente convergente, es decir,
lim sup
k!1
jxk+1 ¡ x¤ j
<1
jxk ¡ x¤ j2
Escoger entre estos dos m¶etodos no es f¶acil, ya que por una parte el m¶etodo de la
bisecci¶on, aunque lento, siempre converge a cero, mientras que el m¶etodo de Newton es
34
r¶apido, pero puede no converger.
Veamos a continuaci¶on el algor¶³tmo de este m¶etodo:
Objetivo: Encontrar f (x) = 0; f : R ! R; f 2 C 1
Paso 1: Escoger "; ± y el punto inicial x0 , es decir, k = 0
Paso 2: Computar la siguiente iteraci¶on xk+1 = xk ¡
f (xk )
f 0 (xk )
Paso 3: Veri¯car las condiciones para detenerse: si jxk ¡ xk+1 j · "(1 + jxk+1 j); pasar
al siguiente paso, en otro caso pasar al paso 1
Paso 4: Reportar los resultados y detenerse: Si jf (xk+1 )j · ±; reportar ¶exito para
encontrar el cero, en otro caso, reportar fracaso
q
>Que pasa si queremos encontrar 34 ? Lo que hacemos es aproximar el cero de y =
f (x) = 4x2 ¡ 3, como se ve en la gr¶a¯ca 3.3. Para esto vamos a necesitar y 0 = f 0 (x) = 8x:
Dado que f (0) = ¡3 y f (1) = 1, usamos como nuestro guess inicial x0 = 0:5, con lo cual
su corespondiente valor de la funci¶on es y0 = ¡2 y el valor de la derivada es y00 = 4:
1
0.2
0.4
x
0.6
0.8
1
0
-1
-2
-3
-4
Gr¶a¯ca 3.3. y = 4x2 ¡ 3 y la aproximaci¶on de la recta tangente en x = 0:5
Por lo tanto nuestra primera aporximaci¶on del cero es
x1 = x0 ¡
¡2
f(x0 )
= 0:5 ¡
=1
0
f (x0 )
4
Continuando con un paso m¶as vamos a tener
x2 = x1 ¡
y1
1
7
= 1 ¡ = = 0: 875
0
y1
8
8
En Matlab podemos computar el m¶etodo de Newton de la siguiente manera:
function [x,y]=Newton(fun,fun pr,x1,tol,max);
35
% Resuelve f(x)=0 usando el m¶
etodo de Newton
% fun Nombre de la funci¶
on entre ap¶
ostrofes
% funpr Especifica la derivada de f
% x0 Estimaci¶
on inicial
% tol Tolerancia permitida al computar el cero
% maxit N¶
umero m¶
aximo de iteraciones permitido
%
% x Vector de aproximaciones de cero
% y Vector de valores de la funci¶
on fun(x)
x(1)=x1;
y(1)=feval(fun,x(1));
ypr(1)=feval(funpr,x(1));
for i=2:max
x(i)=x(i-1)-y(i-1)/ypr(i-1);
y(i)=feval(fun,x(i));
if abs(x(i)-x(i-1))<tol
disp('El metodo de Newton converge');
break;
end
ypr(i)=feval(funpr,x(i));
iter=i;
end
if (iter>=max)
disp('no se encontr¶
o el cero a la tolerancia deseada');
end
n=length(x);
k=1:n;
out=[k' x' y'];
disp(' paso x y')
disp(out)
Para el caso que realizamos manualmente veamos la funci¶on:
function y=fun2(x)
y=4*x^2-3;
36
con su correspondiente derivada:
function y=funpr2(x)
y=8*x;
Con lo cual en la pantalla de Matlab s¶olo escribimos Newton('fun2','funpr2',0.5,0.1,10).
>> Newton('fun2','funpr2',0.5,0.1,10)
El metodo de Newton converge
paso x y
1.0000 0.5000 -2.0000
2.0000 1.0000 1.0000
3.0000 0.8750 0.0625
4.0000 0.8661 0.0003
ans =
0.5000 1.0000 0.8750 0.8661
3.1.3
M¶
etodo de la secante
Un aspecto clave para el m¶etodo de Newton es computar f 0 (x); el cual puede ser costoso.
El m¶etodo de la secante emplea aproximaciones lineales, pero nunca eval¶
ua f 0 ; sino que
aproxima f 0 (xk ) con la pendiente de la secante de f entre xk y xk+1 , resulatndo de esta
manera en la iteraci¶on
xk+1 = xk ¡
f (xk )(xk ¡ xk¡1 )
f (xk ) ¡ f (xk¡1 )
(3.2)
Este m¶etodo es similar al de Newton, s¶olo que en el algor¶³tmo en el paso 2 se utiliza la
ecuaci¶on (3.2) en vez de (3.1) para computar la siguiente iteraci¶on. El m¶etodo de la secante
tiene los mismos problemas de convergencia que el m¶etodo de Newton, s¶olo que puede ahorrar
tiempo debido a que no eval¶
ua f: La tasa de convergencia es entre lineal y cuadr¶atica.
Teorema 3.2. Sea f (x¤ ) = 0; f 0 (x¤ ) 6= 0 y f 0 (x) y f 00 (x) son continuas alrededor de
x¤ ; entonces el m¶etodo de la secante converge a una tasa
lim sup
k!1
jxk+1 ¡ x¤ j
jxk ¡ x¤ j
37
p
(1+ 5)
2
p
(1+ 5)
;
2
<1
es decir,
Usando el ejemplo de y = f (x) = x2 ¡ 3, veamos que tal funciona para el m¶etodo de la
secante. Las extremos que utilizaremos son x0 = 1 y x1 = 2; con lo cual y0 = ¡2 y y1 = 1;
respectivamente.
Nuestra primera aproximaci¶on al cero es
x2 = x1 ¡
2¡1
5
x1 ¡ x0
y1 = 2 ¡
(1) = ¼ 1:667
y1 ¡ y0
1 ¡ (¡2)
3
El m¶etodo de la secante no requiere que los puntos usados para computar la pr¶oxima
aproximaci¶on envuelvan al cero, ni incluso que x0 < x1 : Es conveniente ver los c¶alculos para
ver como est¶an evolucionando los resultados, los cuales los podemos ver en el cuadro 3.2.
Cuadro 3.2.C¶alculos de
Paso
p
3 usando el m¶etodo de la secante
x
y
1
1
-2
2
2
1
3
1.6667 -0.2222
4
1.7273 -0.0165
5
1.7321
6
1.7321 -0.0000
0.0003
Veamos ahora la funci¶on en Matlab para el m¶etodo de la secante:
function [x,y] = secant(fun,a,b,tol,max)
% Resuelve f(x)=0 usando el m¶
etodo de la secante
% fun Nombre de la funci¶
on entre ap¶
ostrofes
% a,b primeros dos estimadores de cero
% tol Tolerancia permitida al computar el cero
% max N¶
umero m¶
aximo de iteraciones permitido
% x Vector de aproximaciones de cero
% y Vector de valores de la funci¶
on fun(x)
x(1)=a;
x(2)=b;
y(1)=feval(fun,x(1));
y(2)=feval(fun,x(2));
for i = 2:max
38
x(i+1)=x(i)-y(i)*(x(i)-x(i-1))/(y(i)-y(i-1));
y(i+1)=feval(fun,x(i+1));
if abs(x(i+1)-x(i)) < tol
disp( ' El m¶
etodo de la secante ha convergido');
break;
end
if y(i) == 0
disp('se encontr¶
o el cero exacto');
break;
end
iter=i;
end
if (iter>=max)
disp('no se encontr¶
o el cero a la tolerancia deseada');
end
n= length(x);
k = 1:n;
out = [k' x' y'];
disp(' paso x y')
disp(out)
Para el caso que realizamos manualmente veamos la funci¶on:
function y=fun1(x)
y=x^2-3;
Con lo cual en la pantalla de Matlab s¶olo escribimos secant('fun1',1,2,0.001,10).
El m¶etodo de la secante ha convergido
paso
x
y
1.0000 1.0000 -2.0000
2.0000 2.0000 1.0000
3.0000 1.6667 -0.2222
4.0000 1.7273 -0.0165
5.0000 1.7321 0.0003
6.0000 1.7321 -0.0000
39
ans =
1.0000 2.0000 1.6667 1.7273 1.7321 1.7321
3.1.4
Iteraci¶
on del punto ¯jo
Como con las ecuaciones lineales, generalmente vamos a poder reescribir las problemas no
lineales en formas que nos sugieran un m¶etodo computacional. Cualquier problema de punto
¯jo x = g(x) sugiere la iteraci¶on xk+1 = g(xk ): Por ejemplo, consideremos x3 ¡ x ¡ 1 = 0; el
1
cual se puede reescribir en la forma del punto ¯jo x = (x + 1) 3 ; que nos sugiere la iteraci¶on
1
xk+1 = (xk + 1) 3 ; la cual converge a una soluci¶on si x0 = 1: Sin embargo, si reescribimos
x3 ¡ x ¡ 1 = 0 como x3 ¡ 1 = x; la forma sugerida ser¶a xk+1 = x3k ¡ 1, la cual diverge a ¡1
cuando x0 = 1:
Veamos un caso sencillo, por ejemplo, queremos sacar
p
3 con la iteraci¶on del punto
¯jo. Escribimos x2 = 3; y la reescribimos como
c
1
x = (x + )
2
x
que nos da la base para la iteraci¶on, es decir que nos va a quedar algo de la forma:
1
c
xk = (xk¡1 +
)
2
xk¡1
Geometricamente hablando lo que queremos es la intersecci¶on entre y = 21 (x+ xc ) = g(x)
p
con y = x: Para el caso de 3 tendr¶³amos x1 = 21 (1 + 31 ) = 2 y xk = 21 (2 + 23 ) = 74 = 1: 75 si
usamos el punto inicial x0 = 1; como se ve en la siguiente gr¶a¯ca.
40
3
2.5
2
1.5
1
0.5
0
0.5
1
1.5
x
2
2.5
3
Gr¶a¯ca 3.4. 1er paso de la iteraci¶on del punto ¯jo de x = 21 (x + x3 )
Hacer la gr¶a¯ca a mano de una divergente, con una media parabola hacia abajo, pero
muy curveada.
La funci¶on en Matlab ser¶³a:
function [x,gval] = fixpoint(g,x,tol,max)
% Computa el punto fijo de una funci¶
on
% g Nombre de la funci¶
on de la forma gval=g(x)
% x Guess inicial para el punto fijo
% x Punto fijo de g
% gval Estimaci¶
on del valor de la funci¶
on
for it=1:max
gval = feval(g,x);
if norm(gval-x)<tol, return, end
x = gval;
end
warning('Fracaso')
Para el caso que realizamos manualmente veamos la funci¶on:
function y=g(x)
y=(x+(3/x))/2;
Con lo cual en la pantalla de Matlab s¶olo escribimos fixpoint('g',1,0.01,10).
>> ¯xpoint('g',1,0.01,10)
41
ans =
1.7321
Este m¶etodo es m¶as o menos utilizado, sin embargo generalmente no se puede con¯ar
mucho en ¶el. Veremos m¶etodos m¶as con¯ables, aunque esto no quiere decir que debemos
ignorar estos esquemas de iteraciones de punto ¯jo, s¶olo hay que tener mucho cuidado y estar
pendientes de sus desventajas.
3.2
Problemas multidimensionales
Muchos problemas tienen varias ecuaciones, por lo que nosotros vamos a querer resolver el
problema f (x) = 0; con f : Rn ! Rn , es decir,
f 1 = (x1 ; x2 ; :::; xn ) = 0
..
.
(3.3)
f n = (x1 ; x2 ; :::; xn ) = 0
Pero esto de extender algunos m¶etodos que vimos en la secci¶on anterior, no es tan
f¶acil para algunos casos. Por ejemplo, para dos funciones de dos variables, z = f (x; y) y
z = g(x; y), encontrar el cero del sistema requiere que encontremos la intersecci¶on de las
curvas f (x; y) = 0 y g(x; y) = 0: Con lo cual es muy importante utilizar toda la informaci¶on
posible para determinar cual es la regi¶on en donde pueda estar el cero.
3.2.1
M¶
etodo de Newton
Para resolver (3.3) con el m¶etodo de Newton, vamos a reemplazar f con una aproximaci¶on
lineal, y posteriormente resolver el problema lineal para generar el nuevo guess; tal como lo
hicimos para el caso unidimensional. Por el teorema de Taylor, la aproximaci¶on lineal de f
alrededor de guess inicial x0 es
:
f(x) = f (x0 ) + J(x0 )(x ¡ x0 )
Podemos resolver la nulidad de esta aproximaci¶on lineal, dando como resultado x1 =
x0 ¡ J(x0 )¡1 f (x0 ): Este cero luego sirve como el nuevo guess sobre el cual vamos a volver a
42
linealizar. La iteraci¶on de Newton es
xk+1 = xk ¡ J(xk )¡1 f (xk )
Veamos c¶omo ser¶³a el algor¶³tmo de este m¶etodo:
Paso 1: Escoger el criterio para detenerse " y ± y el punto inicial x0 , es decir, k = 0
Paso 2: Computar la siguiente iteraci¶on: Computar el jacobiano Ak = J(xk ); resolver
Ak sk = ¡f (xk ) para sk ; y considerar xk+1 = xk + sk
Paso 3: Veri¯car los criterios para detenernos: Si jjxk ¡ xk+1 jj · "(1 + jjxk+1 jj) seguir
al paso 4; en otro caso ir al paso 2.
Paso 4: Detenerse y reportar los resultados: Si jjf (xk+1 )jj · ±; reportar ¶exito al
encontrar el cero; en otro caso reportar fracaso.
function x=Newton sys(F,JF,x0,param,tol,maxit)
% Resuelve el sistema no lineal F(x)=0 usando el m¶
etodo de Newton
% Los vectores x y x0 son vectores fila
% la funci¶
on F da un vector columna [f1(x)...fn(x)]'
% detener cuando la norma del cambio en el vector soluci¶
on sea menor a la
tolerancia
% la siguiente aproximaci¶
on de soluci¶
on es x new=x old+y';
x old=x0;
disp([0 x old]);
iter=1;
while (iter<=maxit)
y=-feval(JF,x old)nfeval(F,x old,param);
x new=x old+y';
dif=norm(x new-x old);
disp([iter x new dif]);
if dif<=tol
x=x new;
disp('El m¶
etodo de Newton ha convergido')
return;
else
x old=x new;
end
43
iter=iter+1;
end
disp('El m¶
etodo de Newton no convergi¶
o')
x=x new;
Usemos por ejemplo el sistema de ecuaciones:
f (x; y; z) = x3 ¡ 10x + y ¡ z + 3 = 0
g(x; y; z) = y 3 + 10y ¡ 2x ¡ 2z ¡ 5 = 0
h(x; y; z) = x + y ¡ 10z + 2 sin(z) + 5 = 0
cuyo Jacobiano es
2
3x2 ¡ 10 1
¡1
1
¡10 + 2 cos(z)
6
J(x) = 6
4 ¡2
3y 2 + 10 ¡2
1
3
7
7
5
En Matlab escribimos el sistema de funciones con su respectivo Jacobiano:
function F=F1(x,param)
F=[x(1)^3-10*x(1)+x(2)-x(3)+3-param(1);...
x(2)^3+10*x(2)-2*x(1)-2*x(3)-5-param(2);...
x(1)+x(2)-10*x(3)+2*sin(x(3))+5-param(3)];
function JF=JF1(x)
JF=[(3*x(1)^2)-10 1 -1; -2 (3*x(2)^2)+10 -2; 1 1 -10+2*cos(x(3))];
As¶³ que en Matlab s¶olo escriben Newton sys('F1','JF1',[1 1 1],[0 0 0],0.01,10)
>> Newton sys('F1','JF1',[1 1 1],[0 0 0],0.01,10)
0111
1.0000 0.1359 0.6699 0.7185 0.9669
2.0000 0.2955 0.6745 0.7305 0.1601
3.0000 0.2970 0.6748 0.7307 0.0015
El m¶etodo de Newton ha convergido
ans =
0.2970 0.6748 0.7307
44
3.2.2
M¶
etodo de la secante
Como vimos para el caso univariado, el m¶etodo de la secante es similar al de Newton, de
hecho, algunos lo llaman el m¶etodo de Newton-Raphson, s¶olo que utiliza aproximaciones
lineales para el Jacobiano, en vez del Jacobiano per se.
Veamos que el ejemplo anterior tambi¶en se puede resolver con este m¶etodo.
function x = secant sys(fun,x0,param,tol,maxit);
% secant sys.m Programa para resolver sistemas de ecuaciones no lineales
% detener cuando la norma del cambio en el vector soluci¶
on sea menor a la
tolerancia
del = diag(max(abs(x0)*1e-4,1e-8));
n = length(x0);
for i=1:maxit
f = feval(fun,x0,param);
for j=1:n
J(:,j) = (f-feval(fun,x0-del(:,j),param))/del(j,j);
end;
x = x0-inv(J)*f;
if norm(x-x0)<tol
break
end
x0 = x;
end
if i>=maxit
sprintf('Se alcanzaron %iteraciones',maxit)
end
Entonces en la ventana de Matlab escribimos
secant sys('F1',[1 1 1]',[0 0 0],0.001,10)
ans =
0.2970
0.6748
0.7307
45
3.2.3
Iteraci¶
on del punto ¯jo
Por analog¶³a al caso unidimensional, el procedimiento m¶as sencillo para resolver x = f (x) es
xk+1 = f (xk ): A este m¶etodo tambi¶en se le llama aproximaci¶on sucesiva, sustituci¶on sucesiva
o iteraci¶on funci¶on. En el caso multivariado, s¶olo se generaliza el proceso univariado.Un
ejemplo del punto ¯jo est¶a en el teorema de la contracci¶on (visto en el curso de Programaci¶on
Din¶amica).
Por ejemplo, si tenemos el sistema de ecuaciones
f1 (x1 ; x2 ) = x31 + 10x1 ¡ x2 ¡ 5 = 0
f2 (x1 ; x2 ) = x1 + x32 ¡ 10x2 + 1 = 0
lo reescribimos de la forma x1 = g1 (x1 ; x2 ) y x2 = g2 (x1 ; x2 ) de la siguiente manera:
x1 = ¡0:1x31 ¡ 0:1x2 ¡ 0:5
x2 = 0:1x1 + 0:1x32 + 0:1
La siguiente gr¶a¯ca ilustra f1 = 0 y f2 = 0, y el cuadro muestra las soluciones aproximadas
de cada iteraci¶on
1
f1=0
0.8
0.6
0.4
f2=0
0.2
0
0.2
0.4
0.6
0.8
Gr¶a¯ca 3.5.Sistema de dos ecauciones no lineales
La funci¶on escrita en lenguaje Matlab ser¶³a:
46
1
function x=fixpoint sys(G,x0,tol,max)
% Resuleve el sistema no lineal x=G(x) usando la iteraci¶
on del punto fijo
% Los vectores x y x0 son vectores fila
% la funci¶
on G da un vector columna [g1(x)...gn(x)]'
% detener cuando la norma del cambio en el vector soluci¶
on sea menor a la
tolerancia
% la siguiente aproximaci¶
on de soluci¶
on es x new=x old+y';
disp([0 x0]);
x old=x0;
iter=1;
while (iter<=max)
y=feval(G,x old)
x new=y';
dif=norm(x new-x old);
disp([iter x new dif]);
if dif<=tol
x=x new;
disp('La iteraci¶
on del punto fijo ha convergido')
return;
else
x old=x new;
end
iter=iter+1;
end
disp('La iteraci¶
on del punto fijo no convirgi¶
o')
x=x new;
La funci¶on que estabamos utilizando la escribimos como
function G=ejem4(x)
G=[(-0.1*x(1)^3+0.1*x(2)+0.5)
(0.1*x(1)+0.1*x(2)^3+0.1)];
En la ventana de Matlab s¶olo escribimos fixpoint sys('ejem4',[1 1],0.01,10)
>> ¯xpoint sys('ejem4',[1 1],0.01,10)
011
47
1.0000 0.5000 0.3000 0.8602
2.0000 0.5175 0.1527 0.1483
3.0000 0.5014 0.1521 0.0161
4.0000 0.5026 0.1505 0.0020
La iteraci¶on del punto ¯jo ha convergido
ans =
0.5026 0.1505
3.3
Ejemplos econ¶
omicos
Veamos en este cap¶³tulo un ejemplo de Microeconom¶³a, en particular de teor¶³a de juegos.
Vamos a ilustrar un ejemplo del modelo de Cournot con el m¶etodo de la secante.
En el modelo vamos a tener que el inverso de la demanda de un bien es P (q) = q
¡1
´
y
las ¯rmas van a producir con las siguientes funciones de producci¶on
1
Ci (qi ) = ci qi2 ;
2
i = 1; 2
Las ganancias de la empresa i son
¼ i (q1 ; q2 ) = P (q1 + q2 )qi ¡ Ci (qi )
Si la empresa i toma como dado el producto de la otra empresa, va a optimizar su nivel
de producci¶on resolviendo
@¼ i
= P (q1 + q2 ) + P 0 (q1 + q2 )qi ¡ Ci0 (qi ) = 0
@qi
Entonces los productos de equilibrio, q1 y q2 ; son las ra¶³ces de las dos ecuaciones no
lineales
fi (q) = (q1 + q2 )
¡1
´
¡
µ
¶
¡1
¡1
(q1 + q2 ) ´ ¡1 qi ¡ ci qi = 0
´
para i = 1; 2:
Vamos a suponer que usamos ´ = 1:6; aunque notemos que ¶este tambi¶en lo podemos
poner como par¶ametro de la funci¶on de Cournot, como hemos hecho para los diferentes
costos.
48
function fval = cournot(q)
c = [param(1);param(2)];
eta = 1.6;
e = -1/eta;
fval = sum(q)^e + e*sum(q)^(e-1)*q - diag(c)*q;
As¶³ que solo escribimos en la ventana de Matlab secant sys('cournot',[0.2 0.2]',[0.6
0.8],0.001,10) ¶o secant sys('cournot',[0.2 0.2]',[0.6 0.6],0.001,10)
>> secant sys('cournot',[0.2 0.2]',[0.6 0.8],0.001,10)
ans =
0.8396
0.6888
>> secant sys('cournot',[0.2 0.2]',[0.6 0.6],0.001,10)
ans =
0.8329
0.8329
49
4
M¶
etodos para sistemas de ecuaciones lineales
Resolver ecuaciones lineales tambi¶en es muy utilizado en la pr¶actica, y a veces se necesita
la combinaci¶on de m¶etodos que resuleven sistemas de ecuaciones lineales con ecuaciones no
lineales, como al ¯nal de este cap¶³tulo veremos un ejemplo. Tambi¶en ya hemos visto que
algunos m¶etodos aproximan con ecuaciones lineales algunos problemas no lineales, por lo
que entender c¶omo se resuleven los problemas lineales y las di¯cultades asociadas es muy
importante para el an¶alisis num¶erico.
Los problemas del sistema de ecuaciones lineales, Ax = b, se divide basicamente en
dos tipos, dependiendo de las entradas de A: i) Decimos que A es denso si aij 6= 0 para
la mayor¶³a de i; j. ii) Decimos que A es dispersa o poco densa (sparse) si aij = 0 para la
mayor¶³a de i; j. N¶otemos que estas de¯niciones no son muy precisas, pero en la pr¶actica
vamos a ver que claramente podemos determinar si la matriz A es densa o no tanto.
Como vimos en el cap¶³tulo anterior, >qu¶e sucede si tenemos ecuaciones no lineales?
Tambi¶en en este caso vamos a tener varias inc¶ognitas, por lo que necesitamos tambi¶en necesitamos m¶etodos multidimensionales. Supongamos de nuevo que f : Rn ! Rn y que queremos
resolver f (x) = 0; con lo cual tenemos lo mismo que en (3.3)
f 1 = (x1 ; x2 ; :::; xn ) = 0
..
.
f n = (x1 ; x2 ; :::; xn ) = 0
Vamos a ver dos m¶etodos que se pueden aplicar tanto a problemas lineales como no
lineales: Gauss-Jacobi y Gauss-Seidel.
4.1
4.1.1
M¶
etodos directos
Eliminaci¶
on Gaussiana
La eliminaci¶on Gaussiana es un m¶etodo directo com¶
un para resolver un sistema lineal. La
eliminaci¶on Gaussiana se basa en el hecho de que si dos ecuaciones tienen un punto en com¶
un,
entonces este punto tambi¶en satisface cualquier combinaci¶on lineal de esas dos ecuaciones.
Si podemos encontrar combinaciones lineales de una forma m¶as o menos simple, entonces
podremos encontrar la soluci¶on del problema original de una forma m¶as sencilla. La forma
50
que estar¶³amos buscando ser¶³a aquella en la que ciertas variables han sido eliminadas de
algunas ecuaciones.
El caso m¶as simple es el de las matrices triangulares, en donde se dice que A es triangular baja si todos los elementos que no son cero est¶an sobre o bajo la diagonal, es decir,
0
B
B
B
A=B
B
@
a11 0
a21
..
.
:::
0
a22 ::: 0
..
...
0
.
an1 an2 :::
ann
1
C
C
C
C
C
A
An¶alogamente para el caso en que A sea triangular alta. Una matriz diagonal es aquella
que s¶olo tienen elementos que no son cero sobre la diagonal. Recordemos que una matriz
triangular es no singular si y s¶olo si todos los elementos de la diagonal no son cero. Los
sistemas lineales en los cuales A es triangular se pueden resolver con sustituci¶on hacia atr¶as.
Supongamos que A es triangular baja y no singular, dado que todos los elementos que est¶an
fuera de la diagonal de la primera ¯la son cero, el problema se reduce a resolver a11 x1 = b1 ;
cuya soluci¶on sabemos que es x1 =
b1
:
a11
Con esta soluci¶on podemos resolver para x2 , que
implica a22 x2 + a21 x1 = b2 ; y as¶³ sucesivamente en donde tendremos la sucesi¶on
x1 =
xk =
b1
a11
bk ¡
k¡1
P
akj xj
j=1
para k = 2; 3; :::; n
akk
que est¶a bien de¯nida para matrices no singulares. Si la matriz fuera triangular alta entonces
empecar¶³amos con xn =
bn
:
ann
Pero lamentablemente en la pr¶actica rara vez nos encontraremos con ecuaciones en esta
forma, pero lo que vamos a hacer ser¶a buscar poner nuestro sistema de ecuaciones de esta
manera.
Por ejemplo, suponer que tenemos el siguiente sistema de ecuaciones:
2x + 6y + 10z = 0
x + 3y + 3z = 2
3x + 14y + 28z = ¡8
51
que en forma matricial la podemos escribir as¶³:
3
2
2 6 10 j 0
7
6
6 1 3 3 j 2 7
5
4
3 14 28 j ¡8
donde recordando un poco de algebra matricial podemos realizar nuestro primer paso de la
eliminaci¶on, con lo cual nos quedar¶³a:
3
2
2 6
10
j 0
7
6
6 0 0
4
j ¡4 7
5
4
0 ¡10 ¡26 j 16
normalizando y cambiando ¯las tenemos
2
1 3 5
3
j 0
7
6
6 0 1 2:6 j ¡1:6 7
5
4
0 0 1
j ¡1
que ya es una matriz triangular alta, por lo que podemos dar los resultados de una vez
sustituyendo en las ecuaciones, o podr¶³amos continuar con este m¶etodo para que nos quede
s¶olo la matriz identidad.
2
1 3 5 j 0
3
7
6
6 0 1 0 j 1 7
5
4
0 0 1 j ¡1
3
2
1 0 0 j 2
7
6
6
)4 0 1 0 j 1 7
5
0 0 1 j ¡1
es decir, x1 = 2; x2 = 1; x3 = ¡1
Veamoslo en la funci¶on de Matlab:
function x = gelim(A,b)
% Soluci¶
on para el sistema lineal de ecuaciones Ax = b
% usando la eliminaci¶
on de Gauss
52
% A es la matriz de coeficientes nxn
% b es el vector del lado derecho de nx1
[n,nl] = size(A);
for i=1:n-1
[pivot,k]= max(abs(A(i:n,i)));
% k es la posici¶
on del pivote, relativo a la fila i
if k > 1
temp1 = A(i,:);
temp2 = b(i,:);
A(i,:)
= A(i+k-1,:);
b(i,:)
= b(i+k-1,:);
A(i+k-1,:)
= temp1;
b(i+k-1,:)
= temp2;
end
for h =i+1:nl
m = A(h,i)/A(i,i);
A(h,:)
= A(h,:)-m*A(i,:);
b(h,:)
= b(h,:)-m*b(i,:);
end;
end;
% sustituci¶
on hacia atr¶
as
x(n,:)=b(n,:)./A(n,n);
for i= n-1:-1:1
x(i,:)
= (b(i ,:)-A(i,i+1:n)*x(i+1:n,:))./A(i,i);
end
As¶³ que en este caso escribimos gelim(A,b), con A; b de¯nidos como en el ejercicio.
4.1.2
Descomposici¶
on LU
Este es un m¶etodo que utiliza matrices triangulares como base para resolver matrices no
singulares en general. Para resolver Ax = b para A general, primero factorizamos A en el
producto de dos matrices triangulares, A = LU, en donde L es triangular baja (lower ) y
U es triangular alta (upper ). A esta factorizaci¶on se le llama la descomposici¶on LU de A:
53
Posteriormente reemplazamos el problema Ax = b con el problema equivalente LUx = b: En
la pr¶actica, es usual que L tenga en la diagonal 1's (m¶etodo de Doolittle), por lo tanto una
matriz de 3£3 queda de la siguiente manera:
2
1
0
0
3 2
u11 u12 u13
3
2
a11 a12 a13
3
7
7 6
6 a21 a22 a23 7
=
u22 u23 7
5
5 4
a31 a32 a33
0
u33
7 6
6
6 l21 1 0 7 ¢ 6 0
5 4
4
0
l31 l32 1
Comenzamos encontrando u11 = a11 y despu¶es resolviendo para los dem¶as elementos
de la primera ¯la de U y la primera columna de L. Como segundo paso, encontramos u22
y despu¶es el resto de la segunda ¯la de U y de la segunda columna de L.
determinamos todos los elementos de U y de L:
2
1 4
5
Siguiendo as¶³
3
7
6
7 escribimos el producto deseado
Si queremos la factorizaci¶on LU de A = 6
4
20
32
5
4
5 32 64
2
1
0
0
3 2
u11 u12 u13
7 6
6
6 l21 1 0 7 ¢ 6 0
5 4
4
0
l31 l32 1
3
2
1 4
3
5
7
7 6
6 4 20 32 7
=
u22 u23 7
5
5 4
5 32 64
0
u33
Empezamos resolviendo para la primera ¯la de U y la primera columna de L :
(1)u11 = a11 = 1
(1)u12 = a12 = 4
(1)u13 = a13 = 5
l21 u11 = a21 = 4
l31 u11 = a31 = 5
Usando estos valores, encontramos la segunda ¯la de U y la columna de L:
2
3
3 2
2
3
1 4 5
1 4
5
1 0 0
7
7 6
7 6
6
6 4 1 0 7 ¢ 6 0 u22 u23 7 = 6 4 20 32 7
5
5 4
5 4
4
5 32 64
0 0
u33
5 l32 1
54
(4)(4) + u22 = 20 ) u22 = 4
(4)(5) + u23 = 32 ) u23 = 12
(5)(4) + l32 u22 = 32 ) l32 = 3
Finalmente s¶olo nos queda una desconocida en U :
2
1 0 0
3 2
1 4 5
3
2
1 4
5
3
7
7 6
7 6
6
6 4 1 0 7 ¢ 6 0 4 12 7 = 6 4 20 32 7
5
5 4
5 4
4
5 32 64
0 0 u33
5 3 1
(5)(5) + (3)(12) + u33 = 64 ) u33 = 3
Por lo que la descomposici¶on queda:
2
1 0 0
3
2
1 4 5
3
7
6
7
U =6
0
4
12
5
4
0 0 3
7
6
7;
L=6
4
1
0
5
4
5 3 1
La funci¶on en Matlab es:
function [L,U] = LU(A)
[n, m] = size(A); % La dimensi¶
on de A
U = zeros(n, n); % Inicializar U
L = eye(n); % Inicializar L
for k = 1 :
n % computar el paso k
U(k, k) = A(k, k) - L(k, 1:k-1)*U(1:k-1, k);
for j = k+1 :
n
U(k, j) = A(k, j) - L(k, 1:k-1)*U(1:k-1, j);
L(j, k) = (A(j, k) - L(j, 1:k-1)*U(1:k-1, k))/U(k, k);
end
end
As¶³ que escribimos en la ventana de Matlab/Octave [L,U]=LU(A).
octave:10> A=[1 4 5;4 20 32;5 32 64]
A=
55
145
4 20 32
5 32 64
octave:15> [L,U]=LU(A)
L=
100
410
531
U=
145
0 4 12
003
4.2
4.2.1
M¶
etodos iterativos
M¶
etodo de Gauss-Jacobi
Los m¶etodos de descomposici¶on para las ecuaciones lineales son m¶etodos directos de soluci¶on.
Sin embargo, dichos m¶etodos pueden ser muy costosos para sistemas grandes. Para estos
casos existen los m¶etodos iterativos, los cuales son r¶apidos y buenos. Dichos m¶etodos son
valiosos, ya que las ideas b¶asicas se generalizan para los sitemas no lineales.
El m¶etodo de Gauss-Jacobi es uno de los m¶as sencillos para resolver ecuaciones lineales.
Comienza con la observaci¶on de que cada ecuaci¶on es una sola ecuaci¶on lineal, un tipo de
ecuaci¶on en la cual podemaos resolver una variable en t¶ermino de las otras.
Consideremos la ecuaci¶on de la primera ¯la de Ax = b:
a11 x1 + a12 x2 + ::: + a1n xn = b1
Podemos resolver para x1 en t¶erminos de (x2 ; :::; xn ) si a11 6= 0; resultando
¡1
(b1 ¡ a12 x2 ¡ ::: ¡ a1nxn )
x1 = a11
56
En general, si aii 6= 0; podemos usar la ¯la i de A para resolver para xi , encontrando
1
xi =
aii
Ã
bi ¡
X
aij xj
j6=i
!
Aparentemente este parece una manera diferente de expresar el sistema lineal, pero
vamos a convertirlo a continuaci¶on en un m¶etodo iterativo. Empecemos con un guess para
x y usemos las soluciones de las ecuaciones solas para computar un nuevo guess de x:
Suponiendo que aii 6= 0 podemos verlo como:
xik+1
1
=
aii
Ã
bi ¡
X
j6=i
aij xkj
!
;
para i = 1; :::; n
(4.1)
Este m¶etodo reduce el problema de resolver n inc¶ognitas simultanemente en n ecuaciones a resolver repetidamente n ecucaciones con una inc¶ognita. Esto s¶olo de¯ne una
sucesi¶on de guesses, por lo que esperamos que (4.1) converja a la verdadera soluci¶on.
Veamos un ejemplo gr¶a¯co de este m¶etodo. Primero consideremos un sistema de 2£2
2x + y = 6
x + 2y = 6
Con el m¶etodo de Gauss-Jacobi, reescribimos el sistema como:
1
x = ¡ y+3
2
1
y = ¡ x+3
2
Empezando con x1 =
1
2
y y1 =
1
2
,la primera ecuaci¶on produce el siguiente estimado
para x (usando y 1 ); y la segundo ecuaci¶on da el siguiente valor de y (usando x1 ):
1
1
x2 = ¡ y 1 + 3 = ¡ + 3 =
2
4
1
1
y 2 = ¡ x1 + 3 = ¡ + 3 =
2
4
11
4
11
4
Notemos que los nuevos valores de las variables no se usan hasta que empieza una
nueva iteraci¶on. A esto se le llama actualizaci¶on simult¶anea. La primera iteraci¶on se ilustra
en la siguiente gr¶a¯ca:
57
y
6
5
4
3
y(2)
2
1
y(1)
0
0.5
1
1.5
x
x(1)
2
2.5
x(2)
3
Gr¶a¯ca 4.1. La primera iteraci¶on del m¶etodo de Gauss-Jacobi
Veamos ahora en Matlab el siguiente sistema de 3£3
2x1 ¡ x2 + x3 = ¡1
x1 + 2x2 ¡ x3 = 6
x1 ¡ x2 ¡ 2x3 = ¡3
el cual se convierte en
x1 = :5x2 ¡ :5x3 ¡ :5
x2 = ¡:5x1 + :5x3 + 3
x3 = ¡:5x1 + :5x2 ¡ 1:5
En notaci¶on matricial, Ax = b; es
2
32
3
2
3
¡1
x1
2 ¡1 1
7
7 6
76
6
7
6
6
7
6 1 2
¡1 5 4 x2 5 = 4 6 7
5
4
¡3
x3
1 ¡1 2
58
que se transforma a
2
6
6
4
3
xk1
xk2
xk3
32
2
0
:5 ¡:5
76
7 6
7 = 6 ¡:5 0 :5 7 6
54
5 4
¡:5 :5 0
3
2
¡:5
3
x1k¡1
x2k¡1
x3k¡1
3
¡:5
7
7 6
7
7+6 3
5
5 4
¡1:5
Empezando con x0 = (0; 0; 0) encontramos
2
x11
3
2
0
:5 ¡:5
7 6
6
6 x1 7 = 6 ¡:5 0 :5
4 2 5 4
¡:5 :5 0
x13
32
0
3
2
2
¡:5
3
7
7 6
76 7 6
7
7=6 3
76 0 7 + 6 3
5
5 4
54 5 4
¡1:5
¡1:5
0
Para la segunda iteraci¶on tendr¶³amos
2
6
6
4
x21
x22
x23
3
3
32
2
2
3
2
3
1:75
¡:5
¡:5
0
:5 ¡:5
7
7 6
7 6
76
7 6
7 = 6 2:5 7
7+6 3
7 = 6 ¡:5 0 :5 7 6 3
5
5 4
5 4
54
5 4
:25
¡1:5
¡1:5
¡:5 :5 0
El m¶etodo de Gauss-Jacobi va a converger en 20 iteraciones al vector
x = [1:0002 2:0001
¡ 0:9997]0
usando " = 0:001; siendo la soluci¶on verdadera x¤ = [1 2 ¡1]0 :
La funci¶on en Matlab es:
function x=gjacobi(A,b,x0,tol,max)
% Soluci¶
on al sistema de ecuaciones lineales Ax = b
% usando el algoritmo iterativo de Gauss-Jacobi
% Inputs:
% A matriz de coeficientes (nxn)
% b lado derecho (nx1)
% x0 soluci¶
on inicial (nx1)
% tol Tolerancia permitida al computar el cero
% max N¶
umero m¶
aximo de iteraciones permitido
% Output:
% x Vector soluci¶
on (nx1)
[n m] = size (A);
59
xold = x0;
C =-A;
for i = 1 :
n
C(i ,i) = 0;
end
for i = 1 :n
C(i ,:)
=C(i,:)/A(i,i);
end
for i = 1 :n
d(i ,1) = b(i)/A(i ,i);
end
i = 1;
disp( ' i x1 x2 x3 ....');
while (i <= max)
xnew = C*xold + d;
if norm(xnew-xold) <=tol
x = xnew;
disp('El m¶
etodo de Gauss-Jacobi ha convergido');
return;
else
xold=xnew;
end
disp([i xnew']);
i=i+1;
end
disp('El m¶
etodo de Gauss-Jacobi NO ha convergido');
disp('Resulatdos despu¶
es del n¶
umero m¶
aximo de iteracione');
x=xnew;
Veamoslo para nuestro ejercicio, en donde escribimos gjacobi(A,b,[0;0;0],0.001,20)
octave:23> gjacobi(A,b,[0;0;0],0.001,20)
i x1 x2 x3 ....
1.00000 -0.50000 3.00000 -1.50000
2.00000 1.75000 2.50000 0.25000
60
3.00000 0.62500 2.25000 -1.12500
4.00000 1.18750 2.12500 -0.68750
5.00000 0.90625 2.06250 -1.03125
6.00000 1.04688 2.03125 -0.92188
7.00000 0.97656 2.01562 -1.00781
8.00000 1.01172 2.00781 -0.98047
9.00000 0.99414 2.00391 -1.00195
10.00000 1.00293 2.00195 -0.99512
11.00000 0.99854 2.00098 -1.00049
12.00000 1.00073 2.00049 -0.99878
13.00000 0.99963 2.00024 -1.00012
El m¶etodo de Gauss-Jacobi ha convergido
ans =
1.00018
2.00012
-0.99969
Generalizaci¶
on para problemas no lineales
Ahora generalicemos este m¶etodo para problemas no lineales. Dado el valor de la
iteraci¶on k, xk , usamos la ecuaci¶on i para computar el componente i de la inc¶ognita xk+1 ,
nuestra siguiente iteraci¶on. Formalemente, xk+1 se de¯ne en t¶erminos de xk por las ecuaciones:
f 1 = (x1k+1 ; xk2 ; xk3 ; :::; xkn ) = 0
(4.2)
f 2 = (xk1 ; x2k+1 ; xk3 ; :::; xkn ) = 0
..
.
f n = (xk1 ; xk2 ; xk3 ; :::; xkn¡1 ; xnk+1 ) = 0
Cada ecuaci¶on en (4.2) es una sola ecuaci¶on no lineal con una inc¶ognita, por lo que
tenemos las mimas ventajas que para el caso lineal.
El m¶etodo de Gauss-Jacobi se ve afectado por el esquema de indexaci¶on de las variables
y de las ecuaciones, por lo que no hay una elecci¶on natural para decidir qu¶e variable ser¶a
la primera y cu¶al ecuaci¶on la primera. Por lo tanto tenemos
n(n¡1)
2
esquemas diferentes de
Gauss-Jacobi, y es muy dif¶³cil determinar cual es el mejor esquema. Sin embargo, hay algunas
cosas que podemos hacer, por ejemplo, si alguna ecuaci¶on s¶olo depende de una inc¶ognita,
61
entonces esa ecuaci¶on deber¶³a ser la primera y esa variable la primera. En general, lo que
vamos a queres es que (4.2) se asemeje a una forma no lineal de sustituci¶on hacia atr¶as.
Cada paso en el m¶etodo de Gauss-Jacobi es una ecuaci¶on no lineal y se resuelve por lo
general con un m¶etodo iterativo. No tienen tanto sentido resolver cada ecuaci¶on exactamente,
dado que queremos resolver cada ecuaci¶on en la siguiente iteraci¶on. Por lo que podr¶³amos
simplemente resolver aproximando cada ecuaci¶on en (4.2), lo cual conlleva al m¶etodo lineal
Gauss-Jacobi, el cual utiliza un s¶olo paso del Newton para aproximar los componentes de
xk+1 ; resultando en el siguiente esquema:
xik+1
=
xki
f i (xk )
¡ i k ;
fxi (x )
i = 1; :::; n
Dadas las desventajas del m¶etodo de Gauss-Jacobi, pasemos a continuaci¶on al siguiente
m¶etodo, donde daremos un ejemplo.
4.2.2
Algor¶³tmo de Gauss-Seidel
Para el caso de un sistema de ecuaciones lineales, el m¶etodo de Gauss-Jacobi usamos un
nuevo guess para xki ; xik+1 ; s¶olo despu¶es de que hemos computado todo el vector de los nuevos
valores xk+1 : Este retraso para usar la nueva informaci¶on puede no ser sensible. Supongamos
que x¤ es la soluci¶on a Ax = b: Si, como esperamos, xik+1 es una mejor estimaci¶on de x¤i
k+1
que xki , entonces usar xik+1 para computar xi+1
en (4.1) parecer¶³a que es mejor que usar xki :
El argumento intuitivo es que la nueva informaci¶on sobre x¤i debe ser usarse lo m¶as pronto
posible.
La idea b¶asica del m¶etodo de Gauss-Seidel es usar una nueva aproximaci¶on de x¤i en
cuanto est¶e disponible. Es decir, si nuestro nuevo guess para x¤ es xk ; entonces compuatmos
nuestro nuevo guess para x1 como en (4.1):
¡
¢
k
k
x1k+1 = a¡1
b
¡
a
x
¡
:::
¡
a
x
)
1
12
1n
11
2
n
pero usamos este nuevo guess para x¤1 inmediantamente cuando computemos los otros nuevos
componentes de xk+1 : En particular, nuestra ecuaci¶on para x2k+1 derivada de la segunda ¯la
de A se volver¶³a
¡
¢
k+1
x2k+1 = a¡1
¡ a23 xk3 ¡ ::: ¡ a2n xkn )
22 b2 ¡ a21 x1
62
En general, para resolver Ax = b con el m¶etodo de Gauss-Seidel, de¯nimos la sucesi¶on
© k ª1
x k=1 por la iteraci¶on
xik+1
1
=
aii
Ã
bi ¡
i¡1
X
j=1
aij xjk+1 ¡
n
X
aij xkj
j=i+1
!
;
para i = 1; :::; n
(4.3)
Podemos observar claramente c¶omo los nuevos componentes de xk+1 son usado inmediatamente de que son computados, dado que xik+1 es usado para computar xjk+1 para j > i:
Not¶emos que en este m¶etodo el orden en el cual resolvemos para los componentes
sucesivos de x importa, mientras que el orden no importa en el m¶etodo de Gauss-Jacobi.
Por lo que esto le da al m¶etodo de Gauss-Seidel m¶as °exibilidad, en particular si un orden
no converge, entonces tratamos otro orden.
Veamos el mismo ejemplo de 2£2 que vimos en la secci¶on anterior, es decir,
2x + y = 6
x + 2y = 6
Con el m¶etodo de Gauss-Seidel, reescribimos el sistema como:
1
x = ¡ y+3
2
1
y = ¡ x+3
2
Empezando con x1 =
1
2
y y1 =
1
2
,la primera ecuaci¶on produce el siguiente estimado
1
para x (usando y ); y la segundo ecuaci¶on da el siguiente valor de y (usando x2 ):
1
1
11
x2 = ¡ y 1 + 3 = ¡ + 3 =
2
4
4
11
13
1
y 2 = ¡ x2 + 3 = ¡ + 3 =
2
8
8
La primera iteraci¶on se ilustra en la siguiente gr¶a¯ca:
63
y 6
5
4
3
2
y(2)
1
y(1)
0
0.5
1
1.5
x
2
2.5
x(2)
3
Gr¶a¯ca 4.2. La primera iteraci¶on del m¶etodo de Gauss-Seidel
Veamos la funci¶on en Matlab de Gauss-Seidel usando el mismo ejercicio que para
Gauss-Jacobi:
function x=gseidel(A,b,x0,tol,max)
% Soluci¶
on al sistema de ecuaciones lineales Ax = b
% usando el algoritmo iterativo de Gauss-Seidel
% Inputs:
% A matriz de coeficientes (nxn)
% b lado derecho (nx1)
% x0 soluci¶
on inicial (nx1)
% tol Tolerancia permitida al computar el cero
% max N¶
umero m¶
aximo de iteraciones permitido
% Output:
% x Vector soluci¶
on (nx1)
[n m] = size (A);
x = x0;
C=-A;
for i = 1 :
n
C(i ,i) = 0;
end
for i = 1 :n
64
C(i,1:n) = C(i,1:n)/A(i,i);
end
for i = 1:n
r(i,1) = b(i)/A(i,i);
end
i = 1;
disp( ' i x1 x2 x3 ....');
while (i <= max)
xold =x;
for j = 1:n
x(j) = C(j, :)*x + r(j);
end
if norm(xold - x) <= tol
disp('El m¶
etodo de Gauss-Seidel ha convergido');
return;
end
disp([i x']);
i=i+1;
end
disp('El m¶
etodo de Gauss-Seidel NO ha convergido');
disp('Resulatdos despu¶
es del n¶
umero m¶
aximo de iteracione');
x=xnew;
En la ventana de Matlab escribimos gseidel(A,b,[0;0;0],0.001,10) y observamos
que nos sale igual que para el caso de Gauss-Jacobi.
octave:31> gseidel(A,b,[0;0;0],0.001,10)
i x1 x2 x3 ....
1.00000 -0.50000 3.25000 0.37500
2.00000 0.93750 2.71875 -0.60938
3.0000 1.1641 2.1133 -1.0254
4.0000 1.0693 1.9526 -1.0583
5.0000 1.0055 1.9681 -1.0187
6.00000 0.99339 1.99395 -0.99972
7.00000 0.99684 2.00172 -0.99756
65
8.00000 0.99964 2.00140 -0.99912
9.00000 1.00026 2.00031 -0.99997
El m¶etodo de Gauss-Seidel ha convergido
ans =
1.0001
1.9999
-1.0001
Generalizaci¶
on para problemas no lineales
En el m¶etodo de Gauss-Jacobi para ecuaciones no lineales usamos otra vez el nuevo
guess de xi ; xik+1 , s¶olo despu¶es de haber computado todo el vector de nuevos valores xk+1 :
La idea b¶asica del m¶etodo de Gauss-Seidel es de nuevo, usar el nuevo guess de xi en cuanto
est¶e disponible. En el caso no lineal, esto implica que dado xk , construimos xk+1 componente
a componente resolviendo los problemas unidimensionales en sucesi¶on:
f 1 = (x1k+1 ; xk2 ; xk3 ; :::; xkn ) = 0
f 2 = (x1k+1 ; x2k+1 ; xk3 ; :::; xkn ) = 0
..
.
k+1
; xkn ) = 0
f n¡1 = (x1k+1 ; x2k+1 ; :::; xn¡1
k+1
f n = (x1k+1 ; x2k+1 ; :::; xn¡1
; xnk+1 ) = 0
Resolvemos f 1 ; f 2 ; :::; f n en sucesi¶on, usando inmediatamente cada nuevo componente.
Notemos que ahora la notaci¶on de ¶³ndices importa, ya que afecta la manera en que los
u
¶ltimos resultados dependen de los primeros. Cabe se~
nalar que tambi¶en podemos dar una
forma lineal del m¶etodo de Gauss-Seidel
µ i¶
f
k+1
k
k+1
xi = xi ¡
(x1k+1 ; :::; xi¡1
; xki ; :::; xkn );
fxi i
donde i = 1; 2; :::; n, con el ¯n de optimizar en los costos de optimizaci¶on de cada iteraci¶on.
Estos m¶etodos Gaussianos son usados con frecuencia, pero tienen algunos problemas.
Como en el caso lineal, son m¶etodos riesgosos para usar si el sistema no es diagonalmente
dominante, es decir, si cada elemento de la diagonal es mayor a la suma de las magnitudes
66
de los elementos no diagonales en su ¯la, matem¶aticamente hablando,
X
j6=i
jaij j < jaii j
i = 1; :::; n
Teorema 4.1 Si A es diagonalmente dominante, tanto el esquema iterativo de GaussJacobi como el de Gauss-Seidel convergen para todos los valores iniciales.
Aplicaci¶
on para sistemas de ecuaciones en diferencia de segundo orden
Recordemos un poco de los sistemas de ecuaciones en diferencia de segundo orden no
lineal, en donde tenemos una ecuaci¶on de la forma g(xt ; xt+1 ; xt+2 ) con una condici¶on inicial,
es decir, x0 dado y una condici¶on terminal limt!1 xt = x¤ : Si suponemos que en un tiempo
¯nito T se llega a x¤ ; entonces vamos a tener el siguiente sistema de ecuaciones no lineal:
g(x0 ; x1 ; x2 ) = 0
g(x1 ; x2 ; x3 ) = 0
..
.
g(xt¡2 ; xt¡1 ; x¤ ) = 0
que ya viene siendo un sistema con T ¡ 2 ecuaciones y T ¡ 2 inc¶ongitas, con lo cual
podemos ahora si aplicar nuestros metodos num¶ericos.
Dado que en econom¶³a nosotros veremos funciones que depender¶an de varias variables,
como el capital, el consumo, etc., y en donde el sistema de funciones ser¶a no lineal, veamos
a continuaci¶on un ejemplo sencillo del modelo de crecimiento neocl¶asico con plani¯cador
social.
El plan¯cador central va a querer maximizar lo siguiente:
max
1
X
¯ t u(ct )
(4.4)
t=0
s:a:
ct + it = f(kt )
kt+1 = (1 ¡ ±)kt + it
8t
8t
La funci¶on en Matlab quedar¶³a as¶³:
% gseide.m Programa para resolver una simple versi¶
on del problema del planificado
social
% usando Gauss-Seidel.
El problema a resolver es:
67
%
%
inf
% max sum beta^t ln(ct)
%
t=0
%
% s.a.
ct+it = A*k^alpha
% kt+1 = (1-delta)*kt+it
%
% con k0 dado
%
%
clear
% Par¶
ametros del modelo
alpha = 0.35;
beta = 0.95;
delta = 0.06;
A = 10;
% Otros par¶
ametros (m¶
aximo de iteraciones, tolerancia deseada y
% n¶
umero de periodos antes de alcanzar el estado estacionario)
maxit = 250;
crit = 1e-2;
T = 100;
% Computa el capital de estado estacionario (kss), el capital inicial (k0)
y
% el guess inicial para la sucesi¶
on de kt (kg)
kss = ((A*beta*alpha)/(1-(1-delta)*beta))^(1/(1-alpha));
k0 = (1/2)*kss;
kg = zeros(T+1,1);
ksol = zeros(T,1);
kg(1) = k0;
ksol(1) = k0;
for t=2:T+1
kg(t) = kg(t-1)+(kss-k0)/T;
end
68
% Iteraciones de Gauss-Seidel
for i=1:maxit
for t=2:T
param = [alpha;beta;delta;A;ksol(t-1);kg(t+1)];
ksol(t) = secant sys('fgs',kss/2,param,1e-5,20);
end
dist = norm(kss-ksol(T));
plot(1:T,ksol(1:T),1:T,kss*ones(T,1))
title(['Trayectoria para el capital en la iteraci¶
on',int2str(i),...
' con error ',num2str(dist)])
pause(0.01)
if dist<crit
break
end
kg(1:T) = ksol(1:T);
end
% Gr¶
aficas de las trayectorias o
¶ptimas para el producto (yt),
% la inversi¶
on (it) y el consumo (ct)
kt = kg(1:T+1);
figure(2)
yt = A*kt(1:T).^alpha;
plot(1:T,yt)
title('Trayectoria o
¶ptima para el producto')
figure(3)
it = kt(2:T+1)-(1-delta)*kt(1:T);
plot(1:T,it)
title('Trayectoria o
¶ptima para la inversi¶
on')
figure(4)
ct = yt-it;
plot(1:T,ct)
title('Trayectoria o
¶ptima para el consumo')
Y la funci¶on de producci¶on ncesesaria para correr este programa la llamaremos fgs.m
que viene dada por:
function y = fgs(kt1,p)
69
alpha = p(1);
beta = p(2);
delta = p(3);
A = p(4);
kt = p(5);
kt2 = p(6);
y = A*kt1^alpha+(1-delta)*kt1-beta*(A*kt^alpha-kt1+ ...
(1-delta)*kt)*(A*alpha*kt1^(alpha-1)+(1-delta))-kt2;
Notemos que estamos incluyendo todos los par¶ametros en nuestro archivo gseide.m,
por lo que en la ventana de Matlab s¶olo escribimos gseide.m. Los resultado del capital con
las 250 iteraciones dadas ser¶³a que kt = 197:7 y las gr¶a¯cas quedar¶³an:
Trayectoria óptima para el consumo
Trayectoria para el capital en la iteración250 con error 0.01755
52
200
50
180
48
46
160
44
140
42
40
120
38
100
36
80
0
10
20
30
40
50
60
70
80
90
100
34
0
20
30
40
50
60
70
80
90
100
Trayectoria óptima para el producto
Trayectoria óptima para la inversión
64
15.5
62
15
14.5
60
14
58
13.5
56
13
54
12.5
52
12
50
11.5
10
0
10
20
30
40
50
60
70
80
90
100
48
0
70
10
20
30
40
50
60
70
80
90
100
5
Programaci¶
on din¶
amica num¶
erica
5.1
Funci¶
on de valor
Esta sesi¶on tiene como objetivo presentar el m¶etodo de iteraci¶on de la funci¶on de valor.
Para ello, empezamos reformulando la de¯niciones de equilibrio en t¶erminos recursivos y
presentando algunos resultados de la literatura sobre programaci¶on din¶amica. Estos resultados constituyen la base te¶orica para el m¶etodo de iteraci¶on de la funci¶on de valor, que se
presenta a continuaci¶on para resolver el problema del plani¯cador social en un contexto determin¶³stico. Posteriormente se ver¶a la extensi¶on para resolver y simular un modelo estoc¶astico
con shocks tecnol¶ogicos.
Recordando la ecuaci¶on (4.4) y resolviendo el problema del plani¯cador central, nos
vamos a encontrar con la siguiente ecuaci¶on de Euler
¯f 0 (kt )u0 [f (kt ) ¡ kt+1 ] = u0 [f (kt¡1 ) ¡ kt ];
t = 1; 2; :::; T
(5.1)
y las condiciones iniciales y ¯nales
kT +1 = 0;
k0 > 0 dado
La ecuaci¶on (5.1) es una ecuaci¶on en diferencia de segundo orden en kt , por lo tanto
su soluci¶on tiene una familia de soluciones de dos par¶ametros. Sin embargo, el u
¶ nico ¶optimo
para el problema de maximizaci¶on es aquella soluci¶on en esta familia que adem¶as satisface
las condiciones iniciales y terminales.5 Una manera de resolver el problema es que dado que
supusimos que el tiempo es ¯nito, resolvemos el problema y la trayectoria de capital ser¶a
kt+1 = ®¯
1 ¡ (®¯)T ¡1 ®
k ;
1 ¡ (®¯)T ¡t+1 t
t = 0; 1; :::; T
si lo queremos convertir en un problema de horizonte in¯nito, s¶olo tomamos l¶³mites y nos
queda
kt+1 = ®¯kt® ;
t = 0; 1; :::
(5.2)
De hecho esta conjetura va a ser la correcta y esta va a ser la u
¶nica soluci¶on del
5
Cap¶itulo 2 de SLP
71
problema. A la ecuaci¶on (5.2) se le llama funci¶on de pol¶³tica y se puede escribir de tambi¶en
como kt+1 = g(kt ) para t = 0; 1; :::, en donde g : R+ ! R+ .
Aunque empezamos nuestro problema escogiendo sucesiones in¯nitas para el consumo
y el capital f(ct ; kt+1 )g1
t=0 ; el problema que en realidad enfrenta el plani¯cador central en el
periodo t = 0 es el de escoger s¶olo el consumo actual, c0 ; y el capital de ma~
nana k1 , el resto
puede esperar hasta ma~
nana. As¶³ que si conocemos las preferencias sobre estos dos bienes,
entonces podemos realizar nuestro sencillo problema de maximizar (c0 ; k1 ) sujerto a nuestra
restricci¶on presupuestal.
Dado k0 podemos de¯nir una funci¶on v : R+ ! R tomando v(k0 ) como el valor de la
funci¶on objetivo maximizada para cada k0 ¸ 0: Una funci¶on de este tipo se le llama funci¶on
de valor. Con v de¯nida de esta manera, v(k1 ) dar¶³a el valor de de la utilidad del periodo
1 usando k1 y ¯v(k1 ) ser¶³a el valor de esta utilidad descontada al periodo 0. Entonces el
problema del planeador central en el periodo 0 usando esta funci¶on de valor ser¶³a
max[u(c0 ) + ¯v(k1 )]
c0 ;k1
s:a: c0 + k1 · f (k0 )
(5.3)
c0 ; k1 ¸ 0
k0 > 0 dado
Si (5.3) provee una soluci¶on para el problema, entonces v(k0 ) debe de ser la funci¶on
objetivo maximizadora de (5.3) tambi¶en. Por lo tanto, v debe de satifacer
v(k0 ) =
max
0·;k1 ·f (k0 )
fu[f (k0 ) ¡ k1 ] + ¯v(k1 )g
Notemos que cuando vemos el problema de una manera recursiva, los sub¶³ndices de
tiempo no importan, por lo que podemos reescribir el problema como
v(k) =
max fu[f (k) ¡ k 0 ] + ¯v(k 0 )g
0·;k'·f (k)
Esta ecuaci¶on en la funci¶on desconocida v se llama ecuaci¶on funcional. El estudio de
los problemas de optimizaci¶on din¶amica a trav¶es del an¶alisis de estas ecuaciones funcionales
se llama programaci¶on din¶amica.
72
5.2
Formulaci¶
on Recursiva y Programaci¶
on Din¶
amica6
Un Equilibrio General Competitivo Recursivo es un conjunto de funciones v (k; K), c (k; K),
i (k; K), k 0 (k; K), precios w (K) y r (K) y ley de movimiento ¡ (K) tales que:
i) Dadas las funciones w, r y ¡, la funci¶on de valor v (k; K) resuelve la ecuaci¶on de
Bellman:
v (k; K) = max0
c;i;k
fu (c) + ¯v (k 0 ; K 0 )g
(5.4)
s:t: c + i = w (K) + r (K) k
k 0 = (1 ¡ ±) k + i
K 0 = ¡ (K)
y las funciones c (k; K), i (k; K) y k 0 (k; K) son reglas de decisi¶on ¶optimas para este problema.
ii) Para todo K, los precios satisfacen las condiciones marginales:
r (K) = f 0 (K)
(5.5)
w (K) = f (K) ¡ f 0 (K) K
(5.6)
iii) Para todo K, hay igualdad entre oferta y demanda:
f (K) = c (K; K) + i (K; K)
(5.7)
iv) Para todo K, ley de movimiento agregada y comportamiento individual son consistentes:
¡ (K) = k 0 (K; K)
6
(5.8)
Esta secci¶
on fue tomada de las notas de Carlos Urrutia "M¶
etodos Num¶
ericos para Resolver
Modelos Macroecon¶
omicos Din¶
amicos".
73
Resolver un Equilibrio Competitivo en su formulaci¶on recursiva consiste entonces en encontrar un conjunto de funciones que satisfagan las condiciones anteriores. En particular,
una vez encontradas las reglas de decisi¶on ¶optimas y partiendo de un k0 dado, podemos
reconstruir las secuencias para las principales variables de manera recursiva:
k1 = k 0 (k0 ; k0 )
k2 = k 0 (k1 ; k1 ) = k 0 (k 0 (k0 ; k0 ) ; k 0 (k0 ; k0 ))
y as¶³ sucesivamente para kt , t = 3; 4; ::: El Principio de Optimalidad asegura que la secuencia
resultante es igual a la que obtendr¶³amos resolviendo el equilibrio en forma de secuencias
que vimos en la sesi¶on anterior.
De manera similar, de¯nimos un Optimo de Pareto Recursivo como un conjunto de
funciones v (k), c (k), i (k), k 0 (k) que resuelven la ecuaci¶on de Bellman del plani¯cador
social:
v (k) = max0
c;i;k
fu (c) + ¯v (k 0 )g
(5.9)
s:t: c + i = f (k)
k 0 = (1 ¡ ±) k + i
La equivalencia entre el equilibrio competitivo y el problema del plani¯cador social se sigue
manteniendo en este contexto.
5.2.1
Programaci¶
on Din¶
amica
Las ecuaciones de Bellman (5.4) y (5.9) son ejemplos de un problema m¶as general. Sea
x un vector columna de n componentes, X un subconjunto de Rn , la funci¶on de retorno
F : X £ X ! R y la correspondencia - : X ! X. Con esos elementos, de¯nimos la funci¶on
de valor v : X ! R como la soluci¶on de la ecuaci¶on de Bellman:
v (x) = max
y
fF (x; y) + ¯v (y)g
74
(5.10)
s:t: y 2 - (x)
para todo x 2 X. Dada la funci¶on de valor, podemos de¯nir la regla de decisi¶on ¶optima
g : X ! X como:
g (x) = arg max
y
fF (x; y) + ¯v (y)g
(5.11)
s:t: y 2 - (x)
es decir, la funci¶on que satisface:
v (x) = F (x; g (x)) + ¯v (g (x))
Las propiedades de las funciones v y g han sido analizadas en la literatura sobre programaci¶on din¶amica7 . En ella, se demuestra que si (i) X es un conjunto convexo; (ii) - (x)
es un conjunto compacto y no-vac¶³o, para todo x 2 X; (iii) la correspondencia - es convexa
y continua; (iv) la funci¶on F es c¶oncava, acotada y continua; y (v) ¯ < 1; entonces:
1. existe una u
¶ nica funci¶on v que satisface (5.10);
2. v es continua, acotada y estrictamente c¶oncava;
3. existe una u
¶ nica funci¶on g que resuelve (5.11);
4. g es continua.
En la mayor parte de ejemplos pr¶acticos que veremos, incluyendo el Modelo de Crecimiento Neocl¶asico, las restricciones impuestas en las preferencias y tecnolog¶³a aseguran que
los supuestos (i)-(v) se cumplen.
Aparte de asegurarnos existencia y unicidad de la soluci¶on, el m¶etodo de programaci¶on
din¶amica nos ofrece un procedemiento para encontrar dicha soluci¶on. De¯namos el operador
T : B(X) ! B (X), en donde B (X) es el espacio de funciones acotadas de¯nidas en X, de
la siguiente manera:
7
Una referencia completa, aunque bastante t¶ecnica, puede encontrarse en Stokey y Lucas (1989).
75
T [v (x)] = max
y
fF (x; y) + ¯v (y)g
(5.12)
s:t: y 2 - (x)
Encontrar la funci¶on de valor que resuelve (5.10) es equivalente a encontrar un punto ¯jo del
operador T , es decir, una funci¶on v tal que T [v] = v.
Bajo ciertas condiones t¶ecnicas, el operador T es una contracci¶on 8 . Por lo tanto,
partiendo de cualquier funci¶on v 0 2 B (X), la secuencia v n de¯nida por:
v1 = T v0
v2 = T v1 = T 2v 0
:::::::::::::
converge hacia el punto ¯jo v.
5.3
Iteraci¶
on de la Funci¶
on de Valor
Volvamos ahora al problema del plani¯cador social (5.9) en su froma recursiva. Reemplazando las restricciones en la funci¶on objetivo, tenemos:
v (k) = max
0
k
fu [f (k) + (1 ¡ ±) k ¡ k 0 ] + ¯v (k 0 )g
(5.13)
s:t: k 0 2 [0; f (k) + (1 ¡ ±) k]
8
El operador T es una contracci¶on con m¶
odulo ¯ si existe un n¶
umero ¯ 2 (0; 1) tal que kT f (x) ¡ T g (x)k ·
¯ kf (x) ¡ g (x)k para todo f; g 2 B (X), x 2 X.
Las llamadas condiciones de Blackwell su¯cientes para una contracci¶
on son: monotonicidad del operador
T , en el sentido que T f · T g si f (x) · g (x) para todo x 2 X, y descuento, en el sentido que existe un
n¶
umero ¯ 2 (0; 1) tal que T [f + a] (x) · T f (x) + ¯a para todo f 2 B (X), x 2 X y a ¸ 0. Nuevamente, los
supuestos (i)-(v) mencionados anteriormente garantizan su cumplimiento.
76
un problema equivalente a la ecuaci¶on de Bellman general (5.10) con x = k, y = k 0 , F (x; y) =
u [f (x) + (1 ¡ ±) x ¡ y] y - (x) = [0; f (x) + (1 ¡ ±) x].
De acuerdo con la teor¶³a de la programaci¶on din¶amica que revisamos brevemente,
partiendo de cualquier funci¶on v 0 (por ejemplo, v 0 (k) = 0) la secuencia vn de¯nida por
vn+1 (k) = max
0
k
fu [f (k) + (1 ¡ ±) k ¡ k 0 ] + ¯v n (k 0 )g
(5.14)
s:t: k 0 2 [0; f (k) + (1 ¡ ±) k]
converge a la soluci¶on v cuando n tiende a in¯nito. Esto constituye de por s¶³ un m¶etodo iterativo en un espacio de funciones; en esta secci¶on veremos como implementarlo num¶ericamente.
5.3.1
Implementaci¶
on Num¶
erica
Empezamos de¯niendo una malla de valores para k, es decir, un vector K = (K1; K2; :::; Kp )
con K1 = Kmin y Kp = Kmax . Por simplicidad, podemos usar puntos igualmente espaciados,
con K2 = Kmin +´, K3 = Kmin +2´ y as¶³ sucesivamente, en donde ´ = (Kmax ¡ Kmin ) = (p ¡ 1)
Mientras m¶as puntos usemos, es decir mientras p sea mayor y ´ m¶as peque~
no, la aproximaci¶on
ser¶a m¶as precisa.
Luego, de¯nimos una matriz M de la siguiente manera:
2
F (K1 ; K1 ) F (K1 ; K2 ) ::::: F (K1 ; Kp )
3
7
6
6 F (K2 ; K1 ) F (K2 ; K2 ) ::::: F (K2 ; Kp ) 7
7
M =6
6 :::::::::::::::: :::::::::::::::: ::::: :::::::::::::::: 7
5
4
F (Kp ; K1 ) F (Kp ; K2 ) ::::: F (Kp ; Kp )
en donde F (x; y) = u [f (x) + (1 ¡ ±) x ¡ y]. La matriz M contiene la funci¶on de retorno
evaluada para cada posible k 2 K y k 0 2 K. Sin embargo, sabemos que no todos esos valores
son posibles. La restricci¶on k 0 2 [0; f (k) + (1 ¡ ±) k] nos permite eliminar algunas celdas
de la matriz que son inalcanzables para el plani¯cador social, imponi¶endole un valor cero9 .
Hacemos entonces:
9
O un valor negativo (como por ejemplo ¡1000). El punto es impedir que el programa escoja esa celda
al calcular la maximizaci¶on.
77
Mij ´ F (Ki ; Kj ) = 0
si Kj > f (Ki ) + (1 ¡ ±) Ki
Una vez constru¶³da la matriz M , el resto del procedimiento es una simple iteraci¶on.
Empezamos con cualquier vector columna V 0 2 Rp (por ejemplo, V 0 = 0) e iniciamos el
siguiente algoritmo con s = 0:
¡
¢
1. Dado el vector V s = V1s ; V2s ; :::; Vps y la matriz M, calcular V s+1 de la siguiente
manera:
V
s+1
2
M11 + ¯V1s
M12 + ¯V2s
::: M1p + ¯Vps
6
6 M21 + ¯V1s M22 + ¯V2s ::: M2p + ¯Vps
= max 6
6 :::::::::::::::::: :::::::::::::::::: ::: ::::::::::::::::::
4
Mp1 + ¯V1s Mp2 + ¯V2s ::: Mpp + ¯Vps
3
7
7
7
7
5
donde el m¶aximo es calculado por ¯la, es decir:
V
s+1
2
ª
©
max M11 + ¯V1s ; M12 + ¯V2s ; :::; M1p + ¯Vps
©
ª
6
6 max M21 + ¯V1s ; M22 + ¯V2s ; :::; M2p + ¯Vps
= 6
6 ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
4
©
ª
max Mp1 + ¯V1s ; Mp2 + ¯V2s ; :::; Mpp + ¯Vps
3
7
7
7
7
5
2. Evaluar kV s+1 ¡ V s k. Si la norma es mayor que el criterio de tolerancia permitido,
regresar al paso 1 con s = s + 1. En caso contrario, el procedimiento termina con la
soluci¶on V = V s+1 .
Una vez obtenida la aproximaci¶on a la funci¶on de valor V (evaluada en cada uno de los
p puntos de la malla), la correspondiente regla de decisi¶on ¶optima G se obtiene calculando:
2
M11 + ¯V1
M12 + ¯V2
::: M1p + ¯Vp
3
7
6
6 M21 + ¯V1 M22 + ¯V2 ::: M2p + ¯Vp 7
7
G = arg max 6
6 :::::::::::::::::: :::::::::::::::::: ::: :::::::::::::::::: 7
5
4
Mp1 + ¯V1 Mp2 + ¯V2 ::: Mpp + ¯Vp
Es decir, G es un vector columna de n componentes, en donde Gi 2 f1; :::; pg nos indica el
n¶
umero de la columna que maximiza la ¯la i.
78
Podemos entonces reconstruir la secuencia ¶optima de capital, partiendo de k0 = Ki de
la siguiente manera:
k1 = Kj
con j = Gi
k2 = Kl
con l = Gj
:::::::::::::::::::
y as¶³ sucesivamente.
5.3.2
Resolviendo el Equilibrio Competitivo
El m¶etodo de iteraci¶on de la funci¶on de valor es particularmente apropiado para resolver el
problema del plani¯cador social. Sin embargo, su utilidad es menor al momento de resolver
directamente el Equilibrio Competitivo. La raz¶on no es que existan dos variables de estado,
el capital individual y el agregado, pues eso ser¶³a una simple extensi¶on del m¶etodo anterior.
El principal problema es que la ley de movimiento del capital agregado (la funci¶on ¡ en la
ecuaci¶on de Bellman (5.4)) debe ser conocida al momento de realizar la optimizaci¶on. Sin
embargo, este es un objeto que solo conocemos luego de resolver el equilibrio.
Una forma de evitar este problema guarda cierta analog¶³a con el metodo de doble
iteraci¶on que vimos en la sesi¶on anterior. Presentaremos solo la idea general, pues su utilizaci¶on es limitada en la pr¶actica. Suponemos que la ley de movimiento tiene cierta forma
funcional, por ejemplo un polinomio de grado n:
K 0 = ¡ (K) = a1 + a2 K 2 + ::: + an K n
y proponemos un vector inicial de par¶ametros (a1 ; a2 ; :::; an ).
Dada esa ley de movimiento, la ecuaci¶on (5.4) puede resolverse iterando la funci¶on de
valor. Obtenemos entonces la regla de decisi¶on ¶optima y, dado k0 , calculamos la secuencia
¶optima resultante k0 ; k1; :::; kT . Con esos datos, corremos la regresi¶on:
kt+1 = a1 + a2 kt2 + ::: + an ktn
y estimamos un nuevo juego de par¶ametros (^a1 ; a
^2 ; :::; a
^n ). Si estos estan razonablemente
79
cercanos a los que propusimos, termina el procedimiento. En caso contrario, volvemos a
resolver la ecuaci¶on de Bellman usando (^
a1 ; a
^2 ; :::; a
^n ) como los nuevos par¶ametros de la ley
de movimiento.
Evidentemente, este m¶etodo de doble iteraci¶on ser¶a m¶as preciso conforme mayor sea el
grado n del polinomio que se use para aproximar la ley de movimiento del capital agregado.
Pero a¶
un con un n grande, nada garantiza la convergencia del algoritmo.
Veamos ahora un ejemplo num¶erico para entender mejor la iteraci¶on de la funci¶on de
valor:
Suponer que queremos encontrar y gra¯car las trayectorias ¶optimas para el capital,
consumo, producto, salario real y tasa de inter¶es en el Modelo de Crecimiento Neocl¶asico
Simple, con las siguientes formas funcionales:
f (k) = Ak ®
u(c) = log(c)
los siguientes valores para los par¶ametros:
¯ = 0:95
A = 10
® = 0:35
± = 0:06
y un stock de capital inicial k0 igual a dos tercios del capital de estado estacionario k ¤ ,
utilizando el m¶etodo de iteraci¶on de la funci¶on de valor. Usar una malla para el stock de
capital de 500 puntos igualmente espaciados desde 0:5k ¤ hasta 1:2k ¤ y gra¯car las funciones
de valor obtenidas en cada una de las iteraciones.
Para esto vamos a escribir la funci¶on vfiterm.m como:
% vfiterm.m Programa para resolver una versi¶
on simple del problema
del planeador social
% usando la iteraci¶
on de la funci¶
on de valor.
La ecuaci¶
on de
Bellman a resolver es:
%
% v(k) = max fln(A*k^alpha+(1-delta)*k-k') + beta v(k')g
% k'
%
80
% con k0 dado
%
% Output:
file vfout.mat con la funci¶
on de valor, regla de decisi¶
on o
¶ptima,
% y sucesiones de las variables principales.
clear
% Parametros del modelo
alpha = 0.35;
beta = 0.95;
delta = 0.06;
A = 10;
% Otros par¶
ametros (m¶
aximo n¶
umero de iteraciones, tama~
no del grid,
% criterio de tolerancia, n¶
umero de periodos para la simulaci¶
on)
maxit = 1000;
p = 500;
crit = 1e-2;
T = 100;
% Computa el capital de estado estacionario (kss), el capital inicial (k0)
% y define el grid para k
kss = ((A*beta*alpha)/(1-(1-delta)*beta))^(1/(1-alpha));
k0 = (1/2)*kss;;
Kmin = 0.5*kss;
Kmax = 1.2*kss;
K(1)=k0;
K = [Kmin+(0:p-2)*(Kmax-Kmin)/(p-1) Kmax]';
% Ponemos el guess inicial para la funci¶
on de valor (V0=0) y construye la
matriz M
V0 = 0*K;
e = ones(1,p);
M = log(max((A*K.^alpha+(1-delta)*K)*e-e'*K',1e-8));
% Iteraci¶
on de la funci¶
on de valor
figure(1)
iconv = 0;
it = 1;
81
while (iconv==0 & it<maxit)
[V1,G] = max((M+beta*(V0*e)')');
V1 = V1';
G = G';
if norm(V0-V1)<crit
iconv = 1;
end
plot(K,V1)
axis([Kmin Kmax 0 100])
title(['Funci¶
on de valor en la iteraci¶
on ',int2str(it),...
' con error ',num2str(norm(V0-V1))]),...
xlabel('k'),...
ylabel('v(k)'),...
grid
pause(0.01)
V0 = V1;
it = it+1;
end
figure(2)
plot(K,K(G),K,K,':')
title('Regla de decisi¶
on o
¶ptima'),...
xlabel('k'),...
ylabel('kp=g(k)'),...
grid,...
text(kss,kss,'o kss')
%Veamos los estados estacionarios
yss=A*kss^alpha
css=yss+[(1-delta)-1]*kss
iss=yss-css
sss=iss/yss %la tasa de ahorro
% Simula trayectorias o
¶ptimas para el capital (kt), consumo (ct),
% inversi¶
on (it) y producto (yt)
figure(3)
kt = zeros(T+1,1);
82
ind = zeros(T,1);
[aux,ind(1)] = min(K<k0);
kt(1) = K(ind(1));
for t=1:T
ind(t+1) = G(ind(t));
kt(t+1) = K(ind(t+1));
end
kss2=kss*ones(T,1)';
z=[1,T,K(1),kss+1];
plot(1:1:T,kt(1:T),'.-',1:1:T,kss2,'r')
title('Trayectorias o
¶ptimas para el Capital'),...
xlabel('t'),...
ylabel('k t'),...
grid,...
axis(z)
legend('Senda', 'SS')
yss2=yss*ones(T-1,1)';
figure(4)
yt = A*kt(1:T).^alpha;
plot(1:1:T,yt,'.-',1:1:T-1,yss2,'r')
title('Trayectorias o
¶ptimas para el Producto'),...
xlabel('t'),...
ylabel('y t'),...
grid,...
legend('Senda', 'SS')
iss2=iss*ones(T-1,1)';
figure(5)
it = kt(2:T+1)-(1-delta)*kt(1:T);
plot(1:1:T,it,'.-',1:1:T-1,iss2,'r')
title('Trayectorias o
¶ptimas para la inversi¶
on'),...
xlabel('t'),...
ylabel('i t'),...
grid,...
legend('Senda', 'SS')
83
css2=css*ones(T-1,1)';
figure(6)
ct = yt-it;
plot(1:1:T,ct,'.-',1:1:T-1,css2,'r')
title('Trayectorias o
¶ptimas para el consumo'),...
xlabel('t'),...
ylabel('c t'),...
grid,...
legend('Senda', 'SS')
% Guarda los resultados principales en el file vfout.mat
save E:/vfoutm V0 G kt ct it yt
En la ventana de Matlab escribimos vfiterm y vamos a obtener:
Regla de decisión óptima
Función de valor en la iteración 178 con error 0.009833
240
100
90
220
80
200
o kss
70
180
kp=g(k)
v(k)
60
50
40
160
140
30
120
20
100
10
0
100
120
140
160
180
200
80
80
220
100
120
140
k
160
k
180
200
220
240
Trayectorias óptimas para el Producto
Trayectorias óptimas para el Capital
64
Senda
Senda
190
SS
SS
62
180
60
170
58
y
t
150
k
t
160
56
140
54
130
52
120
50
110
100
48
10
20
30
40
50
t
60
70
80
90
0
100
10
20
30
40
50
t
84
60
70
80
90
100
Trayectorias óptimas para el consumo
Trayectorias óptimas para la inversión
52
15.5
Senda
Senda
50
SS
15
SS
48
14.5
46
14
t
c
i
t
44
13.5
42
13
40
12.5
38
12
36
34
11.5
0
10
20
30
40
50
60
70
80
90
100
t
85
0
10
20
30
40
50
t
60
70
80
90
100
6
Problemas de Optimizaci¶
on
En econom¶³a lo que los agentes quieren es optimizar, ya sea que las empresas quieran minimizar sus costos y maximizar sus ganacias, los consumidores quieres maximizar su utilidad,
los jugadores maximizar sus pagos y el planeador social maximizar el bienestar social. Los
econ¶ometras tambi¶en utilizan m¶etodos de optimizaci¶on, como procedimientos de m¶³nmos
cuadrados ordinarios, m¶etodos de momentos y estimaci¶on de m¶axima verosimilitud.
Vamos a examinar problemas de minimizaci¶on, ya que si queremos maximizar f , es lo
mismo que minimizar ¡f; adem¶as de que Matlab tiene un comando para minimizar.
El problema de optimizaci¶on m¶as general es el de minimizar una funci¶on objetivo sujeto
a restricciones de igualdad o desigualdad:
min f (x)
x
s:a: g(x) = 0
h(x) · 0
en donde f : Rn ! R es nuestra funci¶on objetivo, g : Rn ! Rm es el vector de m¡
restricciones con igualdad y h : Rn ! Rl es el vector de l¡ restricciones con desigualdad.
Algunos ejemplos son la maximizaci¶on de la utilidad sujeto a una restricci¶on presupuestal
lineal, o problemas m¶as complejos de agente-principal. Vamos a suponer que f; g y h son
continuas, aunque este supuesto veremos que no siempre es tan necesario.
Todos los m¶etodos de optimizaci¶on buscan dentro del espacio de posibles opciones,
generando una sucesi¶on de guesses que deben de converger a la verdadera soluci¶on. Veremos algunos m¶etodos para resolver estos problemas, se~
nalando tambi¶en sus ventajas y
desventajas.
6.1
Minimizaci¶
on unidimensional
El problema de optimizaci¶on m¶as sencilla es el problema sin restricciones escalar
min f(x)
x2R
en donde f : R ! R. Veremos estos problemas unidimensionales porque ilustran las t¶ecnicas
b¶asicas de una manera clara y porque muchos problemas multidimensionales se reducen a
resolver problemas unidimensionales.
86
6.1.1
M¶
etodo de acorchetamiento
Este es el m¶as sencillo dentro de los m¶etodos de comparaci¶on. Suponer que hemos encontrado
los puntos a; b y c tal que
a < b < c; f(a); f (c) > f (b)
(6.1)
:
f(x)
a
m1 d
b
m2
c
x
Gr¶a¯ca 6.1. El m¶etodo de acorchamiento
La gr¶a¯ca 6.1 ilustra este caso, en el cual dados los supuestos sabemos que existe un
m¶³nimo local en alguna parte del intervalo [a; c], es decir que tenemos un m¶³nimo acorchado.
Lo que haremos es hacer otro subintervalo en donde se encuentre el m¶³nimo, y lo hacemos
escogiendo otro punto d 2 (a; c): Supongamos que d 2 (a; b], entonces si f (d) > f (b),
entonces sabemos que existe un m¶³nimo en alguna parte de [d; c]. Si f(d) < f (b), como en la
gr¶a¯ca, entonces existe un m¶³nimo en [a; b]. Analogamente para el caso en que d 2 [b:c): Por
lo que para encontrar nuestro m¶³nimo vamos a proceder de esta manera hasta determinar
un criterio de cuando detenernos. Veamos el algor¶³tmo:
Paso 1: Encontrar a < b < c, tal que f(a); f (c) > f (b); y escoger un criterio para
detenernos "
Paso 2: Escoger d: Si b ¡ a < c ¡ b; entonces poner d =
b+c
;
2
en otro caso d =
a+b
2
Paso 3: Computar f (d)
Paso 4: Escoger los nuevos (a; b; c): Si d < b y f (d) > f (b), entonces reemplazar
(a; b; c) con (d; b; c). Si d < b y f(d) < f (b), entonces reemplazar (a; b; c) con (a; d; b). Si
87
d > b y f (d) < f (b), entonces reemplazar (a; b; c) con (b; d; c). En otro caso, reemplazar
(a; b; c) con (a; b; d).
Paso 5: Detener si c ¡ a < " , en otro caso volver al paso 2.
>Pero c¶omo encontrar a; b y c? Lo que hacemos es escoger un punto x0 ; una contante
® > 1 y un paso ¢, luego computar f (x0 ); f(x0 § ®¢); f (x0 § ®2 ¢); :::;hasta que tengamos
nuestros tres puntos que satisfacen las condiciones (6.1).
Las ventajas de este m¶etodo es que funciona con cualquier funci¶on continua y acotada
en un intervalo ¯nito,pero el m¶etodo es lento, adem¶as de que s¶olo encuentra m¶³nimos locales,
por ejemplo en la gr¶a¯ca vemos que puede converger a m1 en vez de al m¶³nimo global m2 :
6.1.2
M¶
etodo de Newton-Raphson
Para funciones C 2 este m¶etodo es muy usado, y la idea principal es empezar en un punto
a y encontrar el polinomio cuadr¶atico p(x) que aproxima a f (x) en a en segundo grado,
implicando que
p(x) ´ f (a) + f 0 (a)(x ¡ a) +
f 00 (a)
(x ¡ a)2
2!
Si f 00 (a) > 0, entones p es convexo. Ahora vamos a aproximar f econtrando el punto
xm que minimiza p(x): El valor de minimizaci¶on es xm = a ¡
f 0 (a)
:
f 00 (a)
Si p(x) es una buena
aproximaci¶on global de f (x), entonces xm estar¶a cerca del m¶³nimo de f; o al menos, esperamos que xm est¶e m¶as cerca del m¶³nimo que a: En general, vamos a repetir este paso hasta
que estemos muy cerca de nuestro objetivo.
Estas consideraciones motivan el m¶etodo de Newton-Raphson:
xn+1 = xn ¡
f 0 (xn )
f 00 (xn )
Notemos que a este m¶etodo no le importa si la funci¶on es c¶oncava o convexa, es decir
que el m¶etodo de Newton-Raphson est¶a tratando de encontra puntos cr¶³ticos en general, es
decir, soluciones a f 0 (x) = 0: Si este m¶etodo converge a x¤ ; tenemos que comprobar con
f 00 (x¤ ) si x¤ es un m¶aximo o m¶³nimo local.
La ventaja de este m¶etodo es que si converge, lo hace muy r¶apido si tenemos un buen
guess inicial, sin embargo puede no converger o puede ser que f 00 (x) sea muy dif¶³cil de
calcular.
88
El algor¶³tmo del m¶etodo Newton-Raphson es muy sencillo:
Paso 1: Escoger el guess inicial y los par¶ametros para detenernos ± > 0 y " > 0
f 0 (xk )
f 00 (xk )
Paso 2: xk+1 = xk ¡
Paso 3: Si jxk+1 ¡ xk j < "(1 + jxk j) y jf 0 (xk )j
Ejemplo: Restricci¶on presupuestaria del consumidor
Consideremos ahora un sencillo ejemplo de maximizaci¶on de utilidad para ver estos
dos m¶etodos. Suponer que el consumidor tiene $1 para gastar en el bien x y en el bien y:
El precio de x es de $2, el de y es de $3, y su funci¶on de utilidad es x0:5 + 2y 0:5 : Sea µ
la cantidad gastada en x, y 1 ¡ µ la cantidad gastada en y: Por lo tanto el problema del
consumidor se reduce a
µ ¶0:5
¶0:5
µ
µ
1¡µ
max
+2
µ
2
3
cuya soluci¶on es µ ¤ =
3
11
= 0:272727:::
Para usar el m¶etodo de comparaci¶on, necesitamos acotar o acorchar el valor ¶optimo
de µ: En este caso sabemos que µ¤ 2 [0; 1], por lo tanto escogermos a = 0; b = 0:5; c = 1
inicialmente y generamos una sucesi¶on de (a; b; c) triples, en donde en cada iteraci¶on b es el
valor m¶³nimo de µ que se ha probado. Los nuevos m¶³nimos van a formar la sucesi¶on
0:5; 0:5; 0:25; 0:25; 0:25; 0:25; 0:281; 0:281; 0:266; 0:2737
Notemos que la tasa de convergencia es lenta como lo hab¶³amos dicho.
Veamos este ejemplo, pero usando el m¶etodo de Newton-Raphson: Si µ0 = 0:5 es
nuestro guess inicial entonces la iteraci¶on nos dar¶³a:
>> NewRaph('funNR','funNRpr','funNRpr2',.5,0.01,10)
El metodo de Newton-Raphson converge
paso
x
ypr
1.0000
0.5000
0.3165
2.0000
0.2596
0.0229
3.0000
0.2724
0.0005
4.0000
0.2727
0.0000
ans =
0.5000 0.2596 0.2724 0.2727
Claramente vemos que el m¶etodo de Newton-Raphson converge m¶as r¶apido, siempre
y cuando nuestro guess sea cercano a la verdadera soluci¶on.
89
Ahora tratemos en clase de replicar esta soluci¶on escribiendo los alumnos sus propias
funciones de Newton-Raphson y de la funci¶on objetivo con sus derivadas.
function [x,y]=NewRaph(fun,funpr,funpr2,x1,tol,max);
% Resuelve f(x)=0 usando el m¶etodo de Newton
% fun Nombre de la funci¶on entre ap¶ostrofes
% funpr Especi¯ca la derivada de f
% x0 Estimaci¶on inicial
% tol Tolerancia permitida al computar el cero
% maxit N¶
umero m¶aximo de iteraciones permitido
%
% x Vector de aproximaciones de cero
% y Vector de valores de la funci¶on fun(x)
x(1)=x1;
y(1)=feval(fun,x(1));
ypr(1)=feval(funpr,x(1));
ypr2(1)=feval(funpr2,x(1));
for i=2:max
x(i)=x(i-1)-ypr(i-1)/ypr2(i-1);
ypr(i)=feval(funpr,x(i));
if abs(x(i)-x(i-1))<tol
disp('El metodo de Newton-Raphson converge');
break;
end
ypr2(i)=feval(funpr2,x(i));
iter=i;
end
if (iter>=max)
disp('no se encontr¶o el cero a la tolerancia deseada');
end
n=length(x);
k=1:n;
out=[k' x' ypr'];
disp(' paso x ypr')
disp(out)
90
function y=funNR(x)
y=(x/2)^(1/2)+2*((1-x)/3)^(1/2);
function y=funNRpr(x)
y=1/4*2^(1/2)/x^(1/2)-1/3/(1/3-1/3*x)^(1/2);
function y=funNRpr2(x)
y=-1/8*2^(1/2)/x^(3/2)-1/18/(1/3-1/3*x)^(3/2);
6.2
Minimizaci¶
on multidimensional
Vimos la secci¶on anterior para darnos una idea b¶asica de estos m¶etodos, pero pasemos ahora
a verlo con varias variables, como sucede m¶as en la pr¶actica.
El problema sin restricciones ser¶³a
min f(x)
x
para f : Rn ! Rn: Por conveniencia vamos a suponer que f es continua , aunque los
siguientes m¶etodos pueden funcionar a¶
un sin este supuesto.
6.2.1
B¶
usqueda de la malla (grid)
Es el procedimiento m¶as primitivo para encontrar un m¶³nimo de una funci¶on, ya que se
especi¯ca una malla de puntos, de digamos, 100 a 1000 puntos, se evalua f(x) en estos
puntos y se escoge el m¶³nimo. Aunque este m¶etodo parace lento y no muy so¯sticado, a
veces es bueno primero intentar con este m¶etodo, porque tal vez tengamos suerte.
Aunque muchas veces no nos conformemos con estas soluciones, si son muy u
¶tiles, ya
que por ejemplo si estamos tratando de maximizar la funci¶on de m¶axima verosimilitud, los
resultados de estos c¶alculos pueden indicar la curvatura general de la funci¶on de verosimilitud.
Si la b¶
usqueda de la malla indica que la funci¶on es plana en un gran intervalo, no tiene mucho
sentido buscar m¶etodos m¶as so¯sticados. Si la b¶
usqueda de la malla nos indica que hay varios
¶optimos locales, entonces sabes que vamos a tener que trabajar m¶as para encontrar el ¶optimo
global. Adem¶as si el Hessiano en el mejor punto de la malla no es muy bueno, entonces es
91
no es muy probable que los siguientes m¶etodos que veamos sean m¶as efectivos. Sin embargo,
la ventaja de este m¶etodo es que puede dar muy buenos guesses iniciales.
Si sabemos que el objetivo es suave con una buena curvatura y con un u
¶ nico ¶optimo
local, entonces este paso preliminar no es muy u
¶ til. Pero fuera de estos casos especiales,
necesitamos la informaci¶on que nos da esta b¶
usqueda de la malla.
6.2.2
M¶
etodo de Newton-Raphson
Algunas funciones objetivo, como las funciones de utilidad o las funciones de ganacias son
suaves y con varias derivadas, por lo que podemos sacar mucho provecho usando esta informaci¶on acerca de la funci¶on objetivo.
Denotemos
µ
@f
@f
(x); :::;
(x)
rf(x) =
@x1
@xn
µ 2
¶n
@ f
H(x) =
(x)
@xi @xj
i;j=1
¶
como el gradiente y el Hessiano de f en x, respectivamente. Notemos que x es un
vector columna, mientras que el gradiente es un vector ¯la.
El m¶etodo de Newton-Raphson multidimensional es un procedimiento sencillo que nos
da resultados rapidamente a problemas C 2 :
Si f es convexa, pero no cuadr¶atica, generalmente no vamos a poder resolver para el
m¶³nimo. Sin embargo, podemos reemplazar f localmente con una aproximaci¶on cuadr¶atica
y resolver para el m¶³nimo de la aproximaci¶on. En este m¶etodo examinamos una sucesi¶on de
puntos, xk ; en donde en cada iteraci¶on reemplazamos f (x) con una aproximaci¶on cuadr¶atica
de f en xk y escoge xk+1 como el punto cr¶³tico de la aproximaci¶on local. En xk esta aprox:
imaci¶on cuadr¶atica es f (x) = f (xk ) + rf (xk )(x ¡ xk ) + 2!1 (x ¡ xk )0 H(xk )(x ¡ xk ): Si f es
convexa, entonces H(xk ) ser¶a de¯nida positiva y la aproximaci¶on va a tener un m¶³nimo en
xk ¡H(xk )¡1 (rf (xk ))0 : Por lo tanto el m¶etodo de Newton-Raphson ser¶a el siguiente esquema
iterativo:
xk+1 = xk ¡ H(xk )¡1 (rf(xk ))0
El m¶etodo multivariante de Newton-Raphson puede ser muy costoso, porque computar
y guardar el Hessiano es muy espacioso, por lo que se utiliza un nuevo paso, de¯nido como
92
sk = xk+1 ¡ xk ; que es la soluci¶on al problema lineal
H(xk )sk = ¡r(f (xk ))0
(6.2)
Por lo que el procedimiento ser¶³a computar H(xk ); resolver el sistema lineal (6.2) para sk y
poner xk+1 = xk + sk :
Este m¶etodo en su forma multivariante tiene los mismos problemas que el univariante,
aunque tambi¶en tiene la ventaja de que cuando converge, converge cuadr¶aticamente.
6.2.3
Encontrar el cero de un sistema no lineal con minimizaci¶
on
Vamos a ver que los m¶etodos de minimizaci¶on tambi¶en nos ayudan a encontrar el cero de
un sistema de ecuaciones. Empezamos de¯niendo una nueva funci¶on que es la suma de las
cuadrados de las funciones cuyo cero en com¶
un se desea encontrar. Por ejemplo, si queremos
el cero de las funciones
f(x; y) = x2 + y 2 ¡ 1
g(x; y) = x2 ¡ y
de¯nimos la funci¶on h(x; y) como
h(x; y) = (x2 + y 2 ¡ 1)2 + (x2 ¡ y)2
El valor m¶³nimo de h(x; y) = 0, lo cual s¶olo ocurre cuando f (x; y) = 0 y g(x; y) = 0:
Las derivadas para el gradiente son
hx (x; y) = 2(x2 + y 2 ¡ 1)2x + 2(x2 ¡ y)2x
hy (x; y) = 2(x2 + y 2 ¡ 1)2y ¡ 2(x2 ¡ y)
Con la siguiente funci¶on de Matlab obtenemos los siguientes valores:
>> ffmin('ex min','ex min g',[0.5 0.5],0.00001,100)
0
0.5000
0.5000
1.0000
0.8750
0.6250
-0.2683
2.0000
0.7451
0.6113
-0.0360
93
3.0000
0.7925
0.6190
-0.0080
4.0000
0.7845
0.6178
-0.0002
5.0000
0.7866
0.6181
-0.0000
las iteraciones han convergido
ans =
0.7860 0.6180
en donde ffmin, ex min y ex min g est¶an de¯nidas como:
function xmin = ffmin(my func, my func g, x0,tol, max it)
% Encuentra el m¶
³nimo de una sistema de ecuaciones
iter = 1;
x old = x0;
disp([ 0 x0])
while (iter <= max it)
dx = feval(my func g, x old);
x new = x old + dx;
z0 = feval(my func, x old);
z1 = feval(my func, x new);
ch = z1- z0;
while( ch >= 0)
dx = dx/2;
if ( norm(dx) < 0.00001)
break;
else
x new = x old + dx;
z1 = feval(my func, x new);
ch = z1 -z0;
end
end
if (abs(ch) < tol)
disp('las iteraciones han convergido')
break;
end
disp([ iter x new ch])
94
x old = x new;
iter = iter +1;
end
xmin = x new;
function y=ex min(x)
y=(x(1)^2+x(2)^2-1)^2+(x(1)^2-x(2))^2;
function dy=ex min g(x)
dy=[-(2*(x(1)^2+x(2)^2-1)*2*x(1)+2*(x(1)^2-x(2))*2*x(1));
-(4*(x(1)^2+x(2)^2-1)*x(2)-2*(x(1)^2-x(2)))]';
6.2.4
Optimizaci¶
on sin restricciones
En las siguientes secciones vamos a ver las funciones que Matlab ya tiene escritas en su
Optimization Toolbox.
Vamos a comenzar con un sencillo ejemplo de un problema de minimizaci¶on sin restricciones. Por ejemplo,
min f(x1 ; x2 ) = (x1 ¡ 2)4 + (x1 ¡ 2x2 )2
En este ejemplo, claramente f(x1 ; x2 ) ¸ 0 y f (2; 1) = 0, por lo tanto (2; 1) es un ¶optimo
global. El gradiente es:
rf(x1 ; x2 ) =
"
4(x1 ¡ 2)3 + 2(x1 ¡ 2x2 )
¡4(x1 ¡ 2x2 )
y en este caso es muy obvio que rf (2; 1) = 0:
#
En esta caja de herramientas que tenemos de optimizaci¶on en Matlab tenemos dos
opciones para resolver estos problemas sin restricciones:
² fminsearch
² fminunc
95
Ambas funciones necesitan un archivo .m o una funci¶on inline para evaluar la funci¶on
objetivo y un guess inicial de la soluci¶on.
Supongamos para nuestro ejemplo un guess inicial de x0 = 0; as¶³
>>f=inline('(x(1)-2)^4+(x(1)-2*x(2))^2');
>>x=fminsearch(f,[0 0])
Los resultados ser¶an x1 = 2 y x2 = 1; con f(x) = 2:9563e ¡ 017
Probemos ahora fminunc:
>>x=fminunc(f,[0 0])
En este caso los resultados son x1 = 1:9897 y x2 = 0:9948; con f(x) = 1:1274e ¡ 008
Podemos ver que este resultado no es exacto, debido a que la funci¶on es plana alrededor
del minimizador. De hecho la funci¶on objetivo es cercana al cero en la soluci¶on reportada
por Matlab. Podr¶³amos cambiar los par¶ametros de tolerancia para mejorar nuestra soluci¶on,
pero en la pr¶actica no hace mucho sentido hacer esto. Podemos notar que Matlab tambi¶en se
queja por la falta de la informaci¶on del gradiente, por lo que no puede aplicar alg¶
un m¶etodo
de regi¶on de con¯anza10 . Sin embargo, esto no es ning¶
un problema, porque el gradiente se
puede estimar num¶ericamente. Notemos que nosotros tambi¶en podemos darle el gradiente a
Matlab para que lo utilice en la computaci¶on de la soluci¶on. Sin embargo, para incorporar
el gradiente tenemos que crear un archivo .m con la funci¶on objetivo y su gradiente:
function [f,g] = myfun(x)
f=(x(1)-2)^4+(x(1)-2*x(2))^2;
g=[4*(x(1)-2)^3+2*(x(1)-2*x(2)), -4*(x(1)-2*x(2))];
>> options = optimset('GradObj','on','tolfun',1e-13);
>> x = fminunc(@(x) myfun(x),[0 0],options)
El resultado es x1 = 1:9997 y x2 = 0:9998; el cual podemos notar que es un poco m¶as
preciso.
Sin embargo, como hemos visto a trav¶es del curso, computar el gradiente anal¶³ticamente
puede ser complicado, por lo cual le podemos preguntar a Matlab que compare e~
n gradiente
10
Si hay tiempo se ver¶a estos m¶etodos en el curso.
96
que nosotros le demos con una estimaci¶on num¶erica. Lo u
¶nico que tenemos que hacer reset
las opciones y activar en las opciones el comando derivativecheck. Es decir,
>>options=optimset;
>>options=optimset('gradobj','on','derivativecheck','on','tolfun',1e-13);
>>x = fminunc(@(x) myfun(x),[0 0],options)
Podemos observar que la diferencia entre los gradientes es muy chica (3:72e ¡ 007); por
lo que preferir¶³amos en la pr¶actica dejar que Matlab evalue el gradiente, en vez de nosotros
calcularlo num¶ericamente. Veamos los valores si Matlab eval¶
ua el gradiente:
>> options = optimset('tolfun',1e-13);
>> [X,FVAL,EXITFLAG,OUTPUT,GRAD]=fminunc(f,[0 0],options)
GRAD =
1.0e-008 *
0.4604
-0.9182
Mientras que el gradiente evaluado nos da:
>> g=inline('[4*(x(1)-2)^3+2*(x(1)-2*x(2)), -4*(x(1)-2*x(2))]');
>> g(x)
ans =
1.0e-009 *
-0.1256 0.0000
Notemos que tambi¶en podemos darle otro valor del gradiente y checar si las derivadas
est¶an cerca, y cu¶an cerca est¶an. Como de costumbre si queremos obtener m¶as informaci¶on
de la funci¶on ponemos help fminunc.
6.2.5
Optimizaci¶
on con restricciones de igualdad y desigualdad
Consideremos el problema convexo de
97
min x21 + x22
s:a: x1 ¸ 0
x2 ¸ 3
x1 + x2 = 4
Primero vamos a comenzar escribiendo el Lagrangeano, recordando que como es una
minimizaci¶on las restricciones de desigualdad deber¶³an de reescribirse como ¡x1 · 0 y
¡x2 + 3 · 0
L(x; ¹; ¸) = x21 + x22 ¡ ¹1 x1 ¡ ¹2 (x2 ¡ 3) + ¸(x1 + x2 ¡ 4)
Las condiciones de Kuhn-Tucker ser¶³an:
2x1 ¡ ¹1 + ¸ = 0
2x2 ¡ ¹2 + ¸ = 0
x1 ¸ 0
x2 ¸ 3
x1 + x2 = 4
¹1 x1 = 0;
¹1 ¸ 0
¹2 (x2 ¡ 3) = 0;
¹2 ¸ 0
Podemos observar que vamos a tener varios casos, por lo que vamos a ver cada uno
de los casos de manera anal¶³tica y luego procederemos a ver como lo podemos resolver en
Matlab.
Caso 1: ¹1 = ¹2 = 0
En este caso las restricciones de desigualdad no se toman en cuenta en el Lagrangeano,
98
y las condiciones de estacionaridad del sistema son:
2x1 + ¸ = 0
2x2 + ¸ = 0
x1 + x2 = 4
Con la cual x1 = x2 = 2; sin embargo se viola la segunda restricci¶on de desigualdad.
Caso 2: ¹1 ; ¹2 6= 0
Las condiciones de factibilidad y exclusi¶on11 inmediatamente conllevan a x1 = 0; x2 = 3
violando la restricci¶on de igualdad.
Caso 3: ¹1 6= 0; ¹2 = 0
En este caso tenemos:
x1 = 0
x2 = 4
¸ = ¡2x2 = ¡8
¹1 = ¸ = ¡8
Se violan las condiciones de Kuhn-Tucker, porque ¹1 es negativa.
Caso 4: ¹1 = 0; ¹2 6= 0
Vamos a obtener:
x2 = 3
x1 = 1
¸ = ¡2
¹2 = 4
11
Ver notas de Ricard Torres "Resultados de Optimizaci¶
on"
99
El cual satisface todas las condiciones!!!
Dado que es un problema conexo, hemos obtenido el ¶optimo global.
Veamos ahora c¶omo lo podemos resolver con Matlab. El comando quadprog resuleve
problemas cuadr¶aticos tal que
1
min xT Hx + f T x
2
s:a:Ax · b
Aeq x = beq
l · x·u
Para nuestro problema algunas entradas est¶an vac¶³as. De¯namos nuestro problema de
la siguiente manera:
>>H=2*eye(2);
>>f=[0 0];
>>Aeq= [1 1];
>>beq=4;
>>lb=[0 3];
>>[x,f,exitflag,output,lambda]=quadprog(H,f,[],[],Aeq,beq,lb);
>>x
>>lambda.eqlin
>>lambda.lower
Los argumentos de salida incluyen las variables de decisi¶on ¶optimas, el valor ¶optimo
de la funci¶on objetivo e informaci¶on adicional, como el valor de lambda. Podemos ver que
los resultados son los mismos que nosotros calculamos anal¶³ticamente, pero en Matlab lo
resolvimos muy r¶apido.
7
Optimizaci¶
on din¶
amica
Hemos visto lo que es la programaci¶on din¶amica num¶erica para un ejemplo de crecimiento
neocl¶asico. Vamos a ver a continuaci¶on otro ejemplo cl¶asico para estos problemas de opti100
mizaci¶on: Problema de comerse un pastel (cake-eating)
7.1
Cake-eating example12
Supongamos que se nos presenta un pastel de tama~
no W1 : En cada momento del tiempo
t = 1; 2; ::: uno puede comer un poco del pastel, pero debe guadar el resto. Sea ct el consumo
en el periodo t, y sea u(ct ) la utilidad de este consumo. Notemos que la funci¶on de utilidad
no est¶a indexada con el tiempo, es decir, que las preferencias son estacionarias. Podemos
suponer que u(¢) vive en los reales, es diferenciable, estricatmente creciente y estricatmente
c¶oncava. Representemos la utilidad a lo largo del ciclo de vida como
max
1
fct g1
1 ;fWt g2
1
X
¯ t u(ct )
t=1
en donde 0 · ¯ · 1 es el factor de descuento.
Vamos a suponer que el pastel no se deprecia, es decir, que no se hecha a perder, pero
tampoco crece, por lo que la evoluci¶on del pastel a trav¶es del tiempo, es decir, la ecuaci¶on
de transici¶on es Wt+1 = Wt ¡ ct para t = 1; 2; :::
Escribiendo este problema en programaci¶on din¶amica nos quedar¶³a la siguiente ecuaci¶on
de Bellman:
V (W ) = max u(c) + ¯V (W ¡ c)
c2[0;W ]
8W
en donde V (W ) es el valor del problema de horizonte in¯nito empezando con el tama~
no
de un pastel W: Por lo tanto en cada periodo el agente escoge su consumo actual y por lo
tanto reduce el tama~
no del paste a W 0 = W ¡cm como en la ecuaci¶on de transici¶on. Notemos
que en este problema las variables de estado son el tama~
no del pastel W dado al principio
de cada periodo. La variable de control es el consumo actual.
Alternativamente podemos especi¯car el problema para que en vez de escoger el consumo actual, podamos escoger el estado de ma~
nana:
V (W ) = max
u(W ¡ W 0 ) + ¯V (W 0 )
0
W 2[0;W ]
12
Ejemplo tomado de Adda y Cooper (2003).
101
8W
Claramente cualquiera de las dos especi¯caciones nos da el mismo resultado. Sin embargo, esta u
¶ ltima expresi¶on es m¶as f¶acil algegraicamente, por lo que trabajaremos con ¶esta.
A esta expresi¶on se le llama ecuaci¶on funcional. Recordemos que la ecuaci¶on de Bellman
no est¶a indexada por el tiempo, caracter¶³stica de la estacionaridad. Notemos que toda la
informaci¶on pasada ya est¶a registrada en W .
Las condiciones de primer orden son:
u0 (c) = ¯V 0 (W 0 )
Resolviendo el problema13 , vamos a obtener la funci¶on de pol¶³tica:
c = Á(W );
W 0 = '(W ) ´ W ¡ Á(W )
En general, no es muy posible que podamos encontrar las soluciones de las funci¶on de
valor y de pol¶³tica, por lo que vamos a tratar de proponer unas formas espec¶³¯cas y ver que
nuestras conjeturas son correctas. Por ejemplo, vamos a suponer que u(c) = ln(c), por lo
que podemos conjeturar que la ecuaci¶on funcional es de la forma
V (W ) = A + B ln(W )
por lo que s¶olo tenemos que determinar dos par¶ametros. Resolviendo vamos a obtener
que:
c = W (1 ¡ ¯)
W 0 = ¯W
Por lo que la pol¶³tica ¶optima nos dice que ahorrarremos una fracci¶on constante del pastel
y comeremos la fracci¶on restante.El programa propuesto por Adda y Cooper (2003) es el
siguiente:
%%%%%%%%%%% DETERMINISTIC CAKE EATING MODEL %%%%%%%%%%%%%%%%%%%%%
%%% THIS FILE SOLVES THE DETERMINISTIC CAKE EATING MODEL PRESENTED
%%% IN CHAPTER 2 OF THE ADDA/COOPER BOOK (SEE SECTION 2.4.1, PAGE 16)
clear
dimIter=30; % number of iterations
13
Para m¶
as detalles, ver Adda y Cooper (2003).
102
beta=0.9; % discount factor
K=0:0.01:1; % grid over cake size
[rk,dimK]=size(K);
V=zeros(dimK,dimIter); % initialize the value function.
Rows are cake size,
columns the iteration
iter=1;
for iter=1:dimIter % loop over all current cake sizes
iter
aux=zeros(dimK,dimK)+NaN;
for ik=1:dimK % loop over next period cake size
for ik2=1:(ik-1)
aux(ik,ik2)=log(K(ik)-K(ik2))+beta*V(ik2,iter);
end
end
V(:,iter+1)=max(aux')'; % optimizing over size of next period cake
end
% computing the optimal consumption as a function of initial cake size
[Val,Ind]=max(aux');
optK=K(Ind);
optK=optK+Val*0;
optC=K'-optK';
% Plotting the value function after each iterations
figure(1)
plot(K,V);
xlabel('Size of Cake');
ylabel('Value Function');
figure(2)
plot(K,[optC K'],'LineWidth',2)
xlabel('Size of Cake');
ylabel('Optimal Consumption');
text(0.4,0.65,'45 degree line','FontSize',18)
text(0.4,0.13,'Optimal Consumption','FontSize',18)
Cuyas gr¶a¯cas correspondientes son:
103
1
0.9
0.8
Optimal Consumption
0.7
45 degree line
0.6
0.5
0.4
0.3
0.2
Optimal Consumption
0.1
0
0
0.1
0.2
0.3
0.4
0.5
0.6
Size of Cake
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.6
Size of Cake
0.7
0.8
0.9
1
0
-5
-10
Value Function
-15
-20
-25
-30
-35
-40
-45
Para ver cu¶ales fueron los valores ¶optimos, s¶olo escribimos:
>>optK
>>optC
104
8
El Modelo Estoc¶
astico
Consideremos ahora un modelo un poco distinto, en el cual la productividad de la ¯rma
representativa esta sujeta a un shock tecnol¶ogico µ. Es decir, en cada per¶³odo t, tenemos:
Yt = µt f (Kt )
en donde µ t es una variable aleatoria que sigue un determinado proceso estoc¶astico. Las
decisiones se toman en cada per¶³odo antes de que el shock correspondiente sea observado.
Para poder aplicar el m¶etodo de programaci¶on din¶amica, suponemos que los shocks
toman un n¶
umero ¯nito q de valores y siguen un proceso de Markov de primer orden. Es
¯
¤
¢
£
¡ 1 2
decir, µ 2 µ ; µ ; :::; µq y P rob µt+1 = µj ¯ µ t = µ i = ¼ ij . Evidentemente, debe cumplirse
P
que qj=1 ¼ ij = 1. La matriz ¦ de orden q £ q con todas las probabilidades ¼ ij se conoce
como la matriz de transici¶on.
Un Equilibrio General Competitivo, Recursivo y Estoc¶astico para esta econom¶³a es un
conjunto de funciones (o planes contingentes) v (k; K; µ), c (k; K; µ), i (k; K; µ), k 0 (k; K; µ),
precios w (K; µ) y r (K; µ) y ley de movimiento ¡ (K; µ) tales que:
i) Dadas las funciones w, r y ¡, la funci¶on de valor v (k; K; µ) resuelve la ecuaci¶on de
Bellman:
v (k; K; µ) = max0
c;i;k
fu (c) + ¯Eµ v (k 0 ; K 0 ; µ0 )g
(8.1)
s:t: c + i = w (K; µ) + r (K; µ) k
k 0 = (1 ¡ ±) k + i
K 0 = ¡ (K; µ)
y las funciones c (k; K), i (k; K) y k 0 (k; K) son reglas de decisi¶on ¶optimas para este problema.
ii) Para todo K y µ, los precios satisfacen las condiciones marginales:
r (K; µ) = µf 0 (K)
105
(8.2)
w (K; µ) = µf (K) ¡ µf 0 (K) K
(8.3)
iii) Para todo K y µ, hay igualdad entre oferta y demanda:
µf (K) = c (K; K; µ) + i (K; K; µ)
(8.4)
iv) Para todo K y µ, ley de movimiento agregada y comportamiento individual son
consistentes:
¡ (K; µ) = k 0 (K; K; µ)
(8.5)
N¶otese que, dada la forma en que hemos de¯nido el proceso estoc¶astico para los shocks
tecnol¶ogicos, podemos expresar el valor esperado en (8.1) como:
0
0
0
Eµi v (k ; K ; µ ) =
q
X
j=1
¢
¡
¼ ij v k 0 ; K 0 ; µj
(8.6)
N¶otese tambi¶en que las reglas de decisi¶on ¶optimas son planes contigentes. Para obtener las
trayectorias correspondientes es necesario simular una historia completa de realizaci¶on de los
shocks tecnol¶ogicos.
Similarmente, de¯nimos un Optimo de Pareto Recursivo y Estoc¶astico como un conjunto de funciones v (k; µ), c (k; µ), i (k; µ), k 0 (k; µ) que resuelven la ecuaci¶on de Bellman del
plani¯cador social:
v (k; µ) = max0
c;i;k
fu (c) + ¯Eµ v (k 0 ; µ0 )g
(8.7)
s:t: c + i = µf (k)
k 0 = (1 ¡ ±) k + i
y nuevamente, en aquellos casos en que se cumplan los supuestos de los Teoremas del Bi106
enestar, resolveremos este problema como una manera indirecta de calcular el Equilibrio
Competitivo.
8.1
Iteraci¶
on de la Funci¶
on de Valor
La ecuaci¶on de Bellman (8.7) puede resolverse num¶ericamente usando el m¶etodo de iteraci¶on
de la funci¶on de valor. Dada una malla (K1; K2; :::; Kp ) para el stock de capital y la malla
(µ1 ; µ2 ; :::; µq ) para los shocks tecnol¶ogicos, de¯nimos la matriz M de orden pq £ p de la
siguiente manera:
2
6
6
6
6
6
6
6
6
6
6
M =6
6
6
6
6
6
6
6
6
6
4
F (K1 ; K1 ; µ 1 )
F (K1 ; K2 ; µ1 )
::::: F (K1 ; Kp ; µ1 )
3
7
7
7
7
F (Kp ; K1 ; µ1 )
F (Kp ; K2 ; µ1 ) ::::: F (Kp ; Kp ; µ1 ) 7
7
F (K1 ; K1 ; µ 2 )
F (K1 ; K2 ; µ2 ) ::::: F (K1 ; Kp ; µ2 ) 7
7
7
::::::::::::::::::::::::: ::::::::::::::::::::::: ::::: :::::::::::::::::::::::: 7
7
7
F (Kp ; K1 ; µ2 )
F (Kp ; K2 ; µ2 ) ::::: F (Kp ; Kp ; µ2 ) 7
7
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 7
7
7
F (K1 ; K1 ; µ q ) F (K1 ; K2 ; µq ) ::::: F (K1 ; Kp ; µq ) 7
7
:::::::::::::::::::::::: :::::::::::::::::::::::: ::::: :::::::::::::::::::::: 7
5
F (Kp ; K1 ; µq ) F (Kp ; K2 ; µq ) ::::: F (Kp ; Kp ; µq )
::::::::::::::::::::::::: :::::::::::::::::::::::: ::::: ::::::::::::::::::::::
en donde F (x; y; µ) = u [µf (x) + (1 ¡ ±) x ¡ y]. Nuevamente, imponemos las restricciones
tecnol¶ogicas en la matriz M haciendo::
Mp(l¡1)+i;j ´ F (Ki ; Kj ; µl ) = 0
si Kj > µl f (Ki ) + (1 ¡ ±) Ki
Una vez constru¶³da M , escojemos cualquier vector columna V 0 2 Rpq (por ejemplo,
V 0 = 0) e iniciamos el siguiente algoritmo con s = 0:
¡
¢
1. Dado el vector V s = V1s ; V2s ; :::; Vpqs , calcular una matriz W s de q £ p de acuerdo a:
2 Pq
Pq
s
s
j=1 ¼ 1j Vj(p¡1)+2
j=1 ¼ 1j Vj(p¡1)+1
P
P
6 q
q
s
s
6 j=1 ¼2j Vj(p¡1)+1
j=1 ¼ 2j Vj(p¡1)+2
s
6
W =6
4 ::::::::::::::::::::::::::: :::::::::::::::::::::::::::
Pq
Pq
s
s
j=1 ¼ qj Vj(p¡1)+2
j=1 ¼ qj Vj(p¡1)+1
107
3
s
¼
V
1j
j=1
j(p¡1)+p
Pq
7
s
7
:::::
j=1 ¼ 2j Vj(p¡1)+p 7
7
::::: ::::::::::::::::::::::::::: 5
Pq
s
:::::
j=1 ¼ qj Vj(p¡1)+p
:::::
Pq
2. Dados el vector W s y la matriz M , calcular la matriz M ¤ de orden pq £p de la siguiente
manera:
2
6
6
6
6
6
6
6
6
6
6
¤
M = 6
6
6
6
6
6
6
6
6
6
4
s
M11 + ¯W11
s
M12 + ¯W12
s
::::: M1p + ¯W1p
3
7
::::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::::: ::::: ::::::::::::::::::::::::::::::: 7
7
7
s
s
s
Mp1 + ¯W11
Mp2 + ¯W12
::::: Mpp + ¯W1p
7
7
7
s
s
s
Mp+1;1 + ¯W21
Mp+1;2 + ¯W22
::::: Mp+1;p + ¯W2p
7
7
::::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::::: ::::: ::::::::::::::::::::::::::::::: 7
7
7
s
s
s
M2p;1 + ¯W21
M2p;2 + ¯W22
::::: M2p;p + ¯W2p
7
7
::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 7
7
7
s
s
s
7
M(q¡1)p+1;1 + ¯Wq1
M(q¡1)p+1;2 + ¯Wq2
::::: M(q¡1)p+1;p + ¯Wqp
7
::::::::::::::::::::::::::::::: ::::::::::::::::::::::::::::::: ::::: ::::::::::::::::::::::::::::::: 7
5
s
s
s
Mqp;1 + ¯Wq1
Mqp;2 + ¯Wq2
::::: Mqp;p + ¯Wqp
3. Dada la matriz M ¤ , calcular el vector V s+1 como:
V s+1 = max M ¤
donde el m¶aximo es calculado por ¯la.
4. Evaluar kV s+1 ¡ V s k. Si la norma es mayor que el criterio de tolerancia permitido,
regresar al paso 1 con s = s + 1. En caso contrario, el procedimiento termina con la
soluci¶on V = V s+1 .
La regla de decisi¶on ¶optima se calcula como antes:
G = arg max M ¤
en donde M ¤ es la matriz obtenida en el segundo paso de la u
¶ ltima iteraci¶on.
8.2
Simulaci¶
on de las Trayectorias Optimas
Con el m¶etodo anterior, obtenemos la funci¶on de valor V y la regla de decisi¶on ¶optima
G evaluadas para cada stock de capital y cada shock tecnol¶ogico. Estos vectores est¶an
ordenados de la siguiente manera:
108
2
6
6
6
6
6
6
6
6
6
6
6
6
V =6
6
6
6
6
6
6
6
6
6
6
4
v (K1; µ1 )
3
2
7
v (K2; µ1 ) 7
7
7
::::::::::::::: 7
7
v (Kp; µ1 ) 7
7
7
v (K1; µ2 ) 7
7
7
::::::::::::::: 7
7
v (Kp; µ2 ) 7
7
7
::::::::::::::: 7
7
7
v (K1; µq ) 7
7
::::::::::::::: 7
5
6
6
6
6
6
6
6
6
6
6
6
6
G=6
6
6
6
6
6
6
6
6
6
6
4
v (Kp; µq )
g (K1; µ1 )
3
7
g (K2; µ1 ) 7
7
7
::::::::::::::: 7
7
g (Kp; µ1 ) 7
7
7
g (K1; µ2 ) 7
7
7
::::::::::::::: 7
7
g (Kp; µ2 ) 7
7
7
::::::::::::::: 7
7
7
g (K1; µq ) 7
7
::::::::::::::: 7
5
g (Kp; µq )
en donde l = g (Ki ; µj ) nos indica que, dado el stock de capital y la realizaci¶on del shock
tecnol¶ogico, el plani¯cador social escoger¶a un stock de capital para el siguiente per¶³odo igual
a Kl .
Para obtener una secuencia ¶optima de capital k0 ; k1 ; :::; kT debemos entonces simular
primero una historia de realizaci¶on de shocks tecnol¶ogicos µ0 ; µ1 ; :::; µT , usando por supuesto
la informaci¶on contenida en la matriz de transici¶on ¦. Cualquier lenguaje de programaci¶on
escoje n¶
umeros aleatrios de acuerdo a una distribuci¶on espec¶³¯ca. Por ejemplo, podemos
pedir que escoja µ0 asignando la misma probabilidad 1=q a cada uno de los posibles valores
¡
¢
en la malla µ1 ; µ2 ; :::; µq . A continuaci¶on, pedimos que simule recursivamente una secuencia
¢
¡
µ1 ; :::; µT de manera tal que, dado µt = µi , escoja µt+1 entre los valores µ1 ; µ2 ; :::; µq con
una distribuci¶on de probabilidad (¼ i1 ; ¼ i2 ; ::; ¼ iq ).
Una vez conocida la secuencia de shocks tecnol¶ogicos, podemos simular una secuencia
¶optima para el capital k0 ; k1 ; :::; kT y para el resto de variables, usando la regla de decisi¶on
¶optima G. Estas secuencias representan series de tiempo, cuyos principales estad¶³sticos
(media, varianza, patr¶on de correlaciones) pueden ser calculados y comprados con los datos.
N¶otese, sin embargo, que las series de tiempo obtenidas dependen de la particular realizaci¶on
de los shocks tecnol¶ogicos. Por eso, se recomienda hacer un n¶
umero grande de simulaciones
y usar el promedio (entre simulaciones) de cada uno de los estad¶³sticos mencionados.
109
8.3
Limitaciones
El m¶etodo de iteraci¶on de la funci¶on de valor tiene la ventaja de ser bastante preciso cuando
la malla de valores para las variables de estado es su¯cientemente ¯na. Adem¶as, permite
incorporar incertidumbre de una manera natural. Sin embargo, este m¶etodo tiene a su
vez ciertas limitaciones. Primero, su aplicabilidad es mayor para econom¶³as en las que los
Teoremas del Bienestar se cumplen y podemos entonces resolver el problema del plani¯cador
social. Ya hemos visto las di¯cultades que se presentan al tratar de resolver directamente el
Equilibrio Competitivo.
En segundo lugar, al aumentar el n¶
umero de variables de estado y/o el n¶
umero de
puntos en la malla, la dimensi¶on del problema (de la matriz M ) crece exponencialmente
volviendo este m¶etodo muy lento y costoso en t¶erminos computacionales.
Por u
¶ltimo, los m¶etodos basado en la programaci¶on din¶amica solo permiten introducir shocks estoc¶asticos de una determinada estructura, b¶asicamente procesos de Markov
de primer orden con n¶
umero ¯nito de estados. Otros procesos estoc¶asticos m¶as generales
requieren m¶etodos distintos, como se ver¶a m¶as adelante.
110
8.4
Ejemplo
Consideremos la versi¶on simple del modelo del agente representativo, con las formas funcionales:
f (k) = µAk ®
u(c) = log(c)
y los valores de los par¶ametros:
¯ = 0:95
A = 10
® = 0:35
± = 0:06
y un stock de capital inicial k0 igual a la mitad de su valor k ¤ en un equilibrio estacionario.
En cada per¶³odo la funci¶on de producci¶on se enfrenta a un shock tecnol¶ogico µ t . Este shock
puede tomar tres valores: µ 2 f0:5; 1; 1:5g y sigue un proceso de Markov discreto de orden
uno con la siguiente matriz de transici¶on:
2
0:5 0:3 0:2
3
7
6
7
¦=6
0:1
0:7
0:2
5
4
0 0:4 0:6
Tenemos que resolver el problema del planeador usando iteraci¶on de la funci¶on de valor con
una malla de 300 puntos igualmente espaciados de 0 a 1:5k ¤ , en donde k ¤ representa el stock
de capital de estado estacionario en la econom¶³a determinpistica. Comenzando con un stock
de capital inicial k0 = k ¤ y un shock µ0 = 0:5, vamos a simular y gra¯car la realizaci¶on
(series de tiempo) de las trayectorias ¶optimas del capital, consumo, producto, salario real y
tasa de retorno del capital. Finalmente, usando un promedio de 100 simulaciones, computar
la desviaci¶on est¶andar del logaritmo de cada una de las series obtenidas y su correlaci¶on con
el producto.
El programa en un archivo .m ser¶³a:
% vfiter2.m Program to solve a simple version of the planner¶s problem
% with technology shocks using value function iteration.
% equation to solve is:
%
111
The Bellman
% v(k,z) = max fln(th*A*k^alpha+(1-delta)*k-k') + beta*E(v(k',th'))g
% k'
%
% with k0 given
%
% Output:
file vfout2.mat with value function, optimal decision rule,
% and sequences for main variables.
% Programa obtenido de Carlos Urrutia.
%
% Needs function shock.m to simulate shocks
clear
% Parameters of the model
alpha = 0.35;
beta = 0.95;
delta = 0.06;
A = 10;
sh = [0.5 1 1.5];
pi = [0.5 0.3 0.2;0.1 0.7 0.2;0 0.4 0.6];
% Other required parameters (maximum number of iterations, size of the
% grid for k, tolerance criterion, number of periods for simulations, size
of the
% shock vector)
maxit = 1000;
p = 300;
crit = 1e-2;
T = 100;
q = 3;
% Computes steady state capital (kss) for an economy without uncertainty,
% initial capital (k0) and defines a grid for k
kss = ((A*beta*alpha)/(1-(1-delta)*beta))^(1/(1-alpha));
k0 = kss;
k = linspace(0,1.5*kss,p);
% Sets initial guess for value function (V0=0) and constructs matrix M
V0 = zeros(p*q,1);
112
e = ones(1,p);
f = ones(1,q);
g = ones(1,p*q);
state = ones(p*q,2);
aux = f'*k;
aux = aux(:);
state(:,1) = aux;
aux = sh'*e;
aux = aux(:);
state(:,2) = aux;
M = log(max((A*state(:,2).*(state(:,1).^alpha)+...
(1-delta)*state(:,1))*e-g'*k,1e-8));
I = eye(q);
E = I;
for i=1:p-1
E=[E;I];
end
% Value function iteration
iconv=0;
it=1;
while (iconv==0 & it<maxit)
[V1,G] = max((M+beta*(E*(pi*reshape(V0,q,p))))');
V1 = V1';
G = G';
if norm(V0-V1)<crit
iconv=1;
end
disp(['Iteration= ',num2str(it),' Error= ',num2str(norm(V0-V1))])
V0 = V1;
it = it+1;
end
% Simulates a history of shocks (uses function shock.m)
theta = shock(2,pi,T);
% Simulates optimal paths for capital (kt), consumption (ct),
113
% investment (it) and output (yt)
G = reshape(G,q,p);
kt = zeros(T+1,1);
ind = zeros(T,1);
[aux,ind(1)] = min(k<k0);
kt(1) = k(ind(1));
for t = 1:T
ind(t+1) = G(theta(t),ind(t));
kt(t+1) = k(ind(t+1));
end
figure(1)
plot(1:T,sh(theta(1:T)))
title('Realizations of Technology Shock')
xlabel('t')
ylabel('z t')
figure(2)
plot(1:T,kt(1:T))
title('Optimal Path for Capital')
xlabel('t')
ylabel('k t')
figure(3)
yt = A*sh(theta(1:T))'.*kt(1:T).^alpha;
plot(1:T,yt)
title('Optimal Path for Output')
xlabel('t')
ylabel('y t')
figure(4)
it = kt(2:T+1)-(1-delta)*kt(1:T);
plot(1:T,it)
title('Optimal Path for Investment')
xlabel('t')
ylabel('i t')
figure(5)
ct = yt-it;
114
plot(1:T,ct)
title('Optimal Path for Consumption')
xlabel('t')
ylabel('c t')
figure(6)
rt = alpha*A*sh(theta(1:T))'.*kt(1:T).^(alpha-1);
plot(1:T,rt)
title('Optimal Path for the Interest Rate')
xlabel('t')
ylabel('r t')
figure(7)
wt = yt-rt(1:T).*kt(1:T);
plot(1:T,wt)
title('Optimal Path for the Real Wage')
xlabel('t')
ylabel('w t')
% Saves main results in file vfout2.mat
save vfout2 V0 G kt ct it yt
% Computes statistics averaging 100 simulations
theta1 = zeros(T,100);
kt1 = zeros(T+1,100);
yt1 = zeros(T,100);
ct1 = zeros(T,100);
it1 = zeros(T,100);
rt1 = zeros(T,100);
wt1 = zeros(T,100);
ind1 = zeros(T+1,100);
for i=1:100
x = rand; % Draws initial th 0 assigning same probability
if(x<=0.333) % to each possible state
aux=1;
elseif(x>0.333 & x<=0.666)
aux=2;
else
115
aux=3;
end
theta1(:,i) = shock(aux,pi,T)';
[aux,ind1(1,i)] = min(k<k0);
kt1(1,i) = k(ind1(1,i));
for t=1:T
ind1(t+1,i) = G(theta1(t,i),ind1(t,i));
kt1(t+1,i) = k(ind1(t+1,i));
end
yt1(:,i) = A*sh(theta1(1:T,i))'.*kt1(1:T,i).^alpha;
it1(:,i) = kt1(2:T+1,i)-(1-delta)*kt1(1:T,i);
ct1(:,i) = yt1(:,i)-it1(:,i);
rt1(:,i) = alpha*A*sh(theta1(1:T,i))'.*kt1(1:T,i).^(alpha-1);
wt1(:,i) = yt1(:,i)-rt1(:,i).*kt1(1:T,i);
end
% Standard Deviationsof the log and Correlations with Output
dtk = std(log(kt1));
dtc = std(log(ct1));
dti = std(log(it1));
dty = std(log(yt1));
dtr = std(log(rt1));
dtw = std(log(wt1));
aux = zeros(2,2);
for i=1:100
aux = corrcoef(kt1(1:T,i),yt1(:,i));
corrk(i) = aux(1,2);
aux = corrcoef(ct1(:,i),yt1(:,i));
corrc(i) = aux(1,2);
aux = corrcoef(it1(:,i),yt1(:,i));
corri(i) = aux(1,2);
aux = corrcoef(rt1(:,i),yt1(:,i));
corrr(i) = aux(1,2);
aux = corrcoef(wt1(:,i),yt1(:,i));
corrw(i) = aux(1,2);
116
end
desvk = sum(dtk)/100;
desvc = sum(dtc)/100;
desvi = sum(dti)/100;
desvy = sum(dty)/100;
desvr = sum(dtr)/100;
desvw = sum(dtw)/100;
correlk = sum(corrk)/100;
correlc = sum(corrc)/100;
correli = sum(corri)/100;
correlr = sum(corrr)/100;
correlw = sum(corrw)/100;
disp(' ')
disp('Standard Deviation of the log:
')
disp(' ')
disp(['- Output = ',num2str(desvy)])
disp(['- Consumption = ',num2str(desvc)])
disp(['- Investment = ',num2str(desvi)])
disp(['- Capital = ',num2str(desvk)])
disp(['- Interest Rate = ',num2str(desvr)])
disp(['- Real Wage = ',num2str(desvw)])
disp(' ')
disp('Correlation with Output of:
')
disp(' ')
disp(['- Consumption = ',num2str(correlc)])
disp(['- Investment = ',num2str(correli)])
disp(['- Capital = ',num2str(correlk)])
disp(['- Interest Rate = ',num2str(correlr)])
disp(['- Real Wage = ',num2str(correlw)])
disp(' ')
Notemos que necesitamos el archivo shock.m para obtener los shocks:
function th = shock(th0,pi,T)
% Function to simulate a history of T shocks, given an initial value of
% the shock equal to th0 and Markov transition matrix pi
117
th = ones(1,T);
th(1) = th0;
cum = cumsum(pi')';
for i = 2:T
x = find(rand < cum(th(i-1),:));
th(i) = x(1);
end
Por lo tanto los resultados obtenidos ser¶an:
Optimal Path for Capital
Realizations of Technology Shock
300
1.5
1.4
280
1.3
1.2
260
t
1
k
z
t
1.1
240
0.9
220
0.8
0.7
200
0.6
0.5
0
10
20
30
40
50
t
60
70
80
90
180
100
0
10
20
30
40
50
t
60
70
80
90
100
70
80
90
100
Optimal Path for Investment
Optimal Path for Output
110
40
100
30
90
20
80
t
i
y
t
10
70
0
60
-10
50
-20
40
30
0
10
20
30
40
50
t
60
70
80
90
100
118
-30
0
10
20
30
40
50
t
60
Optimal Path for Consumption
Optimal Path for the Interest Rate
95
0.18
90
0.16
85
0.14
80
0.12
r
c
t
t
75
70
0.1
65
0.08
60
0.06
55
50
0
10
20
30
40
50
t
60
70
80
90
100
70
80
90
100
Optimal Path for the Real Wage
80
70
w
t
60
50
40
30
20
0
10
20
30
40
50
t
60
Iteration= 184 Error= 0.0098782
Standard Deviation of the log:
- Output = 0.35472
- Consumption = 0.19007
- Investment = 1.2214
- Capital = 0.21661
- Interest Rate = 0.32783
- Real Wage = 0.35472
Correlation with Output of:
- Consumption = 0.78726
119
0.04
0
10
20
30
40
50
t
60
70
80
90
100
- Investment = 0.88858
- Capital = 0.4203
- Interest Rate = 0.72377
- Real Wage = 1
120
9
Aproximaci¶
on
9.1
Aproximaci¶
on por m¶³nimos cuadrados
9.1.1
Lineal
Como es bien sabido, uno de los m¶etodos m¶as comunes para aproximar datos es el de querer
minimizar alguna medida de la diferencia entre la funci¶on de aproximaci¶on y los datos. El
m¶etodo de m¶³nimos cuadrados intenta minimizar la suma (sobre todos los puntos de la base
de datos) de los cuadrados de las diferencias entre el valor de la funci¶on y los puntos de los
datos.
Las ventajas de los m¶³nimos cuadrados ya son bien conocidas, por ejemplo, que las
diferencias negativas no se cancelan con las positivas, es f¶acil diferenciar la funci¶on y las
diferencias peque~
nas se hacen m¶as chicas y las grandes las magni¯ca.
Veamos primero como podemos hacer con Matlab una aproximaci¶on lineal con algunos
datos, para despu¶es ver una aproximaci¶on cuadr¶atica y una c¶
ubica. Dado que queremos una
aproximaci¶on lineal, vamos a querer encontrar la mejor funci¶on f(x) = ax + b que aproxime
los datos.
Consideremos los puntos (1; 2:1); (2; 2:9); (5; 6:1) y (7; 8:3) como se ven en la siguiente
gr¶a¯ca:
10
8
6
4
2
0
2
4
x
6
8
10
Gr¶a¯ca 9.1. Datos para el ejemplo de la aproximaci¶on de m¶inimos cuadrados lineal
Para encontrar los mejores coe¯cientes de la l¶³nea recta que mejor aproxima tenemos
que minimizar la funci¶on de los errores al cuadrado, es decir,
E = [f(x1 ) ¡ y1 ]2 + [f(x2 ) ¡ y2 ]2 + [f(x3 ) ¡ y3 ]2 + [f(x4 ) ¡ y4 ]2
= [ax1 + b ¡ y1 ]2 + [ax2 + b ¡ y2 ]2 + [ax3 + b ¡ y3 ]2 + [ax4 + b ¡ y4 ]2
121
El m¶³nimo de esta funci¶on diferenciable sabemos que se va a encontrar en donde la
derivada sea cero, es decir,
@E
= 0
@a
@E
= 0
@b
La soluci¶on depende de varias cantidades que se pueden computar con los datos. Las
cuatro sumatorias necesarias son:
Sxx =
Sx =
Sxy =
Sy =
4
X
i=1
4
X
i=1
4
X
i=1
4
X
x2i
xi
x i yi
yi
i=1
Las desconocidas a y b se pueden encontrar resolviendo el siguiente sistema de ecuaciones conocidas como ecuaciones normales:
aSxx + bSx = Sxy
aSx + b(4) = Sy
La soluci¶on a este sistema de ecuaciones es:
4Sxy ¡ Sx Sy
4Sxx¡ Sx Sx
Sxx Sy ¡ Sxy Sx
b =
4Sxx¡ Sx Sx
a =
Para nuestro ejemplo, el sistema de ecuaciones que tenemos que resolver es:
a(79) + b(15) = 96:5
a(15) + b(4) = 19:4
122
Cuya soluci¶on es a = 1:044 y b = 0:9352:
Veamos ahora cual es la funci¶on en Matlab que nos va a ayudar a encontrar la mejor
aproximaci¶on lineal de m¶³nimos cuadrados:
function s = Lin LS(x, y)
% funci¶
on de regresi¶
on lineal
% input x y y como vectores columna o fila
% (se convierten en forma de columnas si es necesario)
m = length(x); x = x(:); y = y(:);
sx = sum(x); sy = sum(y);
sxx = sum(x.*x); sxy = sum(x.*y);
a = ( m*sxy -sx*sy) / (m*sxx - sx^2)
b = ( sxx*sy - sxy*sx) / (m*sxx - sx^2)
table = [x y (a*x+b) (y - (a*x+b))];
disp(' x y (a*x+b) (y - (a*x+b))')
disp(table), err = sum(table( :
, 4) .^2)
s(1) = a; s(2) = b;
n=m;
xx=[x(1):(x(n)-x(1))/100:x(n)];
plot(x, y, 'r*', xx, a*xx+b, 'b-');
grid on;
Por lo que en la ventana de Matlab escribimos:
>> x=[1,2,5,7]
>> y=[2.1,2.9,6.1,8.3]
y nos arroja los siguientes resulatdos:
>> Lin LS(x,y)
a =
1.0440
b =
0.9352
x y (a*x+b) (y - (a*x+b))
1.0000 2.1000 1.9791 0.1209
2.0000 2.9000 3.0231 -0.1231
5.0000 6.1000 6.1549 -0.0549
123
7.0000 8.3000 8.2429 0.0571
err =
0.0360
ans =
1.0440 0.9352
9
8
7
6
5
4
3
2
1
0
0
9.1.2
1
2
3
4
5
6
7
8
Cuadr¶
atica
Usando el mismo enfoque que para el caso lineal, ahora vamos a querer aproximar nuestros
datos a la funci¶on cuadr¶atica f (x) = ax2 + bx + c; cuya funci¶on de error es la misma que
para el caso anterior, es decir,
n
X
E=
[f (xi ) ¡ yi ]2
i=1
Igualando las derivadas de E con respecto a a; b y c a cero, obtenemos las siguientes
ecuaciones normales para a; b y c :
a
a
n
X
x4i
i=1
n
X
+b
x3i + b
x3i
i=1
n
X
+c
x2i + c
n
X
i=1
x2i
+b
n
X
i=1
n
X
x2i
=
xi =
i=1
i=1
i=1
a
n
X
n
X
xi + c(n) =
i=1
n
X
i=1
n
X
i=1
n
X
i=1
124
x2i yi
x i yi
yi
La funci¶on en Matlab es:
function z = Quad LS(x,y)
% regresi¶
on cuadr¶
atica
% input x y y como vectores columna o fila
% (se convierten en forma de columnas si es necesario)
n = length(x); x = x(:); y = y(:);
sx = sum(x); sx2 = sum(x.^2);
sx3 = sum(x.^3); sx4 = sum(x.^4);
sy = sum(y); sxy = sum(x.* y);
sx2y = sum(x.*x.*y);
A = [ sx4 sx3 sx2,
sx3 sx2 sx,
sx2 sx n]
r = [sx2y sxy sy]'
z = Anr;
a = z(1), b = z(2), c = z(3)
table = [x y (a*x.^2 + b*x + c) (y - (a*x.^2 + b*x + c))];
disp(' x y (a*x.^2 + b*x + c) (y - (a*x.^2 + b*x + c))')
disp(table)
err = sum(table( :
, 4) .^2)
xx=[x(1):(x(n)-x(1))/100:x(n)];
plot(x, y, 'r*', xx, a*xx.^2+b*xx+c, 'b-');
grid on;
Consideremos los datos
x = [0; 0:5; 1; 1:5; 2]
y = [0; 0:19; 0:26; 0:29; 0:31]
Cuyos resultados son:
>> Quad LS(x,y)
A =
22.1250 12.5000 7.5000
12.5000 7.5000 5.0000
7.5000 5.0000 5.0000
r =
125
2.2000
1.4100
1.0500
a =
-0.1086
b =
0.3611
c =
0.0117
x y (a*x.^2 + b*x + c) (y - (a*x.^2 + b*x + c))
0 0 0.0117 -0.0117
0.5000 0.1900 0.1651 0.0249
1.0000 0.2600 0.2643 -0.0043
1.5000 0.2900 0.3091 -0.0191
2.0000 0.3100 0.2997 0.0103
err =
0.0012
ans =
-0.1086
0.3611
0.0117
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0
0.2
0.4
0.6
0.8
1
126
1.2
1.4
1.6
1.8
2
9.1.3
C¶
ubica
Por u
¶ltimo veamos c¶omo ser¶³a para el caso c¶
ubica, , en donde el sistema de ecuaciones para
determinar los coe¯cientes de para el mejor ajuste de la funci¶on f (x) = ax3 + bx2 + cx + d
es:
n
X
a
x6i
+b
i=1
n
X
a
a
x5i
+c
i=1
x5i + b
n
X
x4i + b
i=1
n
X
a
x4i + c
x4i
+d
n
X
x3i + c
i=1
n
X
x3i + b
n
X
x3i + d
x3i
=
n
X
x2i + d
i=1
n
X
x2i =
n
X
x3i yi
n
X
x2i yi
i=1
xi =
i=1
n
X
xi + d
i=1
n
X
i=1
i=1
x2i + c
i=1
n
X
i=1
i=1
n
X
i=1
n
X
i=1
i=1
i=1
n
X
n
X
i=1
1 =
n
X
i=1
n
X
xi yi
yi
i=1
Cuya funci¶on en Matlab es:
function z = Cubic LS(x,y)
% regresi¶
on c¶
ubica,
% input x y y como vectores columna o fila
% (se convierten en forma de columnas si es necesario)
n = length(x); x = x( :
); y = y( :
);
sx = sum(x); sx2 = sum(x .^2);
sx3 = sum(x .^3); sx4 = sum(x .^4);
sx5 = sum(x .^5); sx6 = sum(x .^6);
sy = sum(y); syx = sum(y.* x);
syx2 = sum(y.*x.^2); syx3 = sum(y.*x.^3);
A = [ sx6 sx5 sx4 sx3;
sx5 sx4 sx3 sx2;
sx4 sx3 sx2 sx;
sx3 sx2 sx n]
r = [syx3 syx2 syx sy]'
z = A n r ;
a = z(1); b = z (2); c = z(3); d = z(4);
p = a*x.^3 + b*x.^2 + c*x + d;
table = [ x y p (y - p)];
disp(' x y p(x) = a x ^3 + b x ^2 + c x + d y - p(x)')
127
disp( table )
err = sum(table( :
, 4) .^2)
xx=[x(1):(x(n)-x(1))/100:x(n)];
plot(x, y, 'r*', xx, a*xx.^3+b*xx.^2+c*xx+d, 'b-');
grid on;
Usando los siguientes datos:
x = [¡1 : :2 : 1]
y = [0:05; 0:08; 0:14; 0:23; 0:35; 0:5; 0:65; 0:77; 0:86; 0:92; 0:95]
Obtenemos:
>> Cubic LS(x,y)
A =
2.6259 0 3.1328 0.0000
0 3.1328 0.0000 4.4000
3.1328 0.0000 4.4000 0.0000
0.0000 4.4000 0.0000 11.0000
r =
1.5226
2.2000
2.2800
5.5000
x y p(x) = a x ^3 + b x ^2 + c x + d y - p(x)
-1.0000 0.0500 0.0552 -0.0052
-0.8000 0.0800 0.0708 0.0092
-0.6000 0.1400 0.1352 0.0048
-0.4000 0.2300 0.2364 -0.0064
-0.2000 0.3500 0.3621 -0.0121
0 0.5000 0.5000 -0.0000
0.2000 0.6500 0.6379 0.0121
0.4000 0.7700 0.7636 0.0064
0.6000 0.8600 0.8648 -0.0048
0.8000 0.9200 0.9292 -0.0092
1.0000 0.9500 0.9448 0.0052
err =
128
6.4615e-004
ans =
-0.2550
-0.0000
0.6997
0.5000
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-1
9.1.4
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Funci¶
on de Matlab polyfit
La funci¶on de Matlab polyfit encuentra los coe¯cientes de un polinomio de alg¶
un grado en
espec¶³¯co que mejor se ajusta alos datos, en el sentido de m¶³nimos cuadrados. La funci¶on
es del estilo [p,S]=polyfit(x,y,n), en donde x y y son los vectores de los datos y n es el
grado del polinomio deseado.
Veamos un ejemplo paso a paso:
% regresi¶
on polinomial
% generar los datos
x = -3 :
0.1 :
3 ;
yy1 = sin(x);
yy2 = cos(2*x);
yy = yy1 + yy2;
y = 0.01*round(100*yy);
¶nimos cuadrados del polinomio de grado 6
% encontrar mi
129
%para ajustar los datos
z = polyfit(x, y, 6)
% evalua el polinomio
p = polyval(z, x);
% graficar el polinomio
plot(x, p);
grid on;
hold on;
% graficar los datos
plot(x, y, '+');
hold off;
Los resultados son:
>> use pfit
z =
-0.0232 0.0058 0.3819 -0.1567 -1.5712 0.9907 0.8974
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2.5
-3
-2
-1
0
1
2
3
Notemos que en vez de utilizar 61 puntos, podemos hacer una buena aproximaci¶on con
s¶olo 10 puntos, es decir, usemos
x=[-3:2/3:3] ;
y =[.82,-.77,-1.98,-1.26,.46,1.11,.43,.01,.68,1.1];
y obtenemos los siguientes resultado, aunque vemos que no es tan suave la gra¯ca:
130
>> use pfit2
z =
-0.0241 0.0057 0.3983 -0.1562 -1.6272 0.9908 0.9199
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-2.5
-3
10
-2
-1
0
1
2
3
Integraci¶
on de Monte Carlo
Esta es una t¶ecnica de integraci¶on num¶erica muy intuitiva, y la mejor manera de explicarla
es a trav¶es de un ejemplo. Consideremos el caso en el que queremos determinar el ¶area de
un cuarto de c¶³rculo, por lo cual es una integraci¶on de dos variables:
Ahora escogemos un n¶
umero de puntos dentro de esta caja, primero comenzando
eligiendo un n¶
umero aleatorio repreesentando la distancia en el eje-x y posteriormente escogiendo un n¶
umero aleatorio representando la distancia en el eje-y. Estos dos n¶
umeros
131
representan un par (x; y) que localiza un punto en la caja. Realizando este proceso repetidamente conlleva a una imagen como la siguiente:
Si suponemos que el ¶area del cuarto de c¶³rculo relativamente al ¶area de la caja est¶a
determinada por el n¶
umero de puntos que caen dentro del c¶³rculo en relaci¶on al n¶
umero de
puntos del cuadrado, entonces podemos escribir:
Ac¶{rculo = Acuadrado
Nc¶{rculo
Ncuadrado
Entonces podemos aproximar la integral escogiendo un n¶
umero de puntos aleatorios en
el cuadrado, contando cu¶antos caen dentro del ¶area de inter¶es y luego realizar la ecuaci¶on.
La t¶ecnica de Monte Carlo es muy u
¶ til para casos multidimensionales cuyas formas son algo
extra~
nas.
El comando en Matlab para obtener la integral num¶erica de este ejercicio es:
% this is a plot of the function to be solved, x=0 to 1
for i=1:100 %notar que en esta escala, la gr¶
afica de 0 a 1
%se ver¶
a como de 0 a 100
x(i)=i*.01;
f(i)=sqrt(1-x(i)^2);
end
plot(f)
%
T=10;
for i=1:1000
132
x=rand;
y=2*rand;
f=sqrt(1-x^2);
if y<=f
T=T+1;
end
end
integral=2*T/1000
Podemos ver una simulaci¶on en:
http://homepages.cae.wisc.edu/~blanchar/eps/quad/quadframe.htm
133
11
Bibliograf¶³a
Adda, J. & Cooper, R. (2003): Dynamic Economics: Quantitative Methods and Applications, MIT Press.
Fausett, L.V. (1999)Applied Numerical Analysis Using MATLAB. Prentice Hall, L.
Judd, Kenneth (1998): Numerical Methods in Economics, MIT Press.
Lucas, R. Stokey, N. & Prescott, E.(1989) Recursive Methods in Economic Dynamics
Harvard University Press.
Miranda, M. and Fackler, P. (2002) Applied Computational Economics MIT Press
12
Ap¶
endice
En Matlab para trabajar en la carpeta donde est¶an los archivos que nos interesan, s¶olo
vamos a la barra de herramientas que tiene tres puntos ..., y ah¶³ buscamos nuestra carpeta
de inter¶es.
En Octave escribimos:
cd '/cygdrive/c/'
cd..
ls
y ah¶³ vamos buscando nuestra carpeta escribiendo cada vez cd, por ejemplo, si tenemos
algo en nuestro USB, escribimos cd d
134
Related papers
Development Steps of Avionics and Flight Control System of Flight Vehicle
Jurnal RESTI (Rekayasa Sistem dan Teknologi Informasi)
Objeción de conciencia al aborto e igualdad laboral de los profesionales sanitarios
Ius Canonicum, 2024
Potassium-induced alleviation of salinity stress in Brassica campestris L
Open Life Sciences, 2011
The impact of forest management on forest floor vegetation evaluated by species traits
Community Ecology, 2003
Lumbar ontogenetic allometry and dimorphism in humans. A case for comparison between interspecific and intraspecific scaling
European Journal of Morphology, 2005
Direito ao Desenvolvimento no Estado Neoliberal: a Renda Básica Universal e Incondicional Como Alternativa à Redução da Desigualdade Social no Brasil
Direito Público, 2021
Effect of Push-Pull Factors on Change in Behavior of Students' Shuttle by Using Motorcycle Case Study: Gianyar - Bali
Civil engineering and architecture, 2022
Developing a Robust Bioventing Model
Mathematical and computational applications, 2023
Consumer Attitude, Consumer Trust, and Price Perception on Purchace Decision
Journal of Applied Management and Business (JAMB)
- Find new research papers in:
- Physics
- Chemistry
- Biology
- Health Sciences
- Ecology
- Earth Sciences
- Cognitive Science
- Mathematics
- Computer Science