Método de Romberg

De Wikipedia, la enciclopedia libre
Saltar a: navegación, búsqueda

En análisis numérico, el Método de Romberg genera una matriz triangular cuyos elementos son estimaciones numéricas de la integral definida siguiente:

 \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í:

R(0,0) = \frac{1}{2} (b-a) (f(a) + f(b))
R(n,0) = \frac{1}{2} R(n-1,0) + h_n \sum_{k=1}^{2^{n-1}} f(a + (2k-1)h_n)
R(n,m) = R(n,m-1) + \frac{1}{4^m-1} (R(n,m-1) - R(n-1,m-1))

o

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

donde

 n \ge 1
 m \ge 1
 h_n = \frac{b-a}{2^n}.

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

 O\left(h_n^{2^{m+1}}\right).

La extrapolación a orden cero R(n,0) es equivalente a la Regla del trapecio con n+2 puntos. a orden uno R(n,1) es equivalente a la Regla de Simpson con 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 [O,1], esto es la función error evaluada en 1, cuyo valor es {\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.

Referencias[editar]

Enlaces externos[editar]