Método de Romberg

En análisis numérico , el método de Romberg ^{[ 1 ]} se utiliza para estimar la integral definida.$${\displaystyle \int _{a}^{b}f(x)\,dx}$$mediante la aplicación repetida de la extrapolación de Richardson ^{[ 2 ]} en 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 igualmente espaciados. 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 desigualmente espaciados, entonces otros métodos como la cuadratura gaussiana y la cuadratura de Clenshaw-Curtis suelen ser más precisos.

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

Usando${\textstyle h_{n}={\frac {(b-a)}{2^{n+1}}}}$, el método puede definirse inductivamente por $${\displaystyle {\begin{aligned}R(0,0)&=h_{0}(f(a)+f(b))\\R(n,0)&={\tfrac {1}{2}}R(n{-}1,\,0)+2h_{n}\sum _{k=1}^{2^{n-1}}f(a+(2k-1)h_{n-1})\\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}}}$$ dónde${\displaystyle n\geq m}$y${\displaystyle m\geq 1\,}$. Obsérvese que las dos primeras ecuaciones corresponden a la primera columna, y estas son las fórmulas asociadas a la regla trapezoidal compuesta. En notación O grande , el error para R ( n , m ) es: ^{[ 3 ]}${\displaystyle O{\left(h_{n}^{2m+2}\right)}.}$

La extrapolación cero, R ( n , 0) ^, es equivalente a la regla trapezoidal con 2n + 1 puntos; la primera extrapolación, R ( n , 1) , es equivalente a la regla de Simpson con 2n + 1 puntos. La segunda extrapolación, R ( n , 2) , es equivalente a ^la regla de Boole con 2n + 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 forma muy sutil, modificando los pesos en razones 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, llegando finalmente a pesos positivos y negativos grandes. Esto indica que los métodos de Newton-Cotes polinomiales de interpolación de alto grado no convergen para muchas integrales, mientras que la integración de Romberg es más estable.

Al etiquetar nuestro${\textstyle O(h^{2})}$aproximaciones como${\textstyle A_{0}{\big (}{\frac {h}{2^{n}}}{\big )}}$en lugar de${\textstyle R(n,0)}$Podemos realizar la extrapolación de Richardson con la fórmula de error definida a continuación: $${\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 }$$ Una vez que hayamos obtenido nuestro${\textstyle O(h^{2(m+1)})}$aproximaciones${\textstyle A_{m}{\big (}{\frac {h}{2^{n}}}{\big )}}$, podemos etiquetarlos como${\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) .

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

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

• Para la primera iteración, las estimaciones de dos piezas y una pieza se utilizan en la fórmula 4 × (más precisa) − (menos precisa) / 3. La misma fórmula se utiliza luego para comparar la estimación de cuatro piezas y la de dos piezas, y de igual manera 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 utiliza la siguiente potencia de 4: 64 × ( más preciso) − (menos preciso) / 63 sobre los valores derivados por la segunda iteración.
• El proceso 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,842 700 792 949 715 . La matriz triangular se calcula fila por fila y el cálculo se termina si las dos últimas entradas de la última fila difieren menos de 10 ^{− 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 que se muestra en la esquina inferior derecha de la matriz triangular coincide con los dígitos indicados. Es notable que este resultado se derive de las aproximaciones menos precisas obtenidas mediante la regla del trapecio en la primera columna de la matriz triangular.

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

#include <stdio.h>#include <math.h>void print_row ( size_t i , double * R ) {printf ( "R[%2zu] = " , i );para ( tamaño_t j = 0 ; j <= i ; ++ j ) {printf ( "%f " , R [ j ]);}printf ( " \n " );}/*APORTE:(*f): puntero a la función que se va a integrara: límite inferiorb: límite superiormax_steps: número máximo de pasos del procedimientoacc: precisión deseadaPRODUCCIÓN: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 [ pasos_máximos ], R2 [ pasos_máximos ]; // búferesdouble * Rp = & R1 [ 0 ], * Rc = & R2 [ 0 ]; // Rp es la fila anterior, Rc es la fila actualdoble h = b - a ; //tamaño del pasoRp [ 0 ] = ( f ( a ) + f ( b )) * h * 0.5 ; // primer paso trapezoidalprint_row ( 0 , Rp );para ( tamaño_t i = 1 ; i < max_steps ; ++ 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)}// Imprime la i-ésima fila de R, R[i,i] es la mejor estimación hasta el momento.imprimir_fila ( 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.doble * rt = Rp ;Rp = Rc ;Rc = rt ;}return Rp [ max_steps - 1 ]; // devuelve nuestra mejor estimación}

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

import numpy as npfrom math import erf , sqrt , pidef imprimir_fila ( fila ):# Imprime una fila de la tabla de Romberg en un formato numérico de ancho fijo.# Esto hace que el patrón de convergencia sea fácil de leer columna por columna.print ( " " . join ( f " { x : 11.8f } " for x in row ))def romberg ( f , a , b , eps = 1e-8 , max_iter = 20 ):""" Aproxime la integral definida de f de a a b utilizando la integración de Romberg. Parámetros ---------- f: invocable Función a integrar. Idealmente debería aceptar matrices NumPy para que Las evaluaciones en muchos puntos pueden ser vectorizadas. a, b: flotante Límites inferior y superior de integración. eps: flotante, opcional Tolerancia de parada deseada. El algoritmo se detiene cuando dos sucesivos Las estimaciones diagonales de Romberg son suficientemente precisas. max_iter: int, opcional Número máximo de niveles de refinamiento permitidos. Devoluciones ------- flotar La mejor estimación de Romberg de la integral. """# R[n, m] contendrá el n-ésimo refinamiento y la m-ésima extrapolación de Richardson.R = np.zeros ( ( max_iter + 1 , max_iter + 1 ) )# Caso base: regla trapezoidal utilizando solo los puntos extremos a y b.# Esta es la primera y más burda aproximación a la integral.R [ 0 , 0 ] = ( b - a ) * ( f ( a ) + f ( b )) / 2# Imprime la primera fila, que contiene solo la estimación trapezoidal inicial.print_row ( R [ 0 , : 1 ])# Refinar la aproximación nivel por nivel.# Cada nuevo nivel reduce a la mitad el tamaño del paso y agrega los valores de la función en el# Nuevos puntos medios introducidos por la subdivisión.para n en rango ( 1 , max_iter + 1 ):# Tamaño del paso en el nivel de refinamiento n.# Dado que el intervalo se divide en 2^n subintervalos, el ancho de la malla es:h_n1 = ( b - a ) / ( 2 ** n ) # h_{n-1}# Generar los nuevos puntos introducidos en este nivel de refinamiento.# Estos son los múltiplos impares de h en relación con a:# a + h, a + 3h, a + 5h, ..., a + (2^n - 1)h## Hay 2^(n-1) puntos de este tipo, y NumPy los genera de manera eficiente.x = a + h_n1 * np . organizar ( 1 , 2 ** n , 2 )# Actualizar la estimación trapezoidal utilizando la relación recursiva de Romberg:.R [ n , 0 ] = R [ n - 1 , 0 ] / 2 + h_n1 * np . sum ( f ( x ))# Extrapolación de Richardson:para m en rango ( 1 , n + 1 ):R [ n , m ] = R [ n , m - 1 ] + ( R [ n , m - 1 ] - R [ n - 1 , m - 1 ]) / ( 4 ** m - 1 )# Imprime la fila completada hasta el elemento diagonal actual.print_row ( R [ n , : n + 1 ])# Criterio de parada:# Si los dos últimos valores relacionados con la diagonal están suficientemente cerca,# Aceptamos la estimación extrapolada más reciente.si abs ( R [ n , n ] - R [ n , n - 1 ]) < eps :devolver R [ n , n ]# Si el método no converge dentro de max_iter niveles, genere un error.generar RuntimeError ( "El método Romberg no convergió dentro de max_iter." )# Ejemplo de uso:# Integrar la función correspondiente a erf(1):# erf(1) = 2/sqrt(pi) * integral_0^1 exp(-t^2) dt#f = lambda t : 2 / np . sqrt ( np . pi ) * np . exp ( - t * t )# Calcula e imprime el resultado.print ( "Aproximación =" , romberg ( f , 0.0 , 1.0 ))print ( "Valor exacto =" , erf ( 1.0 ))

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