Método de Romberg

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:

${\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

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

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

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.

Referencias

• Richardson, L. F. (1911), «The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differential Equations, with an Application to the Stresses in a Masonry Dam», Philosophical Transactions of the Royal Society of London. Series A 210: pp. 307-357 .
• Romberg, W. (1955), «Vereinfachte numerische Integration», Norske Videnskabers Selskab Forhandlinger (Trondheim) 28 (7): pp. 30-36 .
• Thacher, Jr., Henry C. (julio de 1964), «Remark on Algorithm 60: Romberg integration», Communications of the ACM 7 (7): 420-421 .
• Bauer, F.L.; Rutishauser, H.; Stiefel, E. (1963), «New aspects in numerical quadrature», en Metropolis, N. C., et al., ed., Experimental Arithmetic, high-speed computing and mathematics, Proceedings of Symposia in Applied Mathematics (AMS) (15): pp. 199-218 .
• Bulirsch, Roland; Stoer, Josef (1967), «Handbook Series Numerical Integration. Numerical quadrature by extrapolation», Numerische Mathematik 9: 271-278 .
• Mysovskikh, I.P. (2002), «Romberg method», en Hazewinkel, Michiel, ed., Encyclopaedia of Mathematics, Springer-Verlag, ISBN 1-4020-0609-8 .

Enlaces externos

• ROMBINT -- código para MATLAB (autor: Martin Kacenak)

•   Datos: Q1078449