Error en la regla del trapecio

La regla del trapecio simple comete el siguiente error: $$ \int_{x_0}^{x_1}f(x)\,dx=\frac{h}{2}\left[f(x_0)+f(x_1)\right]-\frac{h^3}{12}f''(\xi) $$ para $h=x_1-x_0$.

Vamos a probar esta fórmula para la función: $$ f(x)=\frac{1}{1+25x^2}. $$

Para vuestra conveniencia, adjuntamos la función f, sus derivadas y su primitiva en código python

f = lambda x:1/(1+25*x**2)
# Primera derivada de f
fp = lambda x:-50*x/(25*x**2 + 1)**2
# Segunda derivada de f
fp2 = lambda x:5000*x**2/(25*x**2 + 1)**3 - 50/(25*x**2 + 1)**2
# Primitiva de f
fi = lambda x:np.arctan(5*x)/5

y en código matlab

f = @(x)(1./(1+25*x.^2));
% Primera derivada de f
fp = @ (x)(-50*x./(25*x.^2 + 1).^2);
% Segunda derivada de f
fp2 = @ (x)(5000*x.^2./(25*x.^2 + 1).^3 - 50./(25*x.^2 + 1).^2);
% Primitiva de f
fi = @ (x)(atan(5*x)/5);

Apartado 1

Fijamos el punto $x_0=0$.

Queremos comparar la integral exacta $$ I(h) = \int_{x0-h}^{x0+h} f(x)\: d x $$ con su aproximación por la regla del trapecio simple, para distintos valores de h, y comprobar la validez del término de error $\frac{h^3}{12}f''(\xi)$.

  • Genere una lista de valores de h que decrece exponencialmente, por ejemplo hs = [1/2, 1/4, 1/8, 1/16 ...].
  • Para cada valor de h, calcule la diferencia entre la integral exacta $I(h)$ y su aproximación por la regla del trapecio.
  • Lo lamento mucho, en el enunciado del examen la función f aparecía definida como f = @(x)(1/(1+25*x.^2)); en vez de f = @(x)(1./(1+25*x.^2));.
In [7]:
f = @(x)(1./(1+25*x.^2));
% Primera derivada de f
fp = @ (x)(-50*x./(25*x.^2 + 1).^2);
% Segunda derivada de f
fp2 = @ (x)(5000*x.^2./(25*x.^2 + 1).^3 - 50./(25*x.^2 + 1).^2);
% Primitiva de f
fi = @ (x)(atan(5*x)/5);
In [8]:
x0 = 0;

%Usamos valores de h de la forma 1/2^n
Nhs = 20;
hs = 2.^(-[1:Nhs]);
trapecio = 0.5*(f(x0+hs) + f(x0-hs)).*(2*hs);
integral = fi(x0+hs) - fi(x0-hs);
errores_trap = abs(trapecio - integral);

title('Error en la regla del trapecio como funcion de h')
loglog(hs, errores_trap, 'g.-')
xlabel("h")
ylabel("error")
legend('error en la regla del trapecio')

Apartado 2

  • Muestre en una gráfica loglog el error cometido como función del valor de h, para x0 fijo.
  • Muestre en esa misma gráfica la expresión del error de truncamiento $|\frac{h^3}{12}f''(\xi)|$, pero reemplazando el punto desconocido $\xi\in(x_0-h,x_0+h)$ por $x_0$.
  • ¿Qué clase de comportamiento tiene el error de redondeo?
  • Comprobamos que la estimación del error no es muy buena cuando h es grande, pero es muy ajustada cuando h es pequeño. Esta gráfica contrasta con la gráfica similar que hicimos para aproximar las derivadas por métodos de diferencias: al aproximar derivadas aparecían errores de redondeo que hacían que el error aumentase al disminuir h, pero aquí la estimación del error de truncamiento aproxima bastante bien al error real sin necesidad de incorporar el error de redondeo en el análisis.
In [9]:
%Usamos valores de h de la forma 1/2^n
Nhs = 20;
hs = 2.^(-[1:Nhs]);
trapecio = 0.5*(f(x0+hs) + f(x0-hs)).*(2*hs);
integral = fi(x0+hs) - fi(x0-hs);
errores_trap = abs(trapecio - integral);

estimacion_error = abs((2*hs).^3*fp2(x0)/12);
set(gca, 'XScale', 'log', 'YScale', 'log');
hold on;
loglog(hs, errores_trap, 'b.-')
loglog(hs, estimacion_error, 'g-')
hold off;

title('Error en la regla del trapecio como funcion de h')
xlabel("h")
ylabel("error")
legend('error en la regla del trapecio', 'estimacion del error')

Apartado 3

Ahora fijamos h=1/10, y comparamos el error cuando cambiamos el punto $x_0$

  • Muestre en una gráfica con escala lineal (lo habitual) el error cometido al aproximar la integral exacta $I(x_0) = \int_{x0-h}^{x0+h} f(x)\: d x$ con su aproximación por la regla del trapecio simple como función del valor de x0, para h=10/10 fijo.
  • Muestre en esa misma gráfica la expresión del error $\frac{h^3}{12}f''(\xi)$, pero reemplazando el punto desconocido $\xi\in(x_0-h,x_0+h)$ por $x_0$.
  • Comprobamos de nuevo que la estimacion del error se parece mucho al error real.
In [10]:
xs = linspace(0,2*pi,100);
fpps = fp2(xs);
h = 0.1;

fi_trapecio = (f(xs + h) + f(xs-h))*(2*h)/2;
fp_integrales = fi(xs+h) - fi(xs-h);

errores_trapecio = fp_integrales - fi_trapecio;
cota_error_trapecio = - ((2*h)^3/12).*fp2(xs);

hold on;
title('Error en la aproximacion de la integral con la regla del trapecio como funcion de x')
plot(xs, errores_trapecio, 'g.')
plot(xs, cota_error_trapecio, 'b-')
xlabel("x")
ylabel("error")
legends = cell(1,2);
legends{1} = 'error en la formula del trapecio';
legends{2} = 'estimacion del error';
legend(legends)
hold off;
  • Probamos la función sin, cuya seguna derivada es -sin.
In [11]:
f = @(x)(sin(x));
fp = @(x)(cos(x));
fp2 = @(x)(-sin(x));
fpi = @(x)(-cos(x));

xs = linspace(0,2*pi,100);
fpps = fp2(xs);
h = 0.1;

fi_trapecio = (f(xs + h) + f(xs-h))*(2*h)/2;
fp_integrales = fpi(xs+h) - fpi(xs-h);

errores_trapecio = fp_integrales - fi_trapecio;
cota_error_trapecio = - ((2*h)^3/12)*fp2(xs);

hold on;
title('Error en la aproximacion de la integral con la regla del trapecio como funcion de x')
plot(xs, errores_trapecio, 'g.')
plot(xs, cota_error_trapecio, 'b-')
xlabel("x")
ylabel("error")
legends = cell(1,2);
legends{1} = 'error en la formula del trapecio';
legends{2} = 'estimacion del error';
legend(legends)
hold off;
In [ ]: