Ejercicios APROXIMACIÓN DE FUNCIONES POR POLINOMIOS

Colección de ejercicios prácticos

Fórmula de cuadratura de Lagrange para \(e^{-2x}\) en \([0,1]\)

Problema:

Consideremos la función \(f(x) = e^{-2x}\), el intervalo \([0,1]\) y los nodos: \[ x_0 = 0,\quad x_1 = \tfrac{1}{3},\quad x_2 = \tfrac{2}{3},\quad x_3 = 1. \] Calcular la fórmula de cuadratura a partir de los polinomios base de Lagrange.

Solución:

Para obtener la fórmula de cuadratura de Newton–Cotes basada en los cuatro nodos \(\{0,\;1/3,\;2/3,\;1\}\), construimos los polinomios base de Lagrange \(\{L_0(x),L_1(x),L_2(x),L_3(x)\}\), y luego integramos cada uno sobre \([0,1]\). La aproximación resultante para \(\displaystyle \int_{0}^{1} f(x)\,dx\) es \[ \int_{0}^{1} f(x)\,dx \;\approx\; \sum_{i=0}^{3} w_i\,f(x_i), \] donde \[ w_i \;=\; \int_{0}^{1} L_i(x)\,dx, \quad L_i(x) = \prod_{\substack{j=0 \\ j \neq i}}^{3} \frac{x - x_j}{x_i - x_j}. \] A continuación se detallan cada uno de los pasos:

1. Polinomios base de Lagrange

Con los nodos \[ x_0 = 0,\quad x_1 = \tfrac{1}{3},\quad x_2 = \tfrac{2}{3},\quad x_3 = 1, \] definimos:

  • \[ L_0(x) = \frac{(x - x_1)(x - x_2)(x - x_3)} {(x_0 - x_1)(x_0 - x_2)(x_0 - x_3)} = \frac{\bigl(x - \tfrac{1}{3}\bigr)\bigl(x - \tfrac{2}{3}\bigr)(x - 1)} {\bigl(0 - \tfrac{1}{3}\bigr)\bigl(0 - \tfrac{2}{3}\bigr)\bigl(0 - 1\bigr)}. \]
  • \[ L_1(x) = \frac{(x - x_0)(x - x_2)(x - x_3)} {(x_1 - x_0)(x_1 - x_2)(x_1 - x_3)} = \frac{x\,(x - \tfrac{2}{3})(x - 1)} {\bigl(\tfrac{1}{3} - 0\bigr)\bigl(\tfrac{1}{3} - \tfrac{2}{3}\bigr)\bigl(\tfrac{1}{3} - 1\bigr)}. \]
  • \[ L_2(x) = \frac{(x - x_0)(x - x_1)(x - x_3)} {(x_2 - x_0)(x_2 - x_1)(x_2 - x_3)} = \frac{x\,\bigl(x - \tfrac{1}{3}\bigr)(x - 1)} {\bigl(\tfrac{2}{3} - 0\bigr)\bigl(\tfrac{2}{3} - \tfrac{1}{3}\bigr)\bigl(\tfrac{2}{3} - 1\bigr)}. \]
  • \[ L_3(x) = \frac{(x - x_0)(x - x_1)(x - x_2)} {(x_3 - x_0)(x_3 - x_1)(x_3 - x_2)} = \frac{x\,\bigl(x - \tfrac{1}{3}\bigr)\bigl(x - \tfrac{2}{3}\bigr)} {\bigl(1 - 0\bigr)\bigl(1 - \tfrac{1}{3}\bigr)\bigl(1 - \tfrac{2}{3}\bigr)}. \]

2. Cálculo de los pesos \(w_i\)

Cada peso se obtiene integrando \(L_i(x)\) en \([0,1]\): \[ w_i = \int_{0}^{1} L_i(x)\,dx. \] Al realizar la integración (por ejemplo, con cálculo simbólico o a mano), se obtienen los valores:

  • \(\displaystyle w_0 = \int_{0}^{1} L_0(x)\,dx = \tfrac{1}{8}\).
  • \(\displaystyle w_1 = \int_{0}^{1} L_1(x)\,dx = \tfrac{3}{8}\).
  • \(\displaystyle w_2 = \int_{0}^{1} L_2(x)\,dx = \tfrac{3}{8}\).
  • \(\displaystyle w_3 = \int_{0}^{1} L_3(x)\,dx = \tfrac{1}{8}\).

Por lo tanto, la fórmula de cuadratura resultante es

\[ \int_{0}^{1} e^{-2x}\,dx \;\approx\; \frac{1}{8}\,f(0) \;+\; \frac{3}{8}\,f\Bigl(\tfrac{1}{3}\Bigr) \;+\; \frac{3}{8}\,f\Bigl(\tfrac{2}{3}\Bigr) \;+\; \frac{1}{8}\,f(1), \] donde \(f(x) = e^{-2x}\).

3. Expresión explícita de la cuadratura

Reemplazando \(f(x) = e^{-2x}\) en los nodos, obtenemos:

\[ \int_{0}^{1} e^{-2x}\,dx \;\approx\; \frac{1}{8}\,e^{0} \;+\; \frac{3}{8}\,e^{-\,\tfrac{2}{3}} \;+\; \frac{3}{8}\,e^{-\,\tfrac{4}{3}} \;+\; \frac{1}{8}\,e^{-2}. \]

Esta es la fórmula de cuadratura basada en los polinomios de Lagrange para los cuatro nodos en \([0,1]\). Se ha obtenido la fórmula de Newton Cotes para \(n=3\), la regla 3/8.

Interpolación en cuadratura: comprobación de fórmulas

Problema:

Comprobar si las siguientes fórmulas de cuadratura son de tipo interpolatorio:

  1. \( \int_{0}^{1} f(x)\,dx \;\approx\; f(0) \;+\; f\Bigl(\tfrac{2}{3}\Bigr) \)
  2. \( \int_{0}^{1} f(x)\,dx \;\approx\; f(0) \;+\; f\Bigl(\tfrac{1}{3}\Bigr) \;+\; f\Bigl(\tfrac{2}{3}\Bigr) \;+\; f(1) \)

Solución:

(a) La fórmula \(\displaystyle \int_{0}^{1}f(x)\,dx \approx f(0) + f\bigl(\tfrac{2}{3}\bigr)\) no es interpolatoria.

Para que una fórmula de cuadratura con dos nodos \(\{x_0=0,\;x_1=\tfrac{2}{3}\}\) sea interpolatoria de orden 1, debe existir un par de pesos \(\{w_0,\,w_1\}\) tal que \[ \int_{0}^{1} p(x)\,dx \;=\; w_0\,p(0)\;+\; w_1\,p\Bigl(\tfrac{2}{3}\Bigr) \] sea exacto para todo polinomio \(p(x)\) de grado \(\le 1\). Es decir, basta verificar la exactitud en \(p(x)=1\) y en \(p(x)=x\).

1. Exactitud para \(p(x)=1\):

  • \(\displaystyle \int_{0}^{1} 1\,dx = 1.\)
  • El lado derecho con pesos \(w_0,\,w_1\) debe ser \[ w_0\cdot 1 \;+\; w_1\cdot 1 \;=\; w_0 + w_1 \;=\; 1. \]

2. Exactitud para \(p(x)=x\):

  • \(\displaystyle \int_{0}^{1} x\,dx = \tfrac{1}{2}.\)
  • El lado derecho impondría \[ w_0\cdot 0 \;+\; w_1\cdot \tfrac{2}{3} \;=\; \tfrac{2}{3}\,w_1 \;=\; \tfrac{1}{2}. \] De aquí \(w_1 = \tfrac{3}{4}.\)
  • Entonces de \(w_0 + w_1 =1\) se deduce \(w_0 = 1 - \tfrac{3}{4} = \tfrac{1}{4}.\)

Por tanto, la única combinación de pesos que haría exacta la fórmula para polinomios de grado 1 sería \[ w_0 = \tfrac{1}{4}, \quad w_1 = \tfrac{3}{4}, \] y la fórmula interpolatoria correspondiente sería \[ \int_{0}^{1} f(x)\,dx \;\approx\; \tfrac{1}{4}\,f(0) \;+\; \tfrac{3}{4}\,f\Bigl(\tfrac{2}{3}\Bigr) \quad\text{(exacto para todo polinomio de grado \(\le1\)).} \] Sin embargo, la fórmula dada en el enunciado usa coeficientes “1” delante de \(f(0)\) y \(f(2/3)\), es decir, \[ f(0) \;+\; f\Bigl(\tfrac{2}{3}\Bigr), \] que correspondería a pesos \(\tilde w_0 = 1,\;\tilde w_1 = 1\). Pero \((\tilde w_0,\tilde w_1) = (1,1)\) no satisface las ecuaciones necesarias \(\tilde w_0 + \tilde w_1 = 1\) y \(\tfrac{2}{3}\,\tilde w_1 = \tfrac12\). Por ello la fórmula \(\displaystyle \int_{0}^{1}f(x)\,dx \approx f(0)+f(2/3)\) no es de tipo interpolatorio.

(b) La fórmula \(\displaystyle \int_{0}^{1} f(x)\,dx \approx f(0) + f\bigl(\tfrac{1}{3}\bigr) + f\bigl(\tfrac{2}{3}\bigr) + f(1)\) no es de tipo interpolatorio.

Consideramos ahora los cuatro nodos \(\;x_0 = 0,\;x_1 = \tfrac{1}{3},\;x_2 = \tfrac{2}{3},\;x_3 = 1\). Una fórmula de cuadratura “de tipo interpolatorio” en estos cuatro nodos debe tener la forma \[ \int_{0}^{1} f(x)\,dx \;\approx\; w_0\,f(0) \;+\; w_1\,f\Bigl(\tfrac{1}{3}\Bigr) \;+\; w_2\,f\Bigl(\tfrac{2}{3}\Bigr) \;+\; w_3\,f(1), \] y ser exacta para todos los polinomios de grado \(\le 3\). Para determinar los pesos \(\{w_0,w_1,w_2,w_3\}\), imponemos exactitud en la base \(\{\,1,\,x,\,x^2,\,x^3\}\).

1. Para \(p(x)=1\):

  • \(\displaystyle \int_{0}^{1} 1\,dx = 1.\)
  • La suma de pesos debe satisfacer \[ w_0 + w_1 + w_2 + w_3 \;=\; 1. \]

2. Para \(p(x)=x\):

  • \(\displaystyle \int_{0}^{1} x\,dx = \tfrac{1}{2}.\)
  • Condición en los pesos: \[ w_0 \cdot 0 \;+\; w_1 \cdot \bigl(\tfrac{1}{3}\bigr) \;+\; w_2 \cdot \bigl(\tfrac{2}{3}\bigr) \;+\; w_3 \cdot 1 \;=\; \tfrac{1}{2}. \]

3. Para \(p(x)=x^2\):

  • \(\displaystyle \int_{0}^{1} x^2\,dx = \tfrac{1}{3}.\)
  • Condición: \[ w_0 \cdot 0^2 \;+\; w_1 \cdot \Bigl(\tfrac{1}{3}\Bigr)^2 \;+\; w_2 \cdot \Bigl(\tfrac{2}{3}\Bigr)^2 \;+\; w_3 \cdot 1^2 \;=\; \tfrac{1}{3}. \]

4. Para \(p(x)=x^3\):

  • \(\displaystyle \int_{0}^{1} x^3\,dx = \tfrac{1}{4}.\)
  • Condición: \[ w_0 \cdot 0^3 \;+\; w_1 \cdot \Bigl(\tfrac{1}{3}\Bigr)^3 \;+\; w_2 \cdot \Bigl(\tfrac{2}{3}\Bigr)^3 \;+\; w_3 \cdot 1^3 \;=\; \tfrac{1}{4}. \]

El sistema lineal resultante en las incógnitas \(w_0,w_1,w_2,w_3\) es:

\[ \begin{pmatrix} 1 & 1 & 1 & 1 \\ 0 & \tfrac{1}{3} & \tfrac{2}{3} & 1 \\ 0 & \bigl(\tfrac{1}{3}\bigr)^2 & \bigl(\tfrac{2}{3}\bigr)^2 & 1 \\ 0 & \bigl(\tfrac{1}{3}\bigr)^3 & \bigl(\tfrac{2}{3}\bigr)^3 & 1 \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ w_2 \\ w_3 \end{pmatrix} = \begin{pmatrix} 1 \\[6pt] \tfrac{1}{2} \\[6pt] \tfrac{1}{3} \\[6pt] \tfrac{1}{4} \end{pmatrix}. \]

Resolviendo este sistema (por ejemplo en MATLAB o a mano), se obtiene el vector de pesos:

\[ w_0 = \tfrac{1}{8}, \quad w_1 = \tfrac{3}{8}, \quad w_2 = \tfrac{3}{8}, \quad w_3 = \tfrac{1}{8}. \]

Por tanto, la fórmula interpolatoria exacta para todos los polinomios de grado \(\le3\) es

\[ \int_{0}^{1} f(x)\,dx \;\approx\; \frac{1}{8}\,f(0) \;+\; \frac{3}{8}\,f\Bigl(\tfrac{1}{3}\Bigr) \;+\; \frac{3}{8}\,f\Bigl(\tfrac{2}{3}\Bigr) \;+\; \frac{1}{8}\,f(1). \]

Observamos que, para que sea de tipo interpolatorio es necesario incluir sus correspondientes pesos \(\{1/8,\,3/8,\,3/8,\,1/8\}\).

Conclusión:

  • En el caso (a), no existen pesos \((1,1)\) que satisfagan las condiciones de exactitud para polinomios de grado \(\le1\). Por tanto, la fórmula dada \(\int_{0}^{1}f(x)\approx f(0)+f(2/3)\) no es interpolatoria.
  • En el caso (b), sí podemos encontrar pesos \(\bigl(1/8,\,3/8,\,3/8,\,1/8\bigr)\) tales que la regla \(\displaystyle \int_{0}^{1} f(x)\,dx = \sum_{i=0}^{3} w_i\,f(x_i)\) sea exacta para todos los polinomios de grado \(\le3\). Por tanto, dicha regla “con pesos \(\tfrac{1}{8},\tfrac{3}{8},\tfrac{3}{8},\tfrac{1}{8}\)” sí es de tipo interpolatorio (sería la regla 3/8).

Aproximación de \(\displaystyle \int_{0}^{\pi} e^{\cos(x)}\,\sin(x)\,dx\)

Problema:

Obtener con MATLAB la aproximación de la integral \[ \int_{0}^{\pi} e^{\cos(x)} \,\sin(x)\,dx \] utilizando:

  1. La regla de los tres octavos (\(n=3\)).
  2. La regla de Hardy (Newton–Cotes cerrado con \(n=6\)).

Solución:

Observación previa: Notemos que la integral puede evaluarse exactamente mediante la sustitución \(u = \cos(x)\), ya que \(du = -\sin(x)\,dx\). Entonces \[ \int_{0}^{\pi} e^{\cos(x)}\,\sin(x)\,dx = \int_{\cos(0)=1}^{\cos(\pi)=-1} e^{u}\,(-\,du) = \int_{-1}^{1} e^{u}\,du = e^{1} - e^{-1} \approx 2.3504. \] Sin embargo, aquí interesa practicar fórmulas de Newton–Cotes, por lo que las aproximaremos numéricamente.

(a) Regla de los tres octavos (\(n=3\))

La regla de los tres octavos es un Newton–Cotes cerrado de grado 3 que utiliza cuatro nodos equiespaciados en \([0,\pi]\). Con \(n=3\) subintervalos, el paso es \[ h = \frac{\pi - 0}{3} = \frac{\pi}{3}. \] Los nodos son: \[ x_0 = 0,\quad x_1 = \frac{\pi}{3},\quad x_2 = \frac{2\pi}{3},\quad x_3 = \pi. \] La fórmula de los tres octavos es \[ \int_{0}^{\pi} f(x)\,dx \;\approx\; \frac{\pi}{8}\Bigl[\,f(x_0) \;+\; 3\,f(x_1) \;+\; 3\,f(x_2) \;+\; f(x_3)\Bigr], \] donde \(f(x) = e^{\cos(x)}\,\sin(x)\). Sustituyendo queda \[ \int_{0}^{\pi} f(x)\,dx \;\approx\; \frac{\pi}{8}\Bigl[\,f(0) + 3\,f\bigl(\tfrac{\pi}{3}\bigr) + 3\,f\bigl(\tfrac{2\pi}{3}\bigr) + f(\pi)\Bigr]. \]

En MATLAB, se implementa así:


% ------------------------------------------------------------
% Regla de los tres octavos
% ------------------------------------------------------------

% Definir la función
f = @(x) exp(cos(x)) .* sin(x);

% Número de nodos n+1 con n=3
n = 3;a=0;b=pi;
h = (b-a) / n;

% Nodos
x = 0:h:pi;

% Evaluar f en los nodos
fx = f(x);

% Fórmula de los tres octavos
I38 = (pi/8) * (fx(1) + 3*fx(2) + 3*fx(3) + fx(4));

% También
a=0;b=pi;
h = (b-a) / n;
nodos=a:h:b;
pesos=[1 3 3 1]/8;
valor=(b-a)*sum(f(nodos).*pesos)

% Mostrar resultado
fprintf('Regla de los tres octavos: %.6f\n', I38);

Al ejecutar el código anterior en MATLAB, se obtiene algo como:


Regla de los tres octavos:  2.300948166133495

Comparado con el valor exacto \(e - e^{-1} \approx 2.350402387287603\), el error absoluto es aproximadamente \(\approx 0.049454221154108\).

(b) Regla de Hardy (Newton–Cotes cerrado con \(n=6\))

La llamada “regla de Hardy” se corresponde con la fórmula de Newton–Cotes cerrada de orden 6 (es decir, 7 nodos). En \([0,\pi]\), tomamos \(n=6\) subintervalos de longitud \[ h = \frac{\pi-0}{6} = \frac{\pi}{6}. \] Los nodos son: \[ x_i = i\,h = i\,\frac{\pi}{6}, \quad i = 0,1,\dots,6. \] El peso asociado a cada nodo viene de la fórmula Newton–Cotes de orden 6 (7 puntos), cuyos coeficientes son: \[ [\,w_0,\,w_1,\,w_2,\,w_3,\,w_4,\,w_5,\,w_6\,] \;=\; \frac{1}{840}\bigl[\,41,\,216,\,27,\,272,\,27,\,216,\,41\,\bigr]. \] Por tanto: \[ \int_{0}^{\pi} f(x)\,dx \;\approx\; \frac{\pi}{840}\Bigl[ 41\,f(x_0)\;+\;216\,f(x_1)\;+\;27\,f(x_2)\;+\;272\,f(x_3) \;+\;27\,f(x_4)\;+\;216\,f(x_5)\;+\;41\,f(x_6) \Bigr]. \]

En MATLAB, se implementa de la siguiente manera:


% ------------------------------------------------------------
% Regla de Hardy (Newton–Cotes de grado 6)
% ------------------------------------------------------------

% Definir la función
f = @(x) exp(cos(x)) .* sin(x);

% Número de subintervalos n=6
n = 6;
h = pi / n;

% Nodos x_i = i*h, i=0..6
x = 0:h:pi;    % vector columna de 7 nodos: [0; π/6; 2π/6; ...; 6π/6]

% Pesos de Newton–Cotes cerrado de orden 6
w = [41, 216, 27, 272, 27, 216, 41]/840;   % coeficientes

% Evaluar f en los 7 nodos
fx = f(x);

% Aproximación con la regla de Hardy
IH = pi*sum(w .* fx);

% Mostrar resultado
fprintf('Regla de Hardy (grado 6): %.6f\n', IH);

Al ejecutar este bloque en MATLAB, se obtiene aproximadamente:


Regla de Hardy (grado 6): 2.344693

Comparado con el valor exacto \(e - e^{-1} \approx 2.350402387287603\), el error absoluto es aproximadamente \(2.35040238728760 - 2.344693 \approx 0.005709794172118\).

4. Resumen de resultados:

Regla Valor aproximado Valor exacto Error absoluto
Tres octavos (n=3) 2.300948166133495 2.350402387287603 ≈ .049454221154108
Hardy (Newton–Cotes n=6) 2.344693 2.350402387287603 ≈ .005709794172118

Se aprecia que la regla de Hardy (grado 6) proporciona una aproximación más precisa que la regla de los tres octavos (grado 3), dado el mismo intervalo \([0,\pi]\).

Trapecios compuestos y estimación de subintervalos

Problema:

Utiliza la fórmula de los trapecios compuesta dividiendo el intervalo \([0, \pi]\) en tres subintervalos para aproximar la integral \[ \int_{0}^{\pi} \sin(x)\,dx. \] Además, calcula cuántos subintervalos serían necesarios para que el error cumpla \[ E \;\le\; 0.5 \times 10^{-8}. \]

Solución:

1. Aproximación con tres subintervalos considerando la fórmula de los trapecios compuestos)

La regla de los trapecios compuesta en \([0,\pi]\) con \(N=3\) subintervalos divide el intervalo en tres partes iguales de longitud \[ h = \frac{\pi - 0}{N} = \frac{\pi}{3}. \] Los nodos son: \[ x_0 = 0, \quad x_1 = \frac{\pi}{3}, \quad x_2 = \frac{2\pi}{3}, \quad x_3 = \pi. \] La fórmula compuesta es \[ \int_{0}^{\pi} \sin(x)\,dx \;\approx\; \frac{h}{2}\Bigl[ f(x_0) + 2\,f(x_1) + 2\,f(x_2) + f(x_3) \Bigr], \] donde \(f(x) = \sin(x)\). Sustituyendo \(h = \pi/3\): \[ \int_{0}^{\pi} \sin(x)\,dx \;\approx\; \frac{\pi/3}{2}\bigl[\sin(0) \;+\; 2\,\sin(\tfrac{\pi}{3}) \;+\; 2\,\sin(\tfrac{2\pi}{3}) \;+\; \sin(\pi)\bigr]. \] Observemos que: \(\sin(0) = 0,\ \sin(\pi) = 0,\ \sin(\tfrac{\pi}{3}) = \tfrac{\sqrt{3}}{2},\ \sin(\tfrac{2\pi}{3}) = \tfrac{\sqrt{3}}{2}.\) Por tanto: \[ \int_{0}^{\pi} \sin(x)\,dx \;\approx\; \frac{\pi}{6}\Bigl[\,0 \;+\; 2\cdot\tfrac{\sqrt{3}}{2} \;+\; 2\cdot\tfrac{\sqrt{3}}{2} \;+\; 0\Bigr] \;=\; \frac{\pi}{6}\,\bigl[\,2\,\tfrac{\sqrt{3}}{2} + 2\,\tfrac{\sqrt{3}}{2}\bigr] \;=\; \frac{\pi}{6}\,\bigl[\,\sqrt{3} + \sqrt{3}\bigr] \;=\; \frac{\pi}{6}\cdot 2\sqrt{3} \;=\; \frac{\pi\,\sqrt{3}}{3}. \] Numéricamente, \(\displaystyle \frac{\pi\,\sqrt{3}}{3} \approx 1.8137993642.\) El valor exacto de \(\int_{0}^{\pi} \sin(x)\,dx\) es \(2\). Por lo tanto, el error con \(N=3\) es \[ E_{N=3} \approx \Bigl\lvert 2 - 1.8137993642 \Bigr\rvert \approx 0.1862006358. \] 


% ------------------------------------------------------------
% Trapecios compuestos con N=3
% ------------------------------------------------------------
f = @(x) sin(x);

N = 3;
h = pi / N; %Longitud de cada subintervalo, en este caso también longitud entre los nodos

x = 0:h:pi;
fx = f(x);

T3 = (h/2) * (fx(1) + 2*fx(2) + 2*fx(3) + fx(4));  %1.813799364234218

fprintf('Aproximación trapecios (N=3): %.10f\n', T3);
fprintf('Valor exacto: 2.0000000000\n');
fprintf('Error: %.10f\n', abs(2 - T3));

2. Estimación del número de subintervalos para \(E \le 0.5\times10^{-8}\)

Para la regla compuesta de los trapecios, el error absoluto está acotado por \[ \bigl\lvert E \bigr\rvert \;\le\; \frac{h^3 \cdot N}{12}\,\max_{x\in[a,b]} \bigl\lvert f''(x)\bigr\rvert, \] donde \([a,b] = [0,\pi]\) y \(h = \frac{b-a}{N}\). En nuestro caso, \(f(x) = \sin(x)\), por lo que \(f''(x) = -\,\sin(x)\) y \(\lvert f''(x)\rvert \le 1\) en \([0,\pi]\). Entonces, \[ \bigl\lvert E \bigr\rvert \;\le\; \frac{\pi^3}{12 N^2} \] Para que \(E \le 0.5\times10^{-8}\), basta exigir \[ \frac{\pi^3}{12\,N^2} \;\le\; 0.5 \times 10^{-8}. \] De aquí: \[ N^2 \;\ge\; \frac{\pi^3}{12 \cdot 0.5\times10^{-8}} \;=\; \frac{\pi^3}{6\times10^{-8}}. \] Por tanto, \[ N \;\ge\; \sqrt{\frac{\pi^3}{6\times10^{-8}}}. \] Numéricamente: \[ \sqrt{\frac{\pi^3}{6\times10^{-8}}} \;\approx\; 2.273260385\times10^{4} \;\approx\; 22733. \] Así, se necesitan al menos \(N = 22733\) subintervalos para garantizar \(E \le 0.5\times10^{-8}\).


% ------------------------------------------------------------
% Cálculo de N para que |E| <= 0.5e-8 en trapecios compuestos
% ------------------------------------------------------------

num = sqrt(pi^3 / (6e-8));
N = ceil(num);

fprintf('Número mínimo de subintervalos N: %d\n', N);
fprintf('Con N = %d, h = %.10e\n', N, pi/N);

Por lo tanto, se requieren \(N = 22733\) subintervalos para que la regla compuesta de los trapecios satisfaga el criterio \(E \le 0.5\times10^{-8}\).

Aproximación de \(\pi\) mediante \(\displaystyle \int_{0}^{1} \frac{4}{1+x^2}\,dx\)

Problema:

Aproximar el valor de \(\pi\) mediante la integral $$\int_0^1 {{4 \over {1 + {x^2}}}} dx = 4\,\,arctg\left( 1 \right) = \pi $$ utilizando las fórmulas compuestas del trapecio y de Simpson dividiendo el intervalo \([0,1]\) respectivamente en 1024 subintervalos y 512 (en ambos casos con 1025 nodos).

Solución:

Consideramos \[ f(x) \;=\; \frac{4}{1 + x^2}, \quad \int_{0}^{1} f(x)\,dx = \pi. \] Para aproximar esta integral dividiendo \([0,1]\) en \(N\) subintervalos, cada intervalo tendrá longigud \(h = \frac{1 - 0}{N} \). Para el caso de la regla de los trapecios, los nodos son \(N+1\) y están separados a distancia \(h\) y son \[ x_i = i\,h, \quad i = 0,1,\dots,N. \] En el caso de la regla de Simpson, si se divide el intervalo \([0,1]\) en \(N\) subintervalos, la distancia entre los nodos es \(h/2\) y son \(2N+1\) nodos.

Si queremos utilizar el mismo número de nodos en ambos casos, 1025, en la regla del trapecio y en la regla de Simpson, se considerará \(N=1024\) para el trapecio y \(N=512\) para la regla de Simpson.

A continuación presentamos las fórmulas compuestas y su implementación en MATLAB.

(a) Trapecio compuesto dividiendo el intervalo en \(N = 1024\) subintervalos

La regla del trapecio compuesta en \([0,1]\) con paso \(h\) es \[ \int_{0}^{1} f(x)\,dx \;\approx\; \frac{h}{2}\Bigl[ f(x_0) + 2\,\sum_{i=1}^{n-1} f(x_i) + f(x_n) \Bigr]. \] Sustituyendo \(N=1024\), \(h=\tfrac{1}{1024}\), y \(f(x)=4/(1+x^2)\). 


% ------------------------------------------------------------
% Trapecio compuesto con N=1024
% ------------------------------------------------------------
f = @(x) 4 ./ (1 + x.^2);
format long

N = 1024;                % número de subintervalos, en cada uno se aplica trapecio
h = 1 / N;               % Longitud de cada subintervalo, también longitud entre nodos

% Generamos los nodos x_i
x = linspace(0,1,N+1);          % vector columna de 1025 puntos: 0, h, 2h, ..., 1
fx = f(x);

% Trapecio compuesto
IT = (h/2) * (fx(1) + 2*sum(fx(2:end-1)) + fx(end));

fprintf('Trapecio compuesto (N=1024): %.10f\n', IT);

Al ejecutar este bloque en MATLAB se obtiene:


Trapecio compuesto (N=1024): 3.141592494644073

El valor exacto de \(\pi\) es aproximadamente \(3.1415926536\). El error absoluto con el trapecio compuesto es \(\approx -1.589457205852796e-07.\)

(b) Simpson compuesto dividiendo [0,1] en \(N = 512\) subintervalos

La longitud de cada intervalo es \(h = \frac{1}{512}\), la longitud entre dos nodos consecutivos \(h/2 =1/1024\) y el número de nodos es \(2N+1=1025\). La fórmula es \[ \int_{0}^{1} f(x)\,dx \;\approx\; \frac{h}{6}\Bigl[ f(x_0) + 4\,\sum_{\substack{i=1 \\ i\,\mathrm{impar}}}^{n-1} f(x_i) + 2\,\sum_{\substack{i=2 \\ i\,\mathrm{par}}}^{n-2} f(x_i) + f(x_n) \Bigr]. \]


% ------------------------------------------------------------
% Simpson compuesto con N=512
% ------------------------------------------------------------
f = @(x) 4 ./ (1 + x.^2);
format long

N = 512;                % número de subintervalos en los que se divide [0,1] para aplicar en cada uno simpson
h = 1 / N;

% Generar nodos
x = linspace(0,1,2*N+1);          % Nodos 2N+1
fx = f(x);

% Sumas para Simpson
sum1 = sum(fx(2:2:end-1));  % f(x_1), f(x_3), ..., f(x_{n-1}) (índices impares)
sum2 = sum(fx(3:2:end-2));  % f(x_2), f(x_4), ..., f(x_{n-2}) (índices pares)

IS = (h/6) * (fx(1) + 4*sum1+ 2*sum2 + fx(end));

fprintf('Simpson compuesto (N=512): %.10f\n', IS);

Al ejecutar este bloque en MATLAB se obtiene:


Simpson compuesto (N=512): 3.141592653589794

Comparado con \(\pi\approx 3.1415926536\), el error absoluto con Simpson compuesto es aproximadamente del orden de \(4.440892098500626e-16\).

4. Resumen de resultados:

Método Aproximación Valor numérico \(\pi\) Error absoluto
Trapecio compuesto (N=1024) 3.141592494644073 3.141592653589793 ≈ 1.58×10−7
Simpson compuesto (N=512) 3.141592653589794 3.141592653589793 ≈ 4.4×10−16 (numérico)

Se observa que la regla de Simpson compuesto brinda una precisión mayor que la regla del trapecio compuesto, usando el mismo número de nodos.

Cota de error para \(\displaystyle \int_{0}^{0.2} \sin(x) e^{x^2}\,dx\) usando regla del trapecio

Problema:

Calcular una cota del error al aproximar \[ \int_{0}^{0.2} \sin(x) e^{x^2}\,dx \] utilizando la regla del trapecio sobre el intervalo \([0,0.2]\).

Solución:

Sea \(f(x) = sin(x) e^{x^2}\). Para la regla del trapecio en un solo subintervalo \([a,b]\), con \(a=0\) y \(b=0.2\), el error exacto se puede expresar como \[ E_T = -\,\frac{(b - a)^3}{12}\,f''(\xi) \quad\text{para algún }\xi\in[a,b]. \] Dado que \(f''(x) = e^{x^2}\sin(x) + 4x^2 e^{x^2}\sin(x) + 4x e^{x^2}\cos(x)\) su valor absoluto está acotado por \(e^{x^2}(1+4x^2+4x)\). En \([0,0.2]\) se cumple \(\max_{x\in[0,0.2]} |f''(x)| = 3.4{e^{0.04}}\). Por tanto, una cota para el error absoluto es \[ |E_T| \;\le\; \frac{(b - a)^3}{12}\,\max_{x\in[0,0.2]} |f''(x)| \;\;\le\; \frac{(0.2)^3}{12}\,3.4 e^{0.04}. \] Calculamos numéricamente: \[ |E_T| \;\approx\; 0.0023591711. \] 


% ------------------------------------------------------------
% Cálculo de la cota de error del trapecio
% ------------------------------------------------------------
fmax = exp(0.2^2)*(1+4*0.2*2+4*0.2);           % máximo de f''(x) en [0,0.2]: 3.538756632254120
a = 0;
b = 0.2;
cota = ((b - a)^3 / 12) * fmax;

fprintf('Cota del error |E_T| <= %.10f\n', cota);

Cálculo del valor aproximado


f = @(x) sin(x).*exp(x.^2);
a = 0;
b = 0.2;
IT = (b-a)*(f(a)+f(b))/2;

% Valor numérico comando integral
I = integral(f,0,0.2)
errorAbsoluto=abs(I-IT)   % 3.407107678857976e-04

Resultado:

La cota del error al usar la regla del trapecio en \([0,0.2]\) es \[ |E_T| \; \approx\; 2.35\times 10^{-3}. \]

Área bajo \(y = \frac{\sin(x)-0.2}{1+x}\) en un intervalo

Problema:

Calcular el área limitada por la gráfica de la función \[ y = \frac{\sin(x)-0.2}{1+x}, \] y la recta \(y=0\) para los valores \(x\) entre \(x_0\) y 1 siendo \(x_0\) el punto de corte de la curva y el eje x en el intervalo \([0,1]\). Para calcular \(x_0\) se debe utilizar el método de Newton con una tolerancia \(10^{-10}\) y para calcular la integral se debe utilizar Simpson compuesta con una tolerancia de \(10^{-8}\).

Solución:

1. Determinación el punto de corte con \(y=0\)

Queremos resolver \[ \frac{\sin(x)-0.2}{1+x} \;=\; 0 \quad\Longleftrightarrow\quad \sin(x) = 0.2. \] En el intervalo \([0,1]\), solo hay una solución \(x_0 \). que se va a calcular numéricamente con el método de Newton con una tolerancia \(10^{-10}\).

Método de Newton

Definimos \[ g(x) \;=\; \sin(x)-0.2, \quad g'(x) = \cos(x). \] A partir de un valor inicial \(x_0 = 0.5\), iteramos \[ x_{k+1} = x_k \;-\; \frac{g(x_k)}{g'(x_k)}, \] hasta \(\lvert x_{k+1} - x_k\rvert < 10^{-10}\).


% ------------------------------------------------------------
% Newton para encontrar la raíz de sin(x)=0.2 en [0,1]
% ------------------------------------------------------------
g   = @(x) sin(x)-0.2;
dg  = @(x) cos(x);

x = 0.5;                % Valor inicial
tol = 1e-10;             % Tolerancia
maxIter = 50;

for k = 1:maxIter
    gx  = g(x);
    dgx = dg(x);
    xk = x - gx/dgx;
    if abs(xk - x) < tol
        break;
    end
    x = xk;
end

raiz = xk;
fprintf('Raíz aproximada (Newton): %.12e (iteraciones: %d)\n', raiz, k);

Después de ejecutar el código, obtenemos:


Raíz aproximada (Newton):  2.013579207903e-01 (iteraciones: 5)

Es decir, la raíz está en \(x_0 \approx 0.201\).

2. Cálculo del área

El área buscada es \[ A = \int_{x_0}^{1} \frac{\sin(x)}{1+x}\,dx, \] donde \(x_0= 2.013579207903e-01\). Debemos aproximar esta integral con error absoluto inferior a \(10^{-8}\) usando la regla compuesta de Simpson con un número de subintervalos \(N\) suficientemente grande. Buscamos este valor,

Estimación del número de subintervalos

La cota de error de Simpson compuesto en \([0,1]\) es \[ \lvert E_S \rvert \;\le\; \frac{N}{90}\,{\left(\frac{b-a}{{2N}}\right)}^5 \max_{x\in[0,1]} \bigl\lvert f^{(4)}(x)\bigr\rvert, \] con \(f(x) = \frac{\sin(x)-0.2}{1+x}\). Calculemos \(f^{(4)}(x)\) de manera simbólica y hallamos su máximo en \([0,1]\).


% ------------------------------------------------------------
% Cálculo simbólico de f4(x) y estimación de su cota en [0,1]
% ------------------------------------------------------------
syms x
f= (sin(x)-0.2) / (1 + x);
f4(x) = diff(f, x, 4);

% Viendo la gráfica de f4, el valor máximo en valor absoluto
% se obtiene en el punto 0
maxAprox = double(max(abs(f4(0))))

% Otra opción: convertir a función numérica
f4 = matlabFunction(f4, 'Vars', x);
% Encontrar máximo de |f4(x)| en [0,1] mediante búsqueda fina
u = linspace(0, 1, 10001);
maxAprox1 = max(abs(f4(u)));

fprintf('Cota |f^{(4)}(x)| ≤ %.6f en [0,1]\n', maxAprox);

Al ejecutar el bloque anterior, se obtiene:


Cota |f^{(4)}(x)| ≤ 24.8 en [0,1]

Por consiguiente, \[ \lvert E_S \rvert \;\le\; \frac{(1-x_0)^5}{90 \cdot 2^5}\,\frac{1}{N^4}\, 24.8 \approx \frac{0.03581}{N^4}. \] Queremos \(0.03581 / N^4 \le 10^{-8}\), de donde \[ N^4 \;\ge\; \frac{0.03582}{10^{-8}} \quad\Longrightarrow\quad N \;\ge\; \bigl(0.03581 \times 10^{8}\bigr)^{1/4} \approx\; 43.50118 \] Tomamos \(N = 44\)

3. Implementación en MATLAB de Simpson compuesto con \(N = 44\)


% ------------------------------------------------------------
% Simpson compuesto con N= 44
% ------------------------------------------------------------
f = @(x) (sin(x)-0.2) ./ (1 + x);

N = 44;                  % número de divisiones del intervalo
h = (1-x0) / N;          % Longitud intervalo donde se aplica Simpson

% Generar nodos x_i y evaluar f
x0=2.013579207903e-01;
x = linspace(x0,1,2*N+1);          % Nodos 2*N+1 : x:0, h/2, h, ..., 1
fx = f(x);

% Sumar impares y pares según Simpson
sum1  = sum(fx(2:2:end-1));  % f(x_1), f(x_3), ..., f(x_{n-1})
sum2 = sum(fx(3:2:end-2));  % f(x_2), f(x_4), ..., f(x_{n-2})

IS = (h/6) * (fx(1) + 4*sum1 + 2*sum2 + fx(end));

fprintf('Simpson compuesto (N=44): %.12f\n', IS);

Al ejecutar este bloque, obtenemos:


Simpson compuesto (N=44): 0.164441612485

Por comparación, por ejemplo con el comando integral en Matlab:


I = integral(f, x0, 1, 'AbsTol', 1e-8);
fprintf('Valor “exacto” numérico: %.12f\n', I);
fprintf('Error absoluto ≈ %.2e\n', abs(I - IS));

Resultado de esa ejecución:


Valor “exacto” numérico: 0.164441612578
Error absoluto ≈ 9.29e-11

El error absoluto es menor que \(10^{-8}\), cumpliendo el requisito.

Resultado final:

  • El punto de corte con el eje \(x\) hallado por Newton es \(x_0 = 2.013579207903e-01\).
  • El área aproximada es \[ A = \int_{x_0}^{1} \frac{\sin(x)}{1+x}\,dx \approx 0.164441612485293, \] con error absoluto menor que \(10^{-8}\).

Regla del trapecio compuesta para \(\log(\sin(x))\) en \([1,\tfrac{\pi}{2}]\)

Problema:

Calcular el valor aproximado de la integral \[ \int_{1}^{\pi/2} \log\bigl(\sin(x)\bigr)\,dx, \] utilizando:

  1. La regla del trapecio compuesta con 10 intervalos.
  2. Determinar cuántos subintervalos serían necesarios para que el error absoluto sea inferior a \(10^{-10}\).

Solución:

1. Parte (a): Trapecio compuesto con \(N=10\) subintervalos

Definimos \[ f(x) \;=\; \log\bigl(\sin(x)\bigr), \quad a = 1,\quad b = \frac{\pi}{2}. \] Con \(N = 10\), el tamaño de cada subintervalo, y la diferencia entre los nodos, es \[ h \;=\; \frac{b - a}{N} \;=\; \frac{\tfrac{\pi}{2} - 1}{10}. \] Los nodos son \[ x_i = a + i\,h, \quad i = 0,1,\dots,10, \] y la regla del trapecio compuesta se escribe: \[ \int_{1}^{\pi/2} \log\bigl(\sin(x)\bigr)\,dx \;\approx\; \frac{h}{2}\Bigl[ f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n) \Bigr]. \] En MATLAB, esto se implementa así:


% ------------------------------------------------------------
% Trapecio compuesto con N=10
% ------------------------------------------------------------
f = @(x) log(sin(x));

N = 10;                               % número de subintervalos
a = 1;
b = pi/2;
h = (b - a) / N;

% Generar nodos x_i y evaluar f
x = a + (0:N)' * h;          % vector columna de 11 puntos: 1, 1+h, ..., π/2
fx = f(x);

% Trapecio compuesto
IT = (h/2) * (fx(1) + 2*sum(fx(2:end-1)) + fx(end));

fprintf('Trapecio compuesto (N=10): %.12f\n', IT);

Al ejecutar este bloque en MATLAB, se obtiene un resultado aproximado:


Trapecio compuesto (N=10): -0.032247145113

Para comparar, un cálculo muy preciso de la integral (por ejemplo usando integral con tolerancia estricta) da:


I = integral(@(x) log(sin(x)), 1, pi/2, 'AbsTol', 1e-12);
fprintf('Valor numérico: %.12f\n', I);
fprintf('Error absoluto ≈ %.2e\n', abs(I - IT));

Este cálculo muestra que el valor “exacto” numérico es aproximadamente:


Valor numérico: -0.032072839160
Error absoluto ≈ 1.74e-04

Por lo tanto, con \(N=10\) el trapecio compuesto da \(\displaystyle -0.032247145113\) y el error absoluto es aproximadamente \(1.74\times10^{-4}\).

2. Parte (b): Número de subintervalos para \(E \le 10^{-10}\)

Considerando \(h=\frac{b-a}{N}\) la cota del error para la regla compuesta del trapecio en \([a,b]\): \[ \bigl\lvert E_T \bigr\rvert \;\le\; \frac{N \;h^3}{12}\,\max_{x\in[a,b]} \bigl\lvert |f''(x)| \bigr\rvert, \] donde \(f(x) = \log(\sin(x))\). Calculamos la segunda derivada:

\[ f'(x) = \frac{d}{dx}\bigl[\log(\sin(x))\bigr] = \cot(x), \] \[ f''(x) = \frac{d}{dx}\bigl[\cot(x)\bigr] = -\,\csc^2(x) = -\,\frac{1}{\sin^2(x)}. \]

En el intervalo \([1,\tfrac{\pi}{2}]\), \(\sin(x)\) es creciente desde \(\sin(1)\approx0.84147\) hasta \(\sin(\pi/2)=1\). Por tanto, \[ \bigl\lvert f''(x)\bigr\rvert = \frac{1}{\sin^2(x)} \;\le\; \frac{1}{\sin^2(1)} \;=\; \frac{1}{(0.84147)^2} \approx 1.412. \] Sea \(M = \max_{x\in[1,\pi/2]}|f''(x)| \approx 1.412.\) Además, \(b - a = \frac{\pi}{2} - 1.\)

La cota del error es: \[ \frac{N}{12} \frac{(b-a)^3}{N^3}\,M = \frac{\bigl(\tfrac{\pi}{2} - 1\bigr)^3}{12\,N^2}\,M. \] Como queremos \[ \frac{\bigl(\tfrac{\pi}{2} - 1\bigr)^3}{12\,N^2}\,M \;\le\; 10^{-10}. \] De aquí, \[ N^2 \;\ge\; \frac{\bigl(\tfrac{\pi}{2} - 1\bigr)^3\,M}{12\,\times\,10^{-10}}. \] Calculando la raiz cuadrada \[ N \approx 14795. \] Por tanto, se requieren al menos \(N = 14795\) subintervalos (redondeando al entero siguiente) para garantizar que el error absoluto sea inferior a \(10^{-10}\).


% ------------------------------------------------------------
% Cálculo de n mínimo para |E| <= 1e-10 en trapecio compuesto
% ------------------------------------------------------------
M = 1 / (sin(1)^2);                  % cota en [1, π/2], ~1.412
a = 1;
b = pi/2;
num = sqrt(((b - a)^3 * M) *10^10/ 12);
N = ceil(num);

fprintf('N = %d', N);

Con \(N = 14795\), el paso es \(h \approx \frac{0.5707963}{14795} \approx 3.858\times10^{-5}\).

3. Verificación rápida con MATLAB (opcional)

Para comprobar que con \(N=14795\) el trapecio compuesto alcanza la precisión deseada, podríamos ejecutar:


% ------------------------------------------------------------
% Verificación numérica de la integral con N = 14795
% ------------------------------------------------------------
f = @(x) log(sin(x));
a = 1; b = pi/2; N=14795;
h = (b - a)/N;
x = linspace(a,b,N+1);
fx = f(x);

IT = (h/2) * (fx(1) + 2*sum(fx(2:end-1)) + fx(end));
I  = integral(@(x) log(sin(x)), a, b, 'AbsTol', 1e-12);
error  = abs(I - IT);

fprintf('Trapecio compuesto (n=14795): %.12f\n', IT);
fprintf('Valor de referencia numerico:   %.12f\n', I);
fprintf('Error absoluto:                %.2e\n', error);

Este código confirma que con \(N=14795\) el error es efectivamente menor que \(10^{-10}\).

Resultado final:

  • Usando \(N = 10\), la regla del trapecio compuesta da \(\displaystyle \int_{1}^{\pi/2} \log(\sin(x))\,dx \approx -0.032072839240\), con un error aproximado de \(1.74\times10^{-4}\).
  • Para garantizar \(\lvert E\rvert \le 10^{-10}\) se requieren al menos \(n = 14795\) subintervalos en la regla del trapecio compuesta.

Trapecio Adaptativa para varias integrales

Problema:

Calcular las siguientes integrales con la función TrapecioAdaptativa e ir comprobando el valor con el comando integral de MATLAB:

  1. \(\displaystyle \int_{0}^{\pi} e^{\cos(x)}\,\sin(x)\,dx,\quad E_a \le 10^{-12}.\)
  2. \(\displaystyle \int_{0}^{4} e^{2x}\,dx,\quad E_a \le 10^{-12}.\)
  3. \(\displaystyle \int_{0}^{\pi} \bigl(1 + x^2\bigr)\,e^{\sin(x)}\,dx,\quad E_a \le 10^{-12}.\)

Solución:

Se utilizaría la función TrapecioAdaptativa siguiente


function valor=TrapecioAdaptativa(f,a,b,tol)
	%Aproxima la integral de f en el intervalo [a,b]
	%usando la regla del trapecio adaptativa
	%en forma recursiva
	xm=(a+b)/2;
	%Trapecios en [a,b]
	I1=(f(a)+f(b))*(b-a)/2;
	%Trapecios en [a,xm] y en [xm,b]
	I2=(f(a)+f(xm))*(xm-a)/2+(f(xm)+f(b))*(xm-a)/2;
	%Estimación del error
	error=abs(I1-I2)/3;
	if error < tol
		valor=I2;
	else
		valor= TrapecioAdaptativa(f,a,xm,tol/2)+
		TrapecioAdaptativa(f,xm,b,tol/2);
	end
end

Apartado a). Cálculo con TrapecioAdaptativa: \(\int_{0}^{\pi} e^{\cos(x)}\,\sin(x)\,dx\)

f = @(x) exp(cos(x)).*sin(x);
a = 0;
b = pi;
Ea = 1e-12;
Ia = TrapecioAdaptativa(f, a, b);

Se obtiene:

Ia ≈ 2.350402387287603

Verificación con integral:

Ia = integral(f, a, b);

El comando integral devuelve:

Ia ≈ 2.3504023872876029

Comparación: Ambos valores coinciden dentro de la tolerancia solicitada. Además, analíticamente:

\[ \int_{0}^{\pi} e^{\cos(x)}\,\sin(x)\,dx \;=\; \Bigl[-\,e^{\cos(x)}\Bigr]_{0}^{\pi} \;=\; \bigl(-e^{-1}\bigr)\;-\;\bigl(-e^{1}\bigr) \;=\; e \;-\; e^{-1} \;\approx\; 2.3504023872876029. \]

Apartado b). Cálculo de \(\displaystyle \int_{0}^{4} e^{2x}\,dx\) con TrapecioAdaptativa (considerando Ea<):

f = @(x) exp(2*x);
a = 0;
b = 4;
Ea = 1e-12;
Er = 1e-10;
Ib= TrapecioAdaptativa(f, a, b, Ea, Er);

Se obtiene:

Ib ≈ 1489.978993520864

Verificación con integral:

Ib = integral(f, a, b);

El comando integral devuelve:

Ib ≈ 1489.9789935208643

Comparación: Ambos valores coinciden dentro de las tolerancias solicitadas. Analíticamente:

\[ \int_{0}^{4} e^{2x}\,dx \;=\; \frac{e^{2x}}{2}\Bigg|_{0}^{4} \;=\; \frac{e^{8} - 1}{2} \;\approx\; 1489.9789935208643. \]

Apartado c) Cálculo de \(\displaystyle \int_{0}^{\pi} \bigl(1 + x^2\bigr)\,e^{\sin(x)}\,dx\) con TrapecioAdaptativa:

f = @(x) (1 + x.^2).*exp(sin(x));
a = 0;
b = pi;
Ea = 1e-12;
Ic= TrapecioAdaptativa(f, a, b, Ea);

Se obtiene:

Ic ≈ 25.366042306612423

Verificación con integral:

Ic = integral(f, a, b);

El comando integral devuelve:

Ic ≈ 25.366042306612424

Comparación: El valor coincide dentro de la tolerancia solicitada.

Integrales impropias en \([0,\infty)\)

Problema:

Calcular el valor aproximado de las siguientes integrales:

  1. \(\displaystyle \int_{0}^{\infty} e^{-x^2}\,\cos(x)\,dx,\quad E_a \le 10^{-12}.\)
  2. \(\displaystyle \int_{0}^{\infty} \log(2+x)\,e^{-x}\,dx,\quad E_a \le 10^{-12}.\)
  3. \(\displaystyle \int_{0}^{\infty} \frac{\sin(x)}{1+x^3}\,dx,\quad E_a \le 10^{-8}.\)
  4. \(\displaystyle \int_{0}^{\infty} \frac{e^{-x^4}}{1+x}\,dx,\quad E_a \le 10^{-12}.\)

Solución:

Para aproximar integrales en \([a,\infty)\) con la precisión exigida, calcularemos en primer lugar \(b\) con \(b \gt a\) para que el valor absoluto de la integral de la función en \([b, \infty]\) sea inferior a la tolerancia dividido entre 2. Con ese valor, se calculará con una fórmula adaptativa el valor de la integral entre \([a, b]\) con error inferior a \(tol/2\). Es decir, dado que

\[\int\limits_a^\infty {f\left( x \right)dx} = \underbrace {\int\limits_a^b {f\left( x \right)dx} }_{{I_1}} + \underbrace {\int\limits_b^\infty {f\left( x \right)dx} }_{{I_2}}\] se tendrá que \[\left| | \int\limits_{a}^{\infty }{f\left( x \right)dx} \right| |\le \underbrace{\left| | \int\limits_{a}^{b}{f\left( x \right)dx} \right| |}_{{{I}_{1}}}+\underbrace{\left| | \int\limits_{b}^{\infty }{f\left( x \right)dx} \right| |}_{{{I}_{2}}}\le \frac{\varepsilon }{2}+\frac{\varepsilon }{2}=\varepsilon \]

Comprobaremos el resultado luego con el comando integral de MATLAB con la opción 'AbsTol' adecuada.

Se puede considerar la siguiente función que implementa el método de Simpson Adaptativo de forma recursiva. Puede aplicarse también el método adaptativo general con otra fórmula de cuadratura.


function I=SimpsonAdaptativa(f,a,b,tol)
 %Simpson en [a,b]
 I1=Simpson(f,a,b);
 m=(a+b)/2;
 %Simpson en [a,puntoMedio] y [puntoMedio,b]
 I2=Simpson(f,a,m)+Simpson(f,m,b);
 %Estimación del error
 E=abs(I2-I1)/15;
 if E < tol
     I=I2;
 else
     I1=SimpsonAdaptativa(f,a,m,tol/2);
     I2=SimpsonAdaptativa(f,m,b,tol/2);
     I=I1+I2;
 end
end

(a) \(\displaystyle I_a = \int_{0}^{\infty} e^{-\,x^2}\,\cos(x)\,dx,\quad E_a \le 10^{-12}.\)

Implementación numérica en MATLAB (división en dos tramos):

Dado que

$$|{I_2}|\, = \left| | {\int\limits_b^\infty {{e^{ - {x^2}}}\cos \left( x \right)dx\,} } \right| | \le \int\limits_b^\infty {\left| | {{e^{ - {x^2}}}\cos \left( x \right)} \right| |dx} \le \int\limits_b^\infty {{e^{ - x}}dx} {\mkern 1mu} = {1 \over {{e^b}}}{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} si{\mkern 1mu} {\mkern 1mu} b > 1$$ $$\left| {{I_2}} \right| \le {1 \over {{e^b}}} \le {{{{10}^{ - 12}}} \over 2}$$

Esta desigualdad de cumple para \(b >28.400000000000134\). Tomamos \(b=28.45\).


% ------------------------------------------------------------
% (a) Calculo con error Ea = 1e-12
% ------------------------------------------------------------
f = @(x) exp(-x.^2) .* cos(x);
tol=10^-12;
b=28.45

% Integración adaptativa [0, b]
Ia=SimpsonAdaptativa(f,0,b,tol/2)

% Valor "exacto" numérico con Matlab
I1= integral(f, 0, Inf, 'AbsTol', 1e-12);

fprintf('Adaptativo en [0,b]: I ≈ %.15f\n', Ia);
fprintf('    Valor numérico            = %.15f\n', I1);
fprintf('    Error absoluto          = %.2e\n', abs(Ia - I1));

Adicional: Conocimiento analítico: Existe la fórmula conocida: \[ \int_{0}^{\infty} e^{-\,x^2}\,\cos(b\,x)\,dx = \frac{\sqrt{\pi}}{2}\,e^{-\tfrac{b^2}{4}}, \quad b=1 \;\Longrightarrow\; I_a = \frac{\sqrt{\pi}}{2}\,e^{-1/4}. \] Numéricamente, \[ I_a = \frac{\sqrt{\pi}}{2}\,e^{-1/4} \approx 0.690194223521571\,. \] Así, el valor “exacto” de referencia es \[ I_{\text{exacto}} = \frac{\sqrt{\pi}}{2}\,e^{-1/4} \approx 0.690194223521571. \]

Resultado para (a):

  • \(\displaystyle Ia \approx 0.690194223521577\)

(b) \(\displaystyle I_b = \int_{0}^{\infty} \log(2 + x)\,e^{-x}\,dx,\quad E_a \le 10^{-12}.\)

Observaciones:

Dado que

$$\left| | {\log \left( {2 + x} \right){e^{ - x}}} \right| | \le x\,{e^{ - x}} \,\,\,\,si\,\,\,\,x \gt x_0 \approx 1.14 $$ $$|{I_2}|\, = \left| | {\int\limits_b^\infty {\log \left( {2 + x} \right){e^{ - x}}dx\,} } \right| | \le \int\limits_b^\infty {x{e^{ - x}}dx} = {{b + 1} \over {{e^b}}}$$ $$\left| {{I_2}} \right| \le {{b + 1} \over {{e^b}}} \le {{{{10}^{ - 12}}} \over 2}$$

Esta desigualdad de cumple para \(b >31.900000000000183\). Tomamos \(b=31.95\).

Implementación en MATLAB:


% ------------------------------------------------------------
% (b) Calculo con error Ea = 1e-12
% ------------------------------------------------------------
f = @(x) log(2 + x) .* exp(-x);
tol=10^-12;
b=31.95

% Integración adaptativa [0, b]
Ib=SimpsonAdaptativa(f,0,b,tol/2)

% Valor "exacto" numérico
I1= integral(f, 0, Inf, 'AbsTol', 1e-12);

fprintf('Adaptativo en [0,b]: I ≈ %.15f\n', Ib);
fprintf('    Error absoluto          = %.2e\n', abs(Ib - I1));

Resultado para (b):

  • \(\displaystyle Ib \approx 1.054475797448241\,.\)

(c) \(\displaystyle I_d = \int_{0}^{\infty} \frac{\sin(x)}{1 + x^3}\,dx,\quad E_a \le 10^{-8}.\)

Cálculo:

Dado que

$$|{I_2}|\, = \left| | {\int\limits_b^\infty {{{\sin \left( x \right)} \over {1 + {x^3}}}dx\,} } \right| | \le \int\limits_b^\infty {\left| | {{{\sin \left( x \right)} \over {1 + {x^3}}}} \right| |dx\,} \le \int\limits_b^\infty {{1 \over {{x^3}}}dx = {1 \over {2{b^2}}}\,} $$ $$\left| {{I_2}} \right| = {1 \over {2{b^2}}}\, \le {{{{10}^{ - 8}}} \over 2}$$

Esta desigualdad de cumple para \(b > 10000\). Tomamos \(b= 10000\).

Implementación en MATLAB:


% ------------------------------------------------------------
% (c) Calculo con error Ea = 1e-10
% ------------------------------------------------------------
f = @(x) sin(x) ./ (1 + x.^3);
tol=10^-8;
b= 10000;

% Integración adaptativa [0, b]
Ic=SimpsonAdaptativa(f,0,b,tol/2)

% Valor "exacto" numérico
I1= integral(f, 0, Inf, 'AbsTol', 1e-12);

fprintf('Adaptativo en [0,b]: I ≈ %.15f\n', Id);
fprintf('    Error absoluto          = %.2e\n', abs(Ic - I1));

Resultado para (c):

  • \(\displaystyle I_c \approx 0.199936086\) (con 8 decimales de precisión).

(d) \(\displaystyle I_e = \int_{0}^{\infty} \frac{e^{-\,x^4}}{1 + x}\,dx,\quad E_a \le 10^{-12}.\)

Cálculo:

Dado que

$$|{I_2}|\, = \left| | {\int\limits_b^\infty {{{{e^{ - {x^4}}}} \over {1 + x}}dx\,} } \right| | \le \int\limits_b^\infty {{{{e^{ - {x^4}}}} \over x}dx\,} = \int\limits_b^\infty {{{{x^3}{e^{ - {x^4}}}} \over {{x^4}}}dx \le \,} \int\limits_b^\infty {{{{x^3}{e^{ - {x^4}}}} \over {{b^4}}}dx = {{4{e^{ - {b^4}}}} \over {{b^4}}}\,} $$ $$\left| | {{I_2}} \right| | = {{4{e^{ - {b^4}}}} \over {{b^4}}}\, \le {{{{10}^{ - 12}}} \over 2}$$

Esta desigualdad de cumple para \(b >2.3\). Tomamos \(b=2.3\).

Implementación en MATLAB:


% ------------------------------------------------------------
% (d) Calculo con error Ea = 1e-12
% ------------------------------------------------------------
f = @(x) exp(-x.^4) ./ (1 + x);
b = 2.3;
tol=10^-12;

% Integración adaptativa [0, b] con tolerancia tol/2
Id=SimpsonAdaptativa(f,0,b,tol/2)

% Valor "exacto" numérico
I1= integral(f, 0, Inf, 'AbsTol', 1e-12);

fprintf('Adaptativo en [0,b]: I ≈ %.15f\n', Id);
fprintf('    Error absoluto          = %.2e\n', abs(Id - I1));

Resultado para (d):

  • \(\displaystyle Id \approx 0.760177910377389\).

Resumen de resultados:

  • (a) \(\displaystyle \int_{0}^{\infty} e^{-x^2}\cos(x)\,dx = \frac{\sqrt{\pi}}{2}\,e^{-1/4} \approx 0.6894131219478915,\) obtenido con error \(<10^{-12}\).
  • (b) \(\displaystyle \int_{0}^{\infty} \log(2+x)\,e^{-x}\,dx \approx 1.099984214827246,\) con error \(<10^{-12}\).
  • (c) \(\displaystyle \int_{0}^{\infty} x^4\,e^{-x}\,dx = 24,\) exacto, verificado numéricamente con error \(<10^{-12}\).
  • (d) \(\displaystyle \int_{0}^{\infty} \frac{\sin(x)}{1 + x^3}\,dx \approx 0.199936086,\) con error \(<10^{-8}\).
  • (e) \(\displaystyle \int_{0}^{\infty} \frac{e^{-x^4}}{1 + x}\,dx \approx 0.760177910377389,\) con error \(<10^{-12}\).