El método de Romberg

En el análisis numérico , se utiliza el método de Romberg ^{[1] 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 espaciados de manera uniforme. El integrando debe tener derivadas continuas, aunque se pueden obtener resultados bastante buenos si solo existen unas pocas derivadas. Si es posible evaluar el integrando en puntos espaciados de manera desigual, entonces otros métodos como la cuadratura gaussiana y la cuadratura de Clenshaw-Curtis son generalmente más precisos. $$estilo de visualización int _{a}^{b}f(x)\,dx}$$

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

Usando , el método puede definirse inductivamente por donde y . En notación O grande , el error para R ( n ,  m ) es: ^[3]${\textstyle h_{n}={\frac {(ba)}{2^{n}}}}$$${\displaystyle {\begin{aligned}R(0,0)&=h_{1}(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))\\&={\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) , es equivalente a la regla de Boole con 2 ⁿ + 1 puntos. Las extrapolaciones posteriores difieren de las fórmulas de Newton-Cotes. En particular, las extrapolaciones de Romberg posteriores amplían la regla de Boole de formas muy leves, modificando los pesos en proporciones similares a las de la regla de Boole. Por el contrario, los métodos de Newton-Cotes posteriores producen pesos cada vez más diferentes, lo que finalmente conduce a pesos positivos y negativos grandes. Esto es indicativo de cómo los métodos de Newton-Cotes de polinomios de interpolación de gran grado no convergen 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 polinomial de Richardson con la interpolación racional propuesta por Bulirsch y Stoer (1967).

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.

Después de obtener las estimaciones de la regla trapezoidal, 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 se utilizan los valores de la primera iteración en la fórmula16 × (más preciso) − (menos preciso)/15⁠
• La tercera iteración utiliza 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 sola estimación.

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 menos de 10 ⁻⁸ .

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 que aparece en la esquina inferior derecha de la matriz triangular es exacto a los dígitos que se muestran. Es notable que este resultado se derive de las aproximaciones menos exactas obtenidas con la regla del trapecio en la primera columna de la matriz triangular.

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 ( tamaño_t i , double * R ) { printf ( "R[%2zu] = " , i ); for ( tamaño_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 
acc : precisión deseada

SALIDA: 
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'. 
*/ 
double romberg ( double ( * f )( double ), double a , double b , size_t max_steps , double acc ) { double R1 [ max_steps ], R2 [ max_steps ]; // buffers double * 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           

     
         
      
       

  imprimir_fila ( 0 , Rp ); 

  para ( tamaño_t i = 1 ; i < máx_pasos ; ++ i ) { h /= 2. ; doble c = 0 ; tamaño_t ep = 1 << ( i -1 ); //2^(n-1) para ( tamaño_t j = 1 ; j <= ep ; ++ j ) { c += f ( a + ( 2 * j -1 ) * h ); } Rc [ 0 ] = h * c + .5 * Rp [ 0 ]; // R(i,0)         
      
       
          
             
            
    
         

    para ( tamaño_t j = 1 ; j <= i ; ++ j ) { doble n_k = pow ( 4 , j ); Rc [ j ] = ( n_k * Rc [ j -1 ] - Rp [ j -1 ]) / ( n_k -1 ); // calcular R(i,j) }         
          
             
    

    // Imprimir la i-ésima fila de R, R[i,i] es la mejor estimación hasta el momento 
print_row ( i , Rc );     

    si ( i > 1 && fabs ( Rp [ i -1 ] - Rc [ i ]) < acc ) { devolver Rc [ i ]; }        
       
    

    // intercambiamos Rn y Rc ya que solo necesitamos la última fila 
double * rt = Rp ; Rp = Rc ; Rc = rt ; } return Rp [ max_steps -1 ]; // devolvemos nuestra mejor estimació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 de Romberg.""" print ( f "R[ { i : 2d } ]=" , end = " ) for j in range ( i + 1 ): print ( f " { R [ j ] : f } =" , end = " ) print ()  
   
       
     
  

def  romberg ( f ,  a ,  b ,  max_steps ,  acc ): 
"""   Calcula la integral de una función utilizando 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   # Buffers para almacenar filas 
  Rp ,  Rc  =  R1 ,  R2   # Punteros a 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(i,j)

    # Imprimir la i-ésima fila 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 : 
      devuelve  Rc [ i ]

    # Intercambia Rn y Rc para la siguiente iteración 
    Rp ,  Rc  =  Rc ,  Rp

  devuelve  Rp [ max_steps  -  1 ]   # Devuelve nuestra mejor estimación

• 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 del método de Romberg en SciPy
• Romberg.jl: implementación de Julia (compatible con factorizaciones arbitrarias, no solo puntos)${\displaystyle 2^{n}+1}$