10 Diferenciación e Integración Numérica

10.1 Diferenciación

Ideas Básicas: La forma más simple y 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) \approx \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) \approx \frac{f(x+h)-f(x-h)}{2h}\]

¿Cuán pequeño debe ser \(h\)? Existe una disyuntiva:

Cuando \(h\) es muy pequeño, errores de redondeo pueden llevar a un resultado sin sentido.
Cuando \(h\) es muy grande, generar una aproximación muy pobre.

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 una función para calcular la primera derivada de una función arbitraria:

function primera_derivada(fun,x)
    
    h = (eps()^(1/3))*max(abs(x), 1.0)
    
    return (fun(x+h)-fun(x-h))/(2*h)
end

primera_derivada (generic function with 1 method)

f(x) = x^2
primera_derivada(f,2)

3.9999999999658735

Sabesmos que:

derf(x) = 2*x
derf(2)

Ahora escribamos la función que calcula el gradiente de una función (adaptada del código MATLAB del libro de Miranda y Fackler).

function gradiente(fun,x)
    
    nv = length(x)
    h = (eps()^(1/3))*max.(abs.(x), 1.0)

    xh1 = x+h
    xh0 = x-h
    hh  = xh1 - xh0  # 2h
    
    fderiv = zeros(nv)

    if nv == 1
        f1 = fun(xh1)
        f0 = fun(xh0)
        fderiv = (f1-f0)/hh
    else
        for j=1:nv
            xx = copy(x)
            xx[j] = xh1[j]
            f1 = fun(xx)
            xx[j] = xh0[j]
            f0 = fun(xx)
            fderiv[j] = (f1-f0)/hh[j]
        end
    end
    
    return fderiv
end

gradiente (generic function with 1 method)

Probemos la función aproximando la derivada de: \[f(x) = 2x^2+x-1\]

# Definimos la funcion
fx(x) = 2*x^2 + x - 1

fx (generic function with 1 method)

# Aplicamos la función gradiente (nota que pasamos una función como argumento)
derv = gradiente(fx,2.0)

8.999999999954165

Ahora probemos la función \[f(x,y) = 2x^2+2y^2+x*y-1\]

fxy(x) = 2*x[1]^2+2*x[2]^2+x[1]*x[2]-1

fxy (generic function with 1 method)

x = [3.0; 2.0]
dervxy = gradiente(fxy,x)

2-element Array{Float64,1}:
 13.99999999995111 
 10.999999999862492

Tomemos ahora la siguiente función: \[f(x,y,z) = xyz + 2xy + 2yz + x + y + z\]

y usemos gradiente() para encontrar el gradiente de la función.

# Definimos la funcion
fxyz(x) = x[1]*x[2]*x[3] + 2*x[1]*x[2] + 2*x[2]*x[3] + x[1] + x[2] + x[3]

fxyz (generic function with 1 method)

# Aplicamos la función fgrad
x0 = [1.0; 1.0; 1.0]
derivxyz = gradiente(fxyz,x0)

3-element Array{Float64,1}:
 4.000000000110005
 5.999999999944997
 4.000000000110005

Calculus

Este ejemplo sencillo ayuda a entender la lógica de las derivada numéricas. Existe un módulo en Julia que permite calcular derivadas, gradientes y jacobinas. También permite el cálculo de derivadas de orden superior a 1. El modelo se llama Calculus y se instala usando:

Pkg.add("Calculus")

Primera Derivada: derivative()
Segunda Derivada: second_derivative()
Gradiente: Calculus.gradient()
Matriz Hessiana: hessian()

Para detalles de todas las opciones ver la documentación de cada función.

using Calculus

Partamos con la función \(fx: R \rightarrow R\)

dx = derivative(fx, 2.0)

8.999999999918636

d2x = second_derivative(fx, 2.0)

3.9999990266181173

Es posible generar una función general con la derivada y evaluarla posteriormente en cualquier punto de interés:

dx = derivative(fx)

#1 (generic function with 1 method)

dx(2.0)

8.999999999918636

Ahora apliquemos el cálculo del gradiente a la función \(fxy: R^2 \rightarrow R\)

dfxy = Calculus.gradient(fxy,[3.0, 2.0])

2-element Array{Float64,1}:
 13.999999999938623
 10.999999999819071

La matriz hessiana (matriz de segundas derivadas) para la misma función se calcula como:

d2fxy = Calculus.hessian(fxy,[3.0, 2.0])

2×2 Array{Float64,2}:
 4.00001   0.999998
 0.999998  3.99999

De la misma forma apliquemos el cálculo del gradiente a la función \(fxy: R^2 \rightarrow R\)

dfxyz = Calculus.gradient(fxyz,[1.0, 1.0, 1.0])

3-element Array{Float64,1}:
 4.000000000094215
 5.999999999921312
 4.000000000094215

y su hessiana:

d2fxyz = Calculus.hessian(fxyz,[1.0, 1.0, 1.0])

3×3 Array{Float64,2}:
 7.6074e-6   2.99999      1.00001   
 2.99999    -6.99636e-7   2.99999   
 1.00001     2.99999     -1.66143e-5

ForwardDiff

Otra alternativa en Julia para la aproximación numérica de las derivadas es el paquete ForwardDiff, el mismo que implementa métodos de diferenciación automática en su modo hacia adelante (forward mode automatic differentiation, AD). Este método es un algoritmo que usa de forma exacta las operaciones básicas que realiza un computador y una sucesión acumulada de derivadas sobre la base de la regla de la cadena. Para mayor detalle sobre el método ver la entrada de forward mode automatic differentiation en Wikipedia.

using ForwardDiff

Usemos ahora ForwarDiff sobre la función \[f(x,y) = 2x^2+2y^2+x*y-1\]

fxy(x) = 2*x[1]^2+2*x[2]^2+x[1]*x[2]-1

fxy (generic function with 1 method)

ForwardDiff.gradient(fxy,[3.0; 2.0])

2-element Array{Float64,1}:
 14.0
 11.0

ForwardDiff.hessian(fxy,[3.0; 2.0])

2×2 Array{Float64,2}:
 4.0  1.0
 1.0  4.0

ForwardDiff implementa también el calculo de la matriz jacobiana de una sistema de ecuaciones (esto es, la matriz conformada por los gradientes de cada ecuación). Consideremos el siguiente sistema definido en \(Fxy: R^2 \rightarrow R^2\):

\[f(x,y) = x^2 + y^2 - 2\] \[g(x,y) = xy\]

function sistema(var)

    f = similar(var)

    x = var[1]
    y = var[2]

    f[1] = x^2 + y^2 - 2
    f[2] = x*y

    return f    
end

sistema (generic function with 1 method)

ForwardDiff.jacobian(sistema,[2.0,2.0])

2×2 Array{Float64,2}:
 4.0  4.0
 2.0  2.0

10.2 Integración

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.

Newton-Cotes: Usa polinomios para aproximar la función entre los nodos.
Cuadraturas Gausianas: Elige los ponderadores y los nodos (puntos en \(x\)) de manera de machear ciertos momentos.
Monte Carlo: Usa ponderadores aleatorios y nodos equidistribuidos.

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\]

Regla basadas en trapezoides:
- \(x_i = a + i \times h\) con \(h=(b-a)/n\) e \(i=0,..,n\).
- En el subintervalo \([x_i,x_{i+1}]\) aproximar: \[\int_{x_i}^{x_{i+1}} f(x) dx \approx \frac{h}{2}[f(x_i)+f(x_{i+1})]\]
- Sumando las subareas: \[\int_a^b f(x)dx \approx \sum^{n}_{i=0}w_if(x_i)\] con \(w_1=w_n=h/2\) y \(w_i\).

function integral_trapezoide(fun,x,n)
   
    a = x[1]
    b = x[2]
    
    h = (b-a)/n
    i = collect(range(0,stop=n,length=n+1))
    x = a*ones(n+1) + i*h
    
    y = map(fun, x)
    
    return 0.5 * h * (2 * sum(y[2:length(y)-1])+y[1]+y[length(y)])
end

integral_trapezoide (generic function with 1 method)

f(x)=2*x
int_res = integral_trapezoide(f,[0,1],10)

1.0000000000000002

Regla de Simpson:
- \(x_i = a + i \times h\) con \(h=(b-a)/n\) e \(i=0,..,n\).
- En el par de subintervalos \([x_{2j-1},x_{2j}]\) y \([x_{2j},x_{2j+1}]\) el área bajo la función cuadrática aproximada que pasa por esos tres puntos es: \[\int_{x_{2j-1}}^{x_{2j+1}} f(x) dx \approx \frac{h}{3}[f(x_{2j-1})+4f(x_{2j})+f(x_{2j+1}]\]
- Sumando las subareas: \[\int_a^b f(x)dx \approx \sum^{n}_{i=0}w_if(x_i)\] con \(w_1=w_n=h/3\), \(w_4h/3\) para \(i\) para y \(w_i=2h/3\) para \(i\) impar.

using LinearAlgebra

function integral_simpson(fun,x,n)
   
    a = x[1]
    b = x[2]
    
    h = (b-a)/n
    i = collect(range(0,stop=n,length=n+1))
    x = a*ones(n+1) + i*h
    
    y = map(fun, x)
    
    return h/3 * (2.0 * dot((ones(length(y)-2)+rem.(1:length(y)-2,2)),(y[2:(length(y)-1)])) +y[1]+y[length(y)])
end

integral_simpson (generic function with 1 method)

f(x)=2*x
int_res = integral_simpson(f,[0,1],10)

1.0

QuadGK

Existen varios paquetes desarrollado en Julia para implementar algoritmos de integración numérica. Uno de ellos es QuadGK, el cual implementa diferentes versiones de Cuadraturas Gausianas para aproximar la integral. Existen tres funciones disponibles en este paqueta: QuadGK.quadgk(), QuadGK.gauss() y QuadGK.kronrod(). La sintaxis general es:

(int,err) = QuadGK.quadgk(f, xmin, xmax)
(int,err) = QuadGK.gauss(f, xmin, xmax)
(int,err) = QuadGK.kronrod(f, xmin, xmax)

El resultado es un par de la forma \((int,err)\), donde \(int\) es la integral calculada numéricamente y \(err\) es la cota superior del error absoluto.

La desventaja de QuadGK() es que solo permite integrar funciones de una dimensión.

using QuadGK

f(x) = 2*x
int_res, err_res = QuadGK.quadgk(f, 0.0, 1.0)

(1.0, 0.0)

Nota: El paquete QuadGK() solo permite integrar la función en límites finitos. Para transformar los límites de la integral usar cambio de variable. 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.

Cubature

Otra alternativa es Cubature. Este paquete está bastante desarrollado y permite aproximar integrales tanto en una como en múltiples dimensiones. La característica más importante del paquete es que el algoritmo usado usa la idea de integración adaptativa, la cual significa que el algoritmo se irá adaptando usando un número creciente de puntos hasta conseguir convergencia en términos de error tolerado de aproximación. Cubature usa dos enfoques que difieren en la forma en que ellos obtienen convergencia: h-adaptivity y p-adaptativity. El primer enfoque (h) debiera ser aplicado cuando se sabe poco respecto de la función a integrar, mientras que el segundo (p) debiera ser usado con funciones bien comportadas y con un número pequeño de dimensiones.

La sintaxis de similar a la usada con quadgk:

(int,err) = hquadrature(f, xmin, xmax)
(int,err) = pquadrature(f, xmin, xmax)

Usemos nuevamente el ejemplo de integrar \(f(x)=2x\) en el intervalo \([0,1]\):

using Cubature

f(x) = 2*x
int, err = hquadrature(f, 0.0, 1.0)

(1.0, 1.1102230246251565e-14)

Cubature también permite calcular integrales múltiples. En particular, permite integrar una función \(f(x)\) con \(x \in R^n\) sobre una caja multidimensional en la cual cada coordenada \(x_i\) está definida en un intervalos \([x_{i,min},x_{i,max}]\). La sintaxis es:

(int,err) = hcubature(f, xmin, xmax)
(int,err) = pcubature(f, xmin, xmax)

Los argumentos son:

f es una función del vector x y retorna un valor real..
xmin and xmax son Arreglos o Tuplas (o cualquier iterable) que especifican los límites xmin[i] y xmax[i] para todas las coordenadas.

Nota: Al igual que antes, Cubature solo permite integrar bajo límites finitos.

Ejemplo:

\[\int^{1}_{0} \int^{1}_{0} f(x,y) dx dy\]

donde: \[f(x,y) = x^3 y\]

funxy(x) =  x[1]^3*x[2]

funxy (generic function with 1 method)

intm, errm = hcubature(funxy, [0,0],[1,1])

(0.12499999999999999, 4.163336342344337e-17)