Integración adaptativa¶
Comenzamos con la regla del trapecio:
$$ T(f, a, b) = \frac{b-a}{2}(f(a) + f(b)) $$T(f, a, b) = (b-a)(f(a) + f(b))/2
Objetivo: integrar una función f en un intervalo [a,b], con un error estimado menor que tol, haciendo el mínimo de cálculos.
Llamamos al resultado
TA(f, a, b, tol): 'integral del trapecio adaptativa de la funciónfen el intervalo [a,b] con toleranciatol'.
Mecanismo:
- Calcula
I1=T(f,a,b)y por otro lado la integración según la regla del trapecio compuestaI2=T(f, a, m) + T(f, m, b), para el punto mediom = (a+b)/2. - Los errores que cometen
I1eI2es del tipo $$ \int_a^b f(x)\:dx = I_1 + \frac{b-a}{12}h^2f''(\xi_1) $$ $$ \int_a^b f(x)\:dx = I_2 + \frac{b-a}{12} \left(\frac{h}{2}\right)^2 f''(\xi_2) $$ Si asumimos que $f''(\xi_1)\approx f''(\xi_2)$, podemos estimar el error: $$ \int_a^b f(x)\:dx - I_2 \approx \frac{1}{3}(I_2 - I_1) $$ - Si esa estimación del error es menor que la tolerancia, reportamos la integral
I2=> TA(f, a, b, tol) = I2 - Si esa estimación del error es mayor que la tolerancia, repetimos el proceso en cada subintervalo $[a,m]$ y $[m,b]$. Atención: nuestra tolerancia para toda la integral es
tol, así que ahora le permitimos a cada subintervalo sólo la mitad de la tolerancia => TA(f, a, b, tol) = TA(f, a, m, tol/2) + TA(f, m, b, tol/2)
Observación: esta definición es recursiva.
Implementación¶
function I = T(f,a,b)
%Regla del trapecio simple
%
%Aproxima la integral de f en el intervalo [a,b] por (b-a)*(f(a) + f(b))/2
%
%El primer argumento es una funcion de variable real a valores reales
I = (b-a)*(f(a) + f(b))/2;
end
function I = TA(f, a, b, tol)
%Regla del trapecio adaptativa
%
% INPUTS:
% - f: integrando
% - a,b: extremos del intervalo de integracion
% - tol: tolerancia para la regla de parada. El algoritmo se detiene cuando
% la estimacion del error es menor que la tolerancia.
I1 = T(f, a, b);
m = (a+b)/2;
I2 = T(f, a, m) + T(f, m, b);
%estimación del error
E = (1/3)*abs(I1 - I2);
if E<tol
I = I2;
else
%llamada recursiva
I = TA(f, a, m, tol/2) + TA(f, m, b, tol/2);
end
end
quad(f,a,b)
f = @(x)(sin(x));
a = 0;
b = pi;
%Evaluamos la integral con mucha precision usando el metodo
% quad de scipy.integrate
[i_quad, e] = quad(f,a,b)
%Imprimimos la integral y el error cometido para distintas tolerancias
% comparamos el error con la integral calculada con mucha precision
tols = zeros(10);
for n = [1:10]
tols(n) = 10^(-(n-1));
end
for n = [1:10]
printf('-----\n')
tol = tols(n)
i_adap = TA(f, a, b, tol)
error_integral = abs(i_adap-i_quad)
end
Ejercicio¶
- Mira en los apuntes, un libro, o en internet cuál es la estimación del error cuando usamos el método de Simpson adaptativo en vez del método del trapecio adaptativo. Implementa el método de Simpson adaptativo, y comprueba el resultado con el obtenido con la regla del trapecio adaptativa.
Comparación de métodos¶
Para poder comparar los métodos del trapecio y de Simpson adaptativos, tenemos que medir cuántas iteraciones realiza cada uno.
El código siguiente no sólo devuelve la estimación de la integral, sino que lleva la cuenta de cuántas iteraciones necesitamos para alcanzar la tolerancia pedida.
Nota: En realidad el código anterior evalúa la función
fmás veces de las imprescindibles (varias veces en cada punto). La función siguiente no mide el número real de llamadas af, sino sólo el número que necesitaría el método del trapecio con una implementación más cuidadosa.No es difícil arreglar este problema, pero he preferido mantener el código simple. Si te interesa cómo solucionarlo con python, puedes mirar el archivo "visualizacion_trapecio.py". Para solucionarlo en matlab, pregunta al profesor.
function [I, E, n] = TA_completo(f, a, b, tol)
%Regla del trapecio adaptativa
%
% I, E, niters = TA_completo(f, a, b, tol)
%
% INPUTS:
% - f: integrando
% - a,b: extremos del intervalo de integracion
% - tol: tolerancia para la regla de parada. El algoritmo se detiene cuando
% la estimacion del error es menor que la tolerancia.
% OUTPUTS:
% - I: aproximación de la integral
% - E: estimación del error
% - niters: número de puntos donde se ha evaluado la función
% (sin contar los extremos del intervalo)
I1 = T(f, a, b);
m = (a+b)/2;
I2 = T(f, a, m) + T(f, m, b);
%estimación del error
E = (1/3)*abs(I1 - I2);
if E<tol
I = I2;
n = 1;
else
%llamada recursiva
[I_am, error_am, niters_am] = TA_completo(f, a, m, tol/2);
[I_mb, error_mb, niters_mb] = TA_completo(f, m, b, tol/2);
I = I_am + I_mb;
E = error_am + error_mb;
n = niters_am + niters_mb;
end
end
f = @(x)(sin(x));
a = 0;
b = pi;
%Evaluamos la integral con mucha precision usando el metodo
% quad de scipy.integrate
[i_quad, e] = quad(f,a,b)
%Imprimimos la integral y el error cometido para distintas tolerancias
% comparamos el error con la integral calculada con mucha precision
tols = zeros(10);
for n = [1:10]
tols(n) = 10^(-(n-1));
end
for n = [1:8]
printf('-----\n')
tol = tols(n)
[i_adap, estimacion_error, niters] = TA_completo(f, a, b, tol)
error_integral = abs(i_adap-i_quad)
end
Ejercicio¶
- Crea una gráfica de tipo
loglogdonde en el ejexmuestres la tolerancia y en el ejeyel número de iteraciones que necesita el método del trapecio adaptativo para alcanzar esa tolerancia. - Adapta esa función al método de Simpson, y compara el número de iteraciones que cada método necesita para obtener la misma tolerancia, en la misma gráfica
loglog.
Ejercicio¶
- Aprenda a usar el método
quad, ointegralcon tolerancia especificada de antemano, y a leer el número de evaluaciones de la función f. Compare el resultado con el número de evaluaciones de la regla del trapecio adaptativo. - Lea la documentación: ¿qué método numérico utiliza?
Ejercicio¶
Aprenda a usar
dblquadointegral2para hacer integrales de funciones de dos variables sobre regiones del planoEvalúe numéricamente la integral doble, para tres valores del radio $R=1,2,3$. $$ \iint_{\{x^2+y^2\leq R\}} \frac {1}{2\pi} e^{\frac{-1}{2}(x^2+y^2)}\:d x $$
Ejercicio¶
- Busque un método numérico para hacer integrales de funciones de
nvariables sobre regiones de $\mathbb{R}^n$. - Evalúe numéricamente la integral triple: $$ \iiint_{\{x^2+y^2+z^2\leq R\}} \sqrt{x^2+y^2+z^2}\:d x $$