6.6: Métodos de integración numérica

Última actualización
Guardar como PDF

Page ID: 110339

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

La sección 5.2 mostró cómo obtener valores exactos para integrales definidas de algunas funciones simples (polinomios de bajo grado) mediante el uso de áreas de rectángulos. Para las funciones sin antiderivada de forma cerrada, el método del rectángulo normalmente produce un valor aproximado de la integral definida: cuanto más rectángulos, mejor será la aproximación.

Por ejemplo, supongamos que quería evaluar

\[\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx\]con el método del rectángulo. Esto significa encontrar el área de la región sombreada en la figura de la derecha. Las computadoras han eliminado la necesidad de hacer este tipo de cálculos a mano. Aunque el método rectangular es simple de implementar en un lenguaje de programación tradicional (por ejemplo, a través de una construcción de bucle), hay formas más fáciles en un lenguaje específico de dominio (DSL) orientado a la computación científica. Uno de esos DSL es MATLAB, o su clon libre de código abierto Octave. ¹⁰

Implementar el método del rectángulo desde cero en Octave es un one-liner. Por ejemplo, supongamos que divide el intervalo\(\ival{0}{\sqrt{\pi}}\) en 100,000 (\(10^5\)) subintervalos de igual longitud, produciendo 100,001 (1 + 1e5 en notación científica) puntos igualmente espaciados en\(\ival{0}{\sqrt{\pi}}\) (incluyendo\(0\) y\(\sqrt{\pi}\,\)). Primero use los extremos izquierdos de los\(10^5\) subintervalos:

octave> sum(sin(linspace(0,sqrt(pi),1+1e5)(1:end-1).^2)*sqrt(pi)/1e5)
ans = 0.8948314693913354

Ahora usa los puntos finales correctos:

octave> sum(sin(linspace(0,sqrt(pi),1+1e5)(2:end).^2)*sqrt(pi)/1e5)
ans = 0.8948314693913354

Finalmente, use los puntos medios de los subintervalos:

octave> sum(sin((linspace(0,sqrt(pi),1+1e5)(1:end-1)+sqrt(pi)/2e5).^2)*
        sqrt(pi)/1e5)
ans = 0.8948314695305527

El verdadero valor de la integral hasta 15 decimales es 0.894831469484145, por lo que las tres aproximaciones son exactas a 9 decimales. La sintaxis en los comandos anteriores se puede explicar con algunos ejemplos. El siguiente comando crea 4 puntos igualmente espaciados en el intervalo\(\ival{1}{7}\) (incluyendo\(x=1\) y\(x=7\)),\(\ival{1}{7}\) dividiendo así en 3 subintervalos cada uno de longitud\((7-1)/3 = 2\):

octave> linspace(1,7,4)
ans =

  1  3  5  7

Obtener todos menos el último número en la lista anterior: ¹¹

octave> linspace(1,7,4)(1:end-1)
ans =

  1  3  5

Ahora cuadre cada número en esa lista de números (el punto antes del operador de exponenciación aplica la operación de cuadratura 2 elementos en la lista):

octave> linspace(1,7,4)(1:end-1).^2
ans =

  1  9  25

Ahora toma el seno de cada uno de esos números cuadrados (medidos en radianes):

octave> sin(linspace(1,7,4)(1:end-1).^2)
ans =

  0.8414709848078965  0.4121184852417566  -0.132351750097773

Ahora multiplica cada uno de esos números (las alturas de los rectángulos) por el ancho (2) de cada rectángulo, luego suma esas áreas:

octave> sum(sin(linspace(1,7,4)(1:end-1).^2)*2)
ans = 2.24247543990376

Antes del advenimiento de la computación moderna, el método del rectángulo se consideraba ineficiente, por lo que se crearon métodos alternativos. Dos de estos métodos son la regla trapezoidal y la regla de Simpson. La idea detrás de ambos métodos es aprovechar la pendiente cambiante de una función no lineal mediante el uso de regiones no rectangulares. Para la regla trapezoidal esas regiones son trapezoides, mientras que la regla de Simpson usa regiones cuasirectangulares cuyos bordes superiores son parábolas, como se muestra en la Figura [fig:nummethods]:

Para una partición\(P= \lbrace a=x_0 < x_1 < \cdots < x_{n-1} < x_n = b \rbrace\) de un intervalo\(\ival{a}{b}\) en\(n \ge 1\) subintervalos de igual ancho\(h = (b-a)/n\), let\(y_i = f(x_i)\) for\(i = 0\),\(1\),\(\ldots\),\(n\). La regla trapezoidal suma las áreas de trapecios en cada subintervalo\(\ival{x_i}{x_{i+1}}\), siendo el borde superior el segmento de línea que une los puntos\((x_i,y_i)\) y\((x_{i+1},y_{i+1})\). La fórmula de aproximación es sencilla de derivar, basada en áreas de trapecios:

La regla de Simpson depende de pares de subintervalos vecinos:\(\ival{x_0}{x_1}\) y\(\ival{x_1}{x_2}\),\(\ival{x_2}{x_3}\) y\(\ival{x_3}{x_4}\),\(\ldots\),\(\ival{x_{n-2}}{x_{n-1}}\) y\(\ival{x_{n-1}}{x_n}\). Así,\(n \ge 2\) debe ser parejo. El borde superior de la región sobre cada par\(\ival{x_{i}}{x_{i+1}}\) y\(\ival{x_{i+1}}{x_{i+2}}\) es la parábola única que une los 3 puntos\((x_i,y_i)\),\((x_{i+1},y_{i+1})\), y\((x_{i+2},y_{i+2})\). La fórmula de aproximación es entonces: ¹²

Ejemplo\(\PageIndex{1}\): trapsimp

Agrega texto aquí.

Solución

Aproximar el valor de\(~\displaystyle\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx~\) usando la regla trapezoidal y la regla de Simpson con\(n=10^5\) subintervalos.

Solución: Desde\(x_0 = 0\) y\(x_n = \sqrt{\pi}\), entonces\(y_0 = \sin\,(x_0^2) = \sin\,0 = 0\) y\(y_n = \sin\,(x_n^2) = \sin\,\pi = 0\). Así,\(y_0\) y no\(y_n\) aportan nada a las fórmulas sumatorias para ambas reglas. En particular, la aproximación de la regla trapezoidal se convierte en

\[\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx ~\approx~ \frac{h}{2}\,(0 ~+~ 2y_1 ~+~ 2y_2 ~+~ \cdots ~+~ 2y_{n-1} ~+~ 0) ~=~ h\,\cdot\,\left(\sum_{k=1}^{n-1} y_k\right)\]que es simple de implementar en Octave:

octave> x = linspace(0,sqrt(pi),1+1e5);
octave> h = sqrt(pi)/1e5;
octave> h*sum(sin(x(2:end-1).^2))
ans = 0.8948314693913405

Asimismo, la aproximación de la regla de Simpson se convierte en

\[\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx ~\approx~ \frac{h}{3}\,(4y_1 ~+~ 2y_2 ~+~ 4y_3 ~+~ 2y_4 ~+~ \cdots ~+~ ~+~ 2y_{n-2} ~+~ 4y_{n-1})\]que se puede implementar fácilmente mediante el uso de las potentes funciones de indexación de Octave:

octave> (h/3)*(4*sum(sin(x(2:2:end-1).^2)) + 2*sum(sin(x(3:2:end-1).^2)))
ans = 0.8948314694841457

En el comando anterior, la instrucción x (3:2:end-1) le permite omitir todos los demás elementos de la lista x después de la posición 3, moviendo hacia arriba la lista en incrementos de 2 posiciones hasta la posición siguiente a la última posición de la lista (end-1). De manera similar para x (2:2:fin-1), que comienza en la posición 2 y luego se mueve hacia arriba en incrementos de 2.

Tenga en cuenta que la regla de Simpson da esencialmente el valor verdadero en este caso, y el valor de la regla trapezoidal es prácticamente el mismo que el valor producido por la función trapz incorporada en Octave/matlab:

octave> trapz(x,sin(x.^2))
ans = 0.8948314693913402

En general, es mejor usar este tipo de funciones integradas en lugar de implementar las suyas propias.

Por lo general, la regla de Simpson es ligeramente más eficiente que la regla trapezoidal, que es un poco más eficiente que el método del rectángulo. Sin embargo, en los ejemplos anteriores todas las aproximaciones fueron exactas a al menos 9 decimales (equivalente a conseguir que la distancia entre Detroit y Chicago sea correcta dentro del grosor de un palillo de dientes). El tiempo de ejecución de cada cálculo fue de sólo unas pocas milésimas de segundo. La computación moderna generalmente ha hecho que las diferencias de eficiencia sean insignificantes. Observe que las aproximaciones en el método rectángulo, la regla trapezoidal y la regla de Simpson pueden escribirse como combinaciones lineales de valores de función\(f(a_i)\) multiplicados por “pesos”\(w_i\):

\[\int_a^b f(x)~\dx ~\approx~ \sum_{i=0}^{n} w_i\,f(a_i)\]Por ejemplo, los pesos en la regla de Simpson son\(w_i = \frac{h}{3}\)\(\frac{2h}{3}\), o\(\frac{4h}{3}\), dependiendo de los puntos\(a_i\) en el intervalo\(\ival{a}{b}\). El método de cuadratura gaussiana transforma una integral sobre cualquier intervalo\(\ival{a}{b}\) en una integral sobre el intervalo específico\(\ival{-1}{1}\) y luego usa un conjunto estándar de puntos\(\ival{-1}{1}\) y pesos conocidos para esos puntos: ¹³

[tbl:gaussquad]

Ejemplo\(\PageIndex{1}\): gaussquad

Agrega texto aquí.

Solución

Aproximar el valor de\(~\displaystyle\int_0^2 \dfrac{\dx}{1 + x^3}~\) usando cuadratura gaussiana con\(n=4\) puntos.

Solución: Para\(a=0\) y\(b=2\), utilizar la sustitución\(u = \frac{1}{b-a}(2x - a - b) = x-1\), para que\(x=u+1\) y\(\dx = \du\). Así,\(g(u) = f(u+1) = \frac{1}{1 + (u+1)^3}\). Usando\(n=4\) en la Tabla [tbl:gaussquad], los puntos\(a_i\) y pesos\(w_i\) son

\[\begin{aligned} {3} a_1 ~&=~ -0.339981 \quad\quad & w_1 ~&=~ 0.652145\\ a_2 ~&=~ 0.339981 \quad\quad & w_2 ~&=~ 0.652145\\ a_3 ~&=~ -0.861136 \quad\quad & w_3 ~&=~ 0.347855\\ a_4 ~&=~ 0.861136 \quad\quad & w_4 ~&=~ 0.347855\end{aligned}\]y así

\[\begin{aligned} \int_0^2 \frac{\dx}{1 + x^3} ~&=~ \frac{2-0}{2}\,\int_{-1}^1 g(u)~\du ~=~ \int_{-1}^1 \frac{\du}{1 + (u+1)^3} ~\approx~ \sum_{i=1}^4 w_i\,g(a_i) ~=~ \sum_{i=1}^4 w_i\,\cdot\,\frac{1}{1 + (a_i+1)^3}\

\ [6pt] &\ approx~\ frac {0.652145} {1 + (-0.339981+1) ^3} ~+~\ frac {0.652145} {1 + (0.339981+1) ^3} ~+~\ frac {0.347855} {1 + (-0.861136+1) ^3} ~+~\ frac {0.347855} {1 + (0.861136+1) ^3}\

\ [4pt] &\ approx~ 1.091621\ end {aligned}\] El verdadero valor de la integral a seis decimales es 1.090002.
Usar más puntos (por ejemplo\(n=7\)) es fácil de implementar en Octave, usando operaciones elementales en matrices:

octave> a = [ 0 -0.405845 0.405845 -0.741531 0.741531 -0.949108 0.949108 ];
octave> w = [ 0.417959 0.381830 0.381830 0.279705 0.279705 0.129485 0.129485 ];
octave> sum(w./(1 + (a+1).^3))
ans = 1.090016688064804

La cuadratura gaussiana se puede aplicar a integrales inadecuadas. Por ejemplo,

\[\int_0^{\infty} f(x)\,e^{-x}\,\dx ~\approx~ \sum_{i=1}^n w_i\,f(a_i)\]usando los puntos\(a_i\) y pesos\(w_i\) de la Tabla [tbl:gaussquadimp] ¹⁴ para\(n=3, 4\), o 5 puntos en\(\lival{0}{\infty}\):

[tbl:gaussquadimp]

Ejemplo\(\PageIndex{1}\): gqinf

Agrega texto aquí.

Solución

Aproximar el valor de\(~\displaystyle\int_0^{\infty} x^5\,e^{-x}\,\dx~\) usando cuadratura gaussiana con\(n=3\) puntos en la Tabla [tbl:gaussquadimp].

Solución: Para\(n=3\), Tabla [tbl:gaussquadimp] da\(a_1=0.415775\),,\(a_2=2.294280\)\(a_3=6.289945\), y\(w_1=0.711093\),,\(w_2=0.278518\),\(w_3=0.010389\). Entonces para\(f(x)=x^5\),

\[\begin{aligned} \int_0^{\infty} x^5\,e^{-x}\,\dx ~&\approx~ \sum_{i=1}^n w_i\,f(a_i) ~=~ \sum_{i=1}^n w_i\,a_i^5\

\ [4pt] &\ approx~ 0.711093\, (0.415775) ^5 ~+~ 0.278518\, (2.294280) ^5 ~+~ 0.010389\, (6.289945) ^5\\ &\ aprox~ 119.9974709727211\ end {alineado}\] El verdadero valor es\(\Gamma\,(6) = 5! = 120\).
Nota: Los puntos\(a_i\) en la Tabla [tbl:gaussquadimp] son las raíces de los polinomios de Laguerre de grado\(n\).

[sec6dot6]

[exer:pendper] Un péndulo simple de longitud\(l\) oscila a través de un ángulo de a cada\(90\Degrees\) lado de la vertical con punto\(P\), dado por

\[P ~=~ 4\,\sqrt{\dfrac{l}{g}}\,\int_0^{\pi/2} \frac{\dtheta}{\sqrt{1 \;-\; 0.5\,\sin^2 \theta}}\]donde\(g = 9.8\) m/s ² es la aceleración debida a la gravedad. Usa el método rectángulo (con extremos izquierdos), la regla trapezoidal y la regla de Simpson para escribir\(P\) como un múltiplo constante de\(\sqrt{l/g}\). Preferiblemente, use una computadora y\(n= 10^5\) subintervalos de igual ancho (o\(n=10\) subintervalos si se calcula a mano).

Repita Ejercicio [exer:pendper] usando cuadratura gaussiana con\(n=5\) puntos.

[exer:gqsinx2] Aproximar el valor de\(~\displaystyle\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx~\) usando cuadratura gaussiana con\(n=7\) puntos.

Repita Ejercicio [exer:gqsinx2] con\(n=9\) puntos.

Repita Ejercicio [exer:gqsinx2] con\(n=10\) puntos.

Los puntos\(a_i\) en Table [tbl:gaussquad] para cuadratura gaussiana son las raíces de los polinomios de Legendre\(P_n(x)\)\(P_0(x) = 1\), con\(P_1(x) = x\), y\(P_n(x)\) definidos para enteros\(n \ge 2\) por la fórmula de recursión

\[n\,P_{n}(x) ~=~ (2n-1)\,x\,P_{n-1}(x) ~-~ (n-1)\,P_{n-2}\,(x)\]

Escriba\(P_n(x)\) explícitamente en forma polinomio estándar para\(n = 2, 3, 4, 5\).
Verificar que las raíces de\(P_n(x)\) emparejar los\(n\) puntos\(a_1\)\(\ldots\),,\(a_n\) en Tabla [tbl:gaussquad] para\(n = 2, 3, 4, 5\).
Sin cálculos, explique por qué\(~\int_{-1}^1 P_2(x)\,P_3(x)\,\dx = \int_{-1}^1 P_3(x)\,P_4(x)\,\dx = \int_{-1}^1 P_4(x)\,P_5(x)\,\dx = 0\).
Para\(n=0, 1, 2, 3\), verificar que
\[\int_{-1}^1 P_n^2(x)~\dx ~=~ \frac{2}{2n+1} ~.\]

Repetir Ejemplo

Ejemplo\(\PageIndex{1}\): gqinf

Agrega texto aquí.

Solución

con\(n=4\) puntos.

Repetir Ejemplo

Ejemplo\(\PageIndex{1}\): gqinf

Agrega texto aquí.

Solución

con\(n=5\) puntos.

Utilice la Tabla [tbl:gaussquadimp] para aproximar el valor de\(~\displaystyle\int_0^{\infty} \ln\,(1+x)\,e^{-x}\,\dx~\) con\(n=5\) puntos.