CÁLCULO

CÁLCULO INTEGRAL

Método de Romberg genera una matriz triangular cuyos elementos son estimaciones numéricas de la integral definida siguiente:

${\displaystyle \int _{a}^{b}f(x)\,dx}$

usando la extrapolación de Richardson de forma reiterada en la regla del trapecio. El método de Romberg evalúa el integrando en puntos equiespaciados del intervalo de integración estudiado. Para que este método funcione, el integrando debe ser suficientemente derivable en el intervalo, aunque se obtienen resultados bastante buenos incluso para integrandos poco derivables. Aunque es posible evaluar el integrando en puntos no equiespaciados, en ese caso otros métodos como la cuadratura gaussiana o la cuadratura de Clenshaw–Curtis son más adecuados.

Método[editar]

El método se define de forma recursiva así:

${\displaystyle R(0,0)={\frac {1}{2}}(b-a)(f(a)+f(b))}$

${\displaystyle R(n,0)={\frac {1}{2}}R(n-1,0)+h_{n}\sum _{k=1}^{2^{n-1}}f(a+(2k-1)h_{n})}$

${\displaystyle R(n,m)=R(n,m-1)+{\frac {1}{4^{m}-1}}(R(n,m-1)-R(n-1,m-1))}$

o

${\displaystyle R(n,m)={\frac {1}{4^{m}-1}}(4^{m}R(n,m-1)-R(n-1,m-1))}$

donde

${\displaystyle n\geq 1}$

${\displaystyle m\geq 1}$

${\displaystyle h_{n}={\frac {b-a}{2^{n}}}.}$

La cota superior asintótica del error de R(n,m) es:

${\displaystyle O\left(h_{n}^{2^{m+1}}\right).}$

La extrapolación a orden cero ${\displaystyle R(n,0)}$ es equivalente a la Regla del trapecio con ${\displaystyle n+2}$ puntos. a orden uno ${\displaystyle R(n,1)}$ es equivalente a la Regla de Simpson con ${\displaystyle n+2}$ puntos.

Cuando la evaluación del integrando es numéricamente costosa, es preferible reemplazar la interpolación polinómica de Richardson por la interpolación racional propuesta por Bulirsch & Stoer.

Implementación en Python del método de Romberg[editar]

La siguiente es una implementación del método de Romberg en Python:

def print_row(lst):
    print ' '.join('%11.8f' % x for x in lst)

def romberg(f, a, b, eps = 1E-8):
    """Approximate the definite integral of f from a to b by Romberg's method.
    eps is the desired accuracy."""
    R = [[0.5 * (b - a) * (f(a) + f(b))]]  # R[0][0]
    print_row(R[0])
    n = 1
    while True:
        h = float(b-a)/2**n
        R.append((n+1)*[None])  # Add an empty row.
        R[n][0] = 0.5*R[n-1][0] + h*sum(f(a+(2*k-1)*h) for k in range(1, 2**(n-1)+1)) # for proper limits
        for m in range(1, n+1):
            R[n][m] = R[n][m-1] + (R[n][m-1] - R[n-1][m-1]) / (4**m - 1)
        print_row(R[n])
        if abs(R[n][n-1] - R[n][n]) < eps:
            return R[n][n]
        n += 1

from math import *

# In this example, the error function erf(1) is evaluated.
print romberg(lambda t: 2/sqrt(pi)*exp(-t*t), 0, 1)

Ejemplo[editar]

Como ejemplo, se le integra la función gaussiana en el intervalo ${\displaystyle [O,1]}$, esto es la función error evaluada en 1, cuyo valor es ${\displaystyle {\rm {erf}}(1)\doteq 0.842700792949715}$. Se calculan los elementos de la matriz triangular fila a fila, terminando los cálculos cuando la diferencia entre las dos últimas filas es menor que 1E-8.

 0.77174333
 0.82526296  0.84310283
 0.83836778  0.84273605  0.84271160
 0.84161922  0.84270304  0.84270083  0.84270066
 0.84243051  0.84270093  0.84270079  0.84270079  0.84270079

El resultado en la esquina inferior derecha de la matriz triangular es el resultado correcto con la precisión deseada. Nótese que este resultado se deriva de aproximaciones mucho peores obtenidas con la regla del trapecio mostradas aquí en la primera columna de la matriz triangular.

 regla o método de Simpson (nombrada así en honor de Thomas Simpson) y a veces llamada regla de Kepler es un método de integración numérica que se utiliza para obtener la aproximación de la integral:

${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{6}}\left[f(a)+4f\left({\frac {a+b}{2}}\right)+f(b)\right]}$.

La función f (x) (azul) es aproximada por una función cuadrática P (x) (rojo).

Introducción[editar]

En integración numérica, una forma de aproximar una integral definida en un intervalo [a,b] es mediante la regla del trapecio, es decir, que sobre cada subintervalo en el que se divide [a,b] se aproxima f por un polinomio de primer grado, para luego calcular la integral como suma de las áreas de los trapecios formados en esos subintervalos . El método utilizado para la regla de Simpson sigue la misma filosofía, pero aproximando los subintervalos de f mediante polinomios de segundo grado.

Derivación de la regla de Simpson[editar]

Consideramos el polinomio interpolador de orden dos ${\displaystyle P_{2}(x)}$, que aproxima a la función integrando ${\displaystyle f(x)}$ entre los nodos x₀ = a, x₁ = b y m = (a+b)/2. La expresión de ese polinomio interpolante, expresado a través de la interpolación polinómica de Lagrange es:

${\displaystyle P_{2}(x)=f(a){\frac {(x-m)(x-b)}{(a-m)(a-b)}}+f(m){\frac {(x-a)(x-b)}{(m-a)(m-b)}}+f(b){\frac {(x-a)(x-m)}{(b-a)(b-m)}}.}$

Así, la integral buscada¹​

${\displaystyle I=\int _{a}^{b}f(x)\,dx}$

es equivalente a

${\displaystyle I=\int _{a}^{b}P_{2}(x)\,dx+{\mbox{término error}}={\frac {b-a}{6}}\left[f(a)+4f(m)+f(b)\right]+E(f),}$

donde E(f) es el término de error; por lo tanto, se puede aproximar como:

${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{6}}\left[f(a)+4f(m)+f(b)\right].}$

Error[editar]

El término error E(f), llamado error global, corresponde a¹​

${\displaystyle E(f)=-{\frac {h^{5}}{90}}\,f^{(4)}(\xi ),}$

donde ${\displaystyle h=(b-a)/2}$ y ${\displaystyle \xi }$ pertenece al intervalo [a,b].

Se puede calcular una estimación del error cometido al aproximar la integral mediante este método. Si las cuatro primeras derivadas de f(x) son continuas en el intervalo, entonces el error (en términos absolutos) está acotado como²​

${\displaystyle |E(f)|=\left|\int _{a}^{b}f(x)\,dx-\int _{a}^{b}P_{2}(x)\,dx\right|\leq {\frac {h^{5}}{90}}\,\max _{a\leq \xi \leq b}\left|f^{(4)}(\xi )\right|,}$

donde, de nuevo ${\displaystyle h=(b-a)/2}$ y ${\displaystyle \xi \in [a,b]}$.

Regla de Simpson compuesta[editar]

En el caso de que el intervalo [a,b] no sea lo suficientemente pequeño, el error al calcular la integral puede ser muy grande. Para ello, se recurre a la fórmula compuesta de Simpson. Se divide el intervalo [a,b] en nsubintervalos iguales (con n par), de manera que ${\displaystyle x_{i}=a+ih}$, donde ${\displaystyle h=(b-a)/n}$ para ${\displaystyle i=0,1,...,n}$.

Aplicando la Regla de Simpson a cada subintervalo ${\displaystyle [x_{j-1},x_{j+1}],\ j=1,3,5,...,n-1,}$ tenemos:

${\displaystyle \int _{x_{j-1}}^{x_{j+1}}f(x)\,dx={\frac {x_{j+1}-x_{j-1}}{3}}\left[f(x_{j-1})+4f(x_{j})+f(x_{j+1})]\right.}$

Sumando las integrales de todos los subintervalos, llegamos a que:

${\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {h}{3}}{\bigg [}f(x_{0})+2\sum _{j=1}^{(n/2)-1}f(x_{2j})+4\sum _{j=1}^{n/2}f(x_{2j-1})+f(x_{n}){\bigg ]},}$

El máximo error viene dado por la expresión ${\displaystyle (b-a)\,{\frac {h^{4}}{180}}\,\max _{a\leq \xi \leq b}\left|f^{(4)}(\xi )\right|.}$

Regla de Simpson 3/8 simple[editar]

Esta forma es muy similar a la regla de Simpson clásica, pero se usa polinomios de Lagrange de tercer orden. Se tiene en consideración que ahora el paso ${\displaystyle h={\tfrac {(b-a)}{3}}}$, ya que la función se tabula con cuatro puntos de igual distancia h y formando tres subintervalos. Si x_n+1= x_n+h con x₀=a, se define de la siguiente manera:

${\displaystyle {\begin{alignedat}{2}I&=\int _{a}^{b}f(x)dx\approx \int _{a}^{b}P_{3}(x)dx\\&\approx {\tfrac {3h}{8}}\left[f(x_{0})+3f(x_{1})+3f(x_{2})+f(x_{3})\right]={\tfrac {3h}{8}}\left[f(a)+3f\left({\tfrac {2a+b}{3}}\right)+3f\left({\tfrac {a+2b}{3}}\right)+f(b)\right]\\\end{alignedat}}}$

El error al usar la regla de Simpson de 3/8 se puede obtener usando:

${\displaystyle E(f)=-{\frac {3}{80}}h^{5}f^{(4)}(\xi )}$

donde ${\displaystyle \xi }$ se encuentra dentro del intervalo [a,b].

Regla de Simpson 3/8 compuesta[editar]

Es más exacta que la regla de Simpson 3/8 simple, ya que divide el intervalo de integración en más subintervalos. Se expresa de la siguiente forma:

${\displaystyle I\approx {\frac {3h}{8}}\left[f(x_{0})+3f(x_{1})+3f(x_{2})+2f(x_{3})+3f(x_{4})+3f(x_{5})+2f(x_{6})+\dots +f(x_{n})\right]}$

tomando ${\displaystyle h={\tfrac {(b-a)}{n}}}$ donde n es el número de subintervalos, con la condición de que n sea múltiplo de 3 y que en cada sumatorio se tomen los valores de ${\displaystyle i=i+3}$.

${\displaystyle I\approx {\frac {3h}{8}}\left[f(x_{0})+3\sum _{i=1,4,7\dots }^{n-2}f(x_{i})+3\sum _{i=2,5,8\dots }^{n-1}f(x_{i})+2\sum _{i=3,6,9\dots }^{n-3}f(x_{i})+f(x_{n})\right]}$

Para el cálculo del error, se obtiene la cuarta derivada de la función y tomando en cuenta que ${\displaystyle \xi }$ debe pertenecer al intervalo de integración, se aplica la siguiente fórmula:

${\displaystyle E(f)={\frac {n}{80}}h^{5}f^{(4)}(\xi )}$

Historia[editar]

La fórmula fue utilizada por primera vez por Evangelista Torricelli, pero debe su nombre al matemático Inglés Thomas Simpson. Corresponde a la regla del tonel que Johannes Kepler ya había formulado en 1615.

Sobre la historia de su surgimiento, Kepler informa en la dedicatoria de su publicación posterior: Después de que la primera esposa de Kepler había muerto en Praga en 1611, Kepler se casó nuevamente - en Linz, donde ahora trabajaba - en 1613. Para la boda compró algunos toneles de vino. Puesto ya el vino en la bodega, el vendedor concurrió con una vara de medir y determinó el contenido para todos los barriles sin pensar o calcular, utilizando un mismo método, consistente en que introducía la punta de metal de la vara de medir a través de la piquera , en diagonal hacia los bordes de ambos fondos y la marca en la piquera arrojaba la medida del volumen del contenido. Kepler se sorprendió con aquello de que una diagonal a través del medio del barril pudiera dar una medida sobre el volumen contenido y puso en duda la exactitud de este método, debido a que, por ejemplo, un barril muy bajo que tuviera una base algo más ancha y por eso un volumen contenido mucho menor podría tener el mismo radio a la vista.³​

A raíz de esto, Kepler formuló en 1615 el escrito Nova Stereometria doliorum vinariorum (Nuevo cálculo del contenido de barriles de vino), en el que buscaba métodos verificables para el cálculo del contenido de los toneles de vino. Uno de estos métodos consistió en aproximar la curvatura del barril por una parábola, dado que los cálculos con ayuda de parábolas ya se podían realizar muy exactamente desde Arquímedes.⁴​

Entre otras cosas, Kepler describió en este texto una fórmula para el cálculo de la capacidad (más precisamente, del volumen) de barriles de vino con formas irregulares. Esta fórmula arroja valores exactos para el tronco de la pirámide (incluida la pirámide), la esfera, el paraboloide elíptico, el hiperboloide de una hoja y todas las demás superficies de un cuerpo que pueden ser generadas por secciones planas perpendiculares al eje del cuerpo.

La función f(x) (en azul) es aproximada por la función lineal (en rojo).

En análisis numérico la regla del trapecio es un método de integración, es decir, un método para calcular aproximadamente el valor de una integral definida. La regla se basa en aproximar el valor de la integral de ${\displaystyle f(x)}$ por el de la función lineal, que pasa a través de los puntos ${\displaystyle (a,f(a))}$ y ${\displaystyle (b,f(b))}$. La integral de ésta es igual al área del trapeciobajo la gráfica de la función lineal.

Regla del trapecio Simple[editar]

Para realizar la aproximación por esta regla es necesario usar un polinomio de primer orden, y esta es representada por:

${\displaystyle P_{1}(x)=f(a)+{\frac {f(b)-f(a)}{b-a}}(x-a)}$

Entonces al sustituir en la integral tenemos:

${\displaystyle {\begin{aligned}I&=\int _{a}^{b}f(x)dx\approx \int _{a}^{b}P_{1}(x)\,dx\\&\approx \int _{a}^{b}\left[f(a)+{\frac {f(b)-f(a)}{b-a}}(x-a)\right]\,dx\end{aligned}}}$

Por último al resolver esa integral nos queda:

${\displaystyle \int _{a}^{b}f(x)dx\approx (b-a){\frac {f(a)+f(b)}{2}}}$

Cálculo del error[editar]

El término de error corresponde a:

${\displaystyle E_{t}=-{\frac {1}{12}}f''(\xi )(b-a)^{3}}$

Siendo ${\displaystyle \xi }$ un número perteneciente al intervalo ${\displaystyle [a,b]}$.

Regla del trapecio compuesta[editar]

Ilustración de la regla del trapecio compuesta

La regla del trapecio compuesta o regla de los trapecios es una forma de aproximar una integral definida utilizando n trapecios. En la formulación de este método se supone que f es continua y positiva en el intervalo [a,b]. De tal modo la integral definida ${\displaystyle \int _{a}^{b}f(x)\,dx}$ representa el área de la región delimitada por la gráfica de f y el eje x, desde x=a hasta x=b. Primero se divide el intervalo [a,b] en n subintervalos, cada uno de ancho ${\displaystyle \Delta x=(b-a)/n}$.

Después de realizar todo el proceso matemático se llega a la siguiente fórmula:

${\displaystyle \int _{a}^{b}f(x)\,dx\sim {\frac {h}{2}}[f(a)+2f(a+h)+2f(a+2h)+...+f(b)]}$

Donde ${\displaystyle \textstyle h={\frac {b-a}{n}}}$ y n es el número de divisiones.

La expresión anterior también se puede escribir como:

${\displaystyle \int _{a}^{b}f(x)dx\sim {\frac {b-a}{n}}\left[{\frac {f(a)+f(b)}{2}}+\sum _{k=1}^{n-1}f\left(a+k{\frac {b-a}{n}}\right)\right]}$

El error en esta aproximación se corresponde con :

${\displaystyle -{\frac {(b-a)^{3}}{12n^{2}}}\,f''(\xi )}$

Siendo n el número de subintervalos

Ejemplo[editar]

${\displaystyle \int _{0}^{2}3x\,dx}$ para n=6

Primero se obtiene h, de los límites de la integral que representan a y b y para n=6 queda: ${\displaystyle h={\frac {b-a}{n}}}$ ${\displaystyle ={\frac {2-0}{6}}={\frac {1}{3}}}$.

Y ahora se sustituye en la fórmula

${\displaystyle \int _{a}^{b}f(x)\,dx}$ = ${\displaystyle {\frac {h}{2}}[f(a)+2f(a+h)+2f(a+2h)+...+f(b)]}$

y queda:

${\displaystyle \int _{0}^{2}3x\,dx}$ = ${\displaystyle {\frac {1}{2}}\cdot {\frac {1}{3}}[3(0)+2[3(0+1\cdot {\frac {1}{3}})]+2[3(0+2\cdot {\frac {1}{3}})]+2[3(0+3\cdot {\frac {1}{3}})]+2[3(0+4\cdot {\frac {1}{3}})]+2[3(0+5\cdot {\frac {1}{3}})]+3(2)]=6}$

En este caso no se comete ningún error en el cálculo (el resultado es exacto) porque la función sujeta a integración es lineal.