Mauricio Tejada
ILADES - Universidad Alberto Hurtado
Agosto 2017
Ideas Básicas: La forma más natural de aproximar una derivada es reemplazarla con una diferencia finita. Recordemos que la definición de derivada:
$$ f'(x) = \lim_{h \rightarrow 0} \frac{f(x+h)-f(x)}{h}$$Por tanto la aproximación será:
$$f'(x) = \frac{f(x+h)-f(x)}{h}$$para un $h$ pequeño. Dado que la anterior ecuación representa una aproximación a la definición de derivada, existirá un error de aproximación asociado.
Una alternativa, con un error de aproximación menor, es construir la aproximación usando diferencias a los dos lados:
$$f'(x) = \frac{f(x+h)-f(x-h)}{2h}$$¿Cuán pequeño debe ser $h$? Existe una disyuntiva:
Miranda y Fackler proveen una regla para la aproximación de dos lados: $$h=max(|x|,1)*\sqrt[3]{\epsilon}$$ con $\epsilon$ es denominado machine precision.
Escribamos la función que calcula la jacobiana de una función (tomada del libro de Miranda y Fackler).
function fjac = fdjac(f,x)
h = eps^(1/3)*max(abs(x),1);
xh1 = x+h;
xh0 = x-h;
hh = xh1- xh0; % 2h
for j=1:length(x);
xx = x;
xx(j) = xh1(j);
f1 = feval(f,xx);
xx(j) = xh0(j);
f0 = feval(f,xx);
fjac(:,j) = (f1-f0)/hh(j);
end
end
Probemos la función aproximando la derivada de la función $$f(x) = 2x^2+x-1$$
Sabemos que: $$f'(x) = 4x + 1$$ y evaluando en $x=2$ tenemos $f'(2) = 9$.
% Definimos la funcion
fx = @(x) 2*x^2 + x - 1;
% Aplicamos la función fdjac (nota que pasamos una función como argumento)
derv = fdjac(fx,2)
derv = 9.0000
Tomemos ahora la siguiente función: $$f(x,y,z) = xyz + 2xy + 2yz + x + y + z$$
y usemos la función fdjac
para encontrar el gradiente de la función.
% Definimos la función
fxyz = @(x) x(1)*x(2)*x(3) + 2*x(1)*x(2) + 2*x(2)*x(3) + x(1) + x(2) + x(3);
% Aplicamos la función fdjac (nota que pasamos una función como argumento)
x0 = [1; 1; 1];
derv = fdjac(fxyz,x0)
derv = 4.0000 6.0000 4.0000
Ahora escribamos el sistema:
$$f(x,y) = x^2 + y^2 - 2$$$$g(x,y) = xy$$en un m-file.
function f = sys(var)
f = zeros(length(var),1);
x = var(1);
y = var(2);
f(1) = x^2 + y^2 - 2;
f(2) = x*y;
end
Ahora usemos la función fdjac
para encontrar la matriz jacobiana del sistema evaluada en el punto $(1,1)$:
clear all;
x0 = [1; 1];
g = fdjac('sys',x0)
g = 2.0000 2.0000 1.0000 1.0000
Un alternativa es usar el toolbox Adaptive Robust Numerical Differentiation. Este toolbox es gratuito y permite, bajo la lógica anterior pero usando métodos adaptativos, calcular la derivada (1er, 2do, 3er y 4to orden), el vector gradiente, la matriz jacobiana y la matriz hessiana. Las opciones básicas son:
derivest(fun,punto,'deriv',orden)
con orden
1,2,3,4 y fun
una función escalargradest(fun,punto)
con fun
una función escalar.jacobianest(fun,punto)
con fun
un sistema de ecuaciones.hessian(fun,punto)
con fun
una función escalar.Para detalles de todas las opciones ver la documentación del toolbox.
Nota: La última versión del toolbox soporta matlab R2014b o superior. Para una versión anterior solicitarla por email al profesor.
Ideas Básicas: En muchas aplicaciones económicas es necesario calcular la integral definida de una función $f(x)$ con respecto a una función de ponderadores $w(x)$ sobre el intervalo $I$ en $R^n$.
$$\int_I f(x)w(x)dx$$En algunos casos la función de ponderadores podría ser la unidad, $w(x)=1$, de tal manera que la integral representa el área bajo la función $f(x)$. En otras aplicaciones, $w(x)$ podría ser una función de densidad de tal manera que la integral representa $E[f(x)]$.
Los métodos conocidos como cuadraturas aproximan la integral de la función con una suma ponderada de valores de la función:
$$\int_I f(x)w(x)dx \approx \sum^{n}_{i=0}w_if(x_i)$$La elección de los ponderadores $w_i$ y de los nodos $x_i$ definen el método.
Para ganar intuición vamos a analizar sólo dos versiones simples de las cuadraturas de Newton-Cotes para una función univariada.
Newton-Cotes
Buscamos calcular: $$\int_a^b f(x) dx$$
Matlab implementa estos dos métodos de aproximación de la integral con los siguientes comandos:
quad(fun,a,b)
donde fun
es una función escalar, implementa la regla de Simpson.trapz(y)
donde $y_i=f(x_i)$, implementa la regla basa de trapezoides.Ejemplo:
% Usando Quad
fx = @(x) 1./(x.^3-2*x+5);
Q = quad(fx,0,2);
disp(Q);
0.4213
x = 0:0.01:2;
f = fx(x);
area(x,f);
% Usando Trapz
Q2 = trapz(x,f);
disp(Q2);
0.4213
Notas:
quad
es evaluada de manera vectorial por lo que las divisiones, las multiplicaciones y las potencias deben explicitarse como operaciones elemento por elemento.quad
sólo permite integrar la función en límites finitos. Para transformar los límites de la integral usar cambio de variable.trapz
es una función vectorial. Evaluamos la función en el intervalo de puntos en que queremos integrar.quad
funciona mejor que trapz
con funciones que tiene un comportamiento suave y sin discontinuidades.Ejemplo: Para integrar
$$\int^{\infty}_{-\infty} f(x) dx$$usar $x=\frac{t}{1-t^2}$ de tal manera que $dx=\frac{1+t^2}{(1-t^2)^2} dt$ y $t=\frac{\sqrt{1+4x^2}-1}{2x}$. Note que:
$$\int^{\infty}_{-\infty} f(x) dx = \int^{1}_{-1} f\left(\frac{t}{1-t^2} \right) \frac{1+t^2}{(1-t^2)^2} dt$$Éste y otro ejemplos pueden ser encontrado en la página de Wikipedia sobre integración numérica.
En versiones más recientes de Matlab se incorporó la función integral
que usa algoritmos más eficientes y precisos. Además permite integrar en límites -Inf
e Inf
.
integral(fun,a,b)
donde fun
es una función escalar.En el futuro, Matlab removerá quad
y sólo quedara la función integral
.
Ejemplo: Resolvamos el siguiente problema:
$$\int^{\infty}_{0} e^{-x^2} (log(x))^2 dx$$f = @(x) exp(-x.^2).*log(x).^2;
Q3 = integral(f,0,Inf);
disp(Q3);
1.9475
Integrales Dobles
Para integrales dobles usamos el función integral2
. La sintaxis es como sigue:
integral2(fun,xmin,xmax,ymin,ymax)
Por ejemplo, integremos la siguiente función:
$$ \int^1_0 \int^{1-x}_{0}\frac{1}{\sqrt{(x + y)} (1 + x + y)^2} dy dx$$Para implementar el problema anterior debemos definir dos funciones:
fun = @(x,y) 1./(sqrt(x + y).*(1 + x + y).^2)
ymax = @(x) 1 - x
fun = function_handle with value: @(x,y)1./(sqrt(x+y).*(1+x+y).^2) ymax = function_handle with value: @(x)1-x
Q4 = integral2(fun,0,1,0,ymax);
disp(Q4)
0.2854
Q4 = integral2(fun,0,1,0,@(x) 1-x)
Q4 = 0.2854