El método de Romberg.

Método de integración numérica

En análisis numérico , el método de Romberg ^[1] se utiliza para estimar la integral definida aplicando la extrapolación de Richardson ^[2] repetidamente sobre la regla del trapecio o la regla del rectángulo (regla del punto medio). Las estimaciones generan una matriz triangular . El método de Romberg es una fórmula de Newton-Cotes : evalúa el integrando en puntos equidistantes. El integrando debe tener derivadas continuas, aunque se pueden obtener resultados bastante buenos si sólo existen unas pocas derivadas. Si es posible evaluar el integrando en puntos espaciados desigualmente, entonces otros métodos, como la cuadratura gaussiana y la cuadratura de Clenshaw-Curtis, son generalmente más precisos. $${\displaystyle \int _{a}^{b}f(x)\,dx}$$

El método lleva el nombre de Werner Romberg , quien lo publicó en 1955.

Método

Usando , el método se puede definir inductivamente por dónde y . En notación O grande , el error para R ( n ,  m ) es: ^[3] ${\textstyle h_{n}={\frac {(b-a)}{2^{n}}}}$ $${\displaystyle {\begin{aligned}R(0,0)&=h_{1}(f(a)+f(b))\\R(n,0)&={\tfrac {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)+{\tfrac {1}{4^{m}-1}}(R(n,m-1)-R(n-1,m-1))\\&={\frac {1}{4^{m}-1}}(4^{m}R(n,m-1)-R(n-1,m-1))\end{aligned}}}$$ ${\displaystyle n\geq m}$ ${\displaystyle m\geq 1\,}$ ${\displaystyle O\left(h_{n}^{2m+2}\right).}$

La extrapolación cero, R ( n , 0) , es equivalente a la regla trapezoidal con 2 ⁿ + 1 puntos; la primera extrapolación, R ( n , 1) , es equivalente a la regla de Simpson con 2 ⁿ + 1 puntos. La segunda extrapolación, R ( n , 2) , equivale a la regla de Boole con 2 ⁿ + 1 puntos. Las extrapolaciones adicionales difieren de las fórmulas de Newton-Cotes. En particular, nuevas extrapolaciones de Romberg amplían la regla de Boole de manera muy leve, modificando los pesos en proporciones similares a las de la regla de Boole. Por el contrario, otros métodos de Newton-Cotes producen ponderaciones cada vez más diferentes, lo que eventualmente conduce a ponderaciones positivas y negativas grandes. Esto es indicativo de cómo los métodos de Newton-Cotes con polinomios de interpolación de gran grado no logran converger para muchas integrales, mientras que la integración de Romberg es más estable.

Al etiquetar nuestras aproximaciones como en lugar de , podemos realizar la extrapolación de Richardson con la fórmula de error definida a continuación: Una vez que hayamos obtenido nuestras aproximaciones , podemos etiquetarlas como . ${\textstyle O(h^{2})}$ ${\textstyle A_{0}{\big (}{\frac {h}{2^{n}}}{\big )}}$ ${\textstyle R(n,0)}$ $${\displaystyle \int _{a}^{b}f(x)\,dx=A_{0}{\bigg (}{\frac {h}{2^{n}}}{\bigg )}+a_{0}{\bigg (}{\frac {h}{2^{n}}}{\bigg )}^{2}+a_{1}{\bigg (}{\frac {h}{2^{n}}}{\bigg )}^{4}+a_{2}{\bigg (}{\frac {h}{2^{n}}}{\bigg )}^{6}+\cdots }$$ ${\textstyle O(h^{2(m+1)})}$ ${\textstyle A_{m}{\big (}{\frac {h}{2^{n}}}{\big )}}$ ${\textstyle R(n,m)}$

Cuando las evaluaciones de funciones son costosas, puede ser preferible reemplazar la interpolación polinómica de Richardson con la interpolación racional propuesta por Bulirsch y Stoer (1967).

Un ejemplo geométrico

Para estimar el área bajo una curva, se aplica la regla del trapezoide primero a una pieza, luego a dos, luego a cuatro, y así sucesivamente.

Una pieza. Tenga en cuenta que, dado que comienza y termina en cero, esta aproximación produce un área cero.

Dos piezas

cuatro piezas

Ocho piezas

Una vez obtenidas las estimaciones de la regla del trapezoide, se aplica la extrapolación de Richardson .

• Para la primera iteración, se utilizan las estimaciones de dos piezas y de una pieza en la fórmula ⁠4 × (más preciso) − (menos preciso)/3⁠ . Luego se utiliza la misma fórmula para comparar la estimación de cuatro piezas y la de dos piezas, y lo mismo para las estimaciones más altas.
• Para la segunda iteración, los valores de la primera iteración se utilizan en la fórmula ⁠16 × (más preciso) − (menos preciso)/15⁠
• La tercera iteración usa la siguiente potencia de 4: ⁠64 × (más preciso) − (menos preciso)/63⁠ sobre los valores derivados de la segunda iteración.
• El patrón continúa hasta que haya una estimación.

Ejemplo

Como ejemplo, la función gaussiana se integra de 0 a 1, es decir, la función de error erf(1) ≈ 0,842700792949715. La matriz triangular se calcula fila por fila y el cálculo finaliza si las dos últimas entradas de la última fila difieren en menos de 10 ⁻⁸ .

0.771743330,82526296 0,843102830,83836778 0,84273605 0,842711600,84161922 0,84270304 0,84270083 0,842700660,84243051 0,84270093 0,84270079 0,84270079 0,84270079

El resultado en la esquina inferior derecha de la matriz triangular tiene una precisión de los dígitos que se muestran. Es notable que este resultado se derive de las aproximaciones menos precisas obtenidas por la regla del trapecio en la primera columna de la matriz triangular.

Implementación

A continuación se muestra un ejemplo de una implementación informática del método Romberg (en el lenguaje de programación C ):

#incluir <stdio.h> #incluir <math.h>  void print_row ( size_t i , double * R ) { printf ( "R[%2zu] =" , i ); for ( size_t j = 0 ; j <= i ; ++ j ) { printf ( "%f" , R [ j ]); } printf ( " \n " ); }                     /* ENTRADA: (*f) : puntero a la función a integrar a : límite inferior b : límite superior max_steps: pasos máximos del procedimiento según : precisión deseadaSALIDA: Rp[max_steps-1]: valor aproximado de la integral de la función f para x en [a,b] con precisión 'acc' y pasos 'max_steps'. */ doble romberg ( doble ( * f )( doble ), doble a , doble b , tamaño_t max_steps , doble acc ) { doble R1 [ max_steps ], R2 [ max_steps ]; // buffers dobles * Rp = & R1 [ 0 ], * Rc = & R2 [ 0 ]; // Rp es la fila anterior, Rc es la fila actual double h = b - a ; //tamaño del paso Rp [ 0 ] = ( f ( a ) + f ( b )) * h * 0.5 ; // primer paso trapezoidal                                   print_row ( 0 , Rp );  for ( tamaño_t i = 1 ; i < max_steps ; ++ i ) { h /= 2 .; doble c = 0 ; tamaño_t ep = 1 << ( yo -1 ); //2^(n-1) for ( size_t j = 1 ; j <= ep ; ++ j ) { c += f ( a + ( 2 * j -1 ) * h ); } Rc [ 0 ] = h * c + .5 * Rp [ 0 ]; // R(i,0)                                                for ( size_t j = 1 ; j <= i ; ++ j ) { double n_k = pow ( 4 , j ); Rc [ j ] = ( n_k * Rc [ j -1 ] - Rp [ j -1 ]) / ( n_k -1 ); // calcula R(i,j) }                        // Imprime la fila i de R, R[i,i] es la mejor estimación hasta el momento print_row ( i , Rc );   if ( i > 1 && fabs ( Rp [ i -1 ] - Rc [ i ]) < acc ) { return Rc [ i ]; }            // intercambia Rn y Rc ya que solo necesitamos la última fila double * rt = Rp ; Rp = Rc ; Rc = rt ; } devolver Rp [ max_steps -1 ]; // devolvemos nuestra mejor suposición }              

Aquí hay una implementación del método Romberg (en el lenguaje de programación Python ):

def  print_row ( i ,  R ): """Imprime una fila de la tabla Romberg.""" print ( f "R[ { i : 2d } ] = " , end = "" ) para j en el rango ( i + 1 ): imprimir ( f " { R [ j ] : f } " , fin = "" ) imprimir ()            def  romberg ( f ,  a ,  b ,  max_steps ,  acc ): """  Calcula la integral de una función usando la integración de Romberg.  Args:  f: La función a integrar.  a: Límite inferior de integración.  b: Límite superior de integración.  max_steps: Número máximo de pasos.  acc: Precisión deseada. Devuelve:  El valor aproximado de la integral.  """  R1 ,  R2  =  [ 0 ]  *  max_steps ,  [ 0 ]  *  max_steps  # Búfers para almacenar filas  Rp ,  Rc  =  R1 ,  R2  # Punteros a las filas anteriores y actuales h  =  b  -  a  # Tamaño del paso  Rp [ 0 ]  =  0,5  *  h  *  ( f ( a )  +  f ( b ))  # Primer paso trapezoidal imprimir_fila ( 0 ,  Rp ) para  i  en el  rango ( 1 ,  max_steps ):  h  /=  2.  c  =  0  ep  =  1  <<  ( i  -  1 )  # 2^(i-1)  para  j  en el  rango ( 1 ,  ep  +  1 ):  c  + =  f ( a  +  ( 2  *  j  -  1 )  *  h )  Rc [ 0 ]  =  h  *  c  +  0.5  *  Rp [ 0 ]  # R(i,0) para  j  en  el rango ( 1 ,  i  +  1 ):  n_k  =  4 ** j  Rc [ j ]  =  ( n_k  *  Rc [ j  -  1 ]  -  Rp [ j  -  1 ])  /  ( n_k  -  1 )  # Calcular R( yo,j) # Imprimir la fila i de R, R[i,i] es la mejor estimación hasta el momento  print_row ( i ,  Rc ) si  i  >  1  y  abs ( Rp [ i  -  1 ]  -  Rc [ i ])  <  acc :  return  Rc [ i ] # Intercambiar Rn y Rc para la siguiente iteración  Rp ,  Rc  =  Rc ,  Rp return  Rp [ max_steps  -  1 ]  # Devuelve nuestra mejor suposición

Referencias

Citas

1. Romberg 1955
2. Richardson 1911
3. Mysovskikh 2002

Bibliografía

• Richardson, LF (1911), "La solución aritmética aproximada mediante diferencias finitas de problemas físicos que involucran ecuaciones diferenciales, con una aplicación a las tensiones en una presa de mampostería", Philosophical Transactions of the Royal Society A , 210 (459–470): 307 –357, doi :10.1098/rsta.1911.0009, JSTOR  90994
• Romberg, W. (1955), "Vereinfachte numerische Integration", Det Kongelige Norske Videnskabers Selskab Forhandlinger , 28 (7), Trondheim: 30–36
• Thacher Jr., Henry C. (julio de 1964), "Observación sobre el algoritmo 60: integración de Romberg", Comunicaciones de la ACM , 7 (7): 420–421, doi : 10.1145/364520.364542
• Bauer, Florida; Rutishauser, H.; Stiefel, E. (1963), Metropolis, Carolina del Norte; et al. (eds.), "Nuevos aspectos en cuadratura numérica", Aritmética experimental, computación de alta velocidad y matemáticas, Actas de simposios en matemáticas aplicadas (15), AMS : 199–218
• Bulirsch, Roland; Stoer, Josef (1967), "Handbook Series Numerical Integration. Cuadratura numérica por extrapolación", Numerische Mathematik , 9 : 271–278, doi :10.1007/bf02162420
• Mysovskikh, IP (2002) [1994], "Método Romberg", en Hazewinkel, Michiel (ed.), Encyclopedia of Mathematics , Springer-Verlag , ISBN 1-4020-0609-8
• Prensa, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), "Sección 4.3. Integración de Romberg", Recetas numéricas: el arte de la informática científica (3.ª ed.), Nueva York: Cambridge University Press, ISBN 978-0-521-88068-8

enlaces externos

• ROMBINT – código para MATLAB (autor: Martin Kacenak)
• Herramienta de integración en línea gratuita que utiliza Romberg, Fox–Romberg, Gauss–Legendre y otros métodos numéricos
• Implementación SciPy del método de Romberg
• Romberg.jl: implementación de Julia (que admite factorizaciones arbitrarias, no solo puntos) ${\displaystyle 2^{n}+1}$