En el análisis numérico , el método de Romberg ( Romberg 1955 ) se utiliza para estimar la integral definida
aplicando la extrapolación de Richardson ( Richardson 1911 ) 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 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 de Gauss y la cuadratura de Clenshaw-Curtis son generalmente más precisos.
El método lleva el nombre de Werner Romberg (1909-2003), quien publicó el método en 1955.
Método
Utilizando
el método puede definirse inductivamente por
o
dónde y . En notación O grande , el error para R ( n , m ) es ( Mysovskikh 2002 ) :
La extrapolación cero, R ( n , 0), es equivalente a la regla trapezoidal con 2 n + 1 puntos; la primera extrapolación, R ( n , 1), es equivalente a la regla de Simpson con 2 n + 1 puntos. La segunda extrapolación, R ( n , 2), es equivalente a la regla de Boole con 2 n + 1 puntos. Otras extrapolaciones difieren de las fórmulas de Newton-Cotes. En particular, otras 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 pesos cada vez más diferentes, lo que eventualmente conduce a grandes pesos positivos y negativos. Esto es indicativo de cómo los métodos polinomiales de Newton-Cotes de interpolación de gran grado no logran converger para muchas integrales, mientras que la integración de Romberg es más estable.
Cuando las evaluaciones de funciones son caras, 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, la regla del trapezoide se aplica primero a una pieza, luego a dos, luego a cuatro, y así sucesivamente.
Una vez obtenidas las estimaciones de la regla trapezoidal, se aplica la extrapolación de Richardson .
- Para la primera iteración, las estimaciones de dos piezas y de 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 la misma 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 precisa) - menos precisa)) / 15
- La tercera iteración usa la siguiente potencia de 4: (64 (más precisa) - menos precisa)) / 63 en los valores derivados de la segunda iteración.
- El patrón continúa hasta que haya una estimación.
Numero de piezas | Estimaciones trapezoidales | Primera iteración | Segunda iteración | Tercera iteración |
---|---|---|---|---|
(4 MA - LA) / 3 * | (16 MA - LA) / 15 | (64 MA - LA) / 63 | ||
1 | 0 | (4 × 16 - 0) / 3 = 21,333 ... | (16 × 34,667 - 21,333) / 15 = 35,556 ... | (64 × 42.489 - 35.556) / 63 = 42.599 ... |
2 | dieciséis | (4 × 30 - 16) / 3 = 34,666 ... | (16 × 42 - 34,667) / 15 = 42,489 ... | |
4 | 30 | (4 × 39 - 30) / 3 = 42 | ||
8 | 39 |
- MA significa más precisa, LA significa menos precisa
Ejemplo
Como ejemplo, la función gaussiana está integrada 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 −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 exacto a 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 ).
#include #include void dump_row ( size_t i , doble * R ) { printf ( "R [% 2zu] =" , i ); para ( tamaño_t j = 0 ; j <= i ; ++ j ) { printf ( "% f" , R [ j ]); } printf ( " \ n " ); }romberg doble ( doble ( * f / * función para integrar * / ) ( doble ), doble / * límite inferior * / a , doble / * límite superior * / b , size_t max_steps , doble / * precisión deseada * / acc ) { doble R1 [ max_steps ], R2 [ max_steps ]; // almacena en búfer el doble * Rp = & R1 [ 0 ], * Rc = & R2 [ 0 ]; // Rp es la fila anterior, Rc es la fila actual double h = ( b - a ); // tamaño de paso Rp [ 0 ] = ( f ( a ) + f ( b )) * h * .5 ; // primer paso trapezoidal dump_row ( 0 , Rp ); para ( size_t i = 1 ; i < max_steps ; ++ i ) { h / = 2 .; doble c = 0 ; size_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 + 0,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) } // Volcado ith fila de R, R [i, i] es la mejor estimación hasta ahora dump_row ( i , Rc ); if ( i > 1 && fabs ( Rp [ i -1 ] - Rc [ i ]) < acc ) { return Rc [ i -1 ]; } // intercambia Rn y Rc ya que solo necesitamos la última fila double * rt = Rp ; Rp = Rc ; Rc = rt ; } return Rp [ max_steps -1 ]; // devuelve nuestra mejor suposición }
________________________________________________________
A continuación se muestra un ejemplo de una implementación informática del método Romberg (en el lenguaje de programación javascript ).
function auto_integrator_trap_romb_hnm ( func , a , b , nmax , tol_ae , tol_rae ) // ENTRADAS // func = integrando // a = límite inferior de integración // b = límite superior de integración // nmax = número de particiones, n = 2 ^ nmax // tol_ae = error máximo absoluto aproximado aceptable (debe ser> = 0) // tol_rae = máximo error absoluto absoluto aproximado aceptable (debe ser> = 0) // SALIDAS // integ_value = valor estimado de integral{ if ( typeof a ! == 'number' ) { throw new TypeError ( ' debe ser un número' ); } if ( typeof b ! == 'number' ) { throw new TypeError ( ' debe ser un número' ); } if (( ! Number . isInteger ( nmax )) || ( nmax < 1 )) { lanzar nuevo TypeError ( ' debe ser un número entero mayor o igual que uno.' ); } if (( typeof tol_ae ! == 'number' ) || ( tol_ae < 0 )) { throw new TypeError ( ' debe ser un número mayor o igual a cero' ); } if (( typeof tol_rae ! == 'number' ) || ( tol_rae <= 0 )) { throw new TypeError ( ' debe ser un número mayor o igual a cero' ); } var h = b - a // inicializa la matriz donde se almacenan los valores de la integralvar Romb = []; // filas para ( var i = 0 ; i < nmax + 1 ; i ++ ) { Romb . empujar ([]); para ( var j = 0 ; j < nmax + 1 ; j ++ ) { Romb [ i ]. empujar ( matemáticas . bignumber ( 0 )); } }// calculando el valor con la regla trapezoidal de 1 segmento Romb [ 0 ] [ 0 ] = 0.5 * h * ( func ( a ) + func ( b )) var integ_val = Romb [ 0 ] [ 0 ]for ( var i = 1 ; i <= nmax ; i ++ ) // actualizando el valor con el doble del número de segmentos // usando solo los valores donde deben calcularse // Ver https://autarkaw.org / 2009/02/28 / una-fórmula-eficiente-para-un-integrador-automático-basado-en-la-regla-trapezoidal / { h = 0.5 * h var integ = 0 for ( var j = 1 ; j <= 2 ** i - 1 ; j + = 2 ) { var integ = integ + func ( a + j * h ) }Romb [ i ] [ 0 ] = 0.5 * Romb [ i - 1 ] [ 0 ] + integ * h // Usando el método Romberg para calcular el siguiente valor extrapolable // Ver https://young.physics.ucsc.edu/115/ romberg.pdf para ( k = 1 ; k <= i ; k ++ ) { var addterm = Romb [ i ] [ k - 1 ] - Romb [ i - 1 ] [ k - 1 ] addterm = addterm / ( 4 * * k - 1.0 ) Romb [ i ] [ k ] = Romb [ i ] [ k - 1 ] + addterm// Calculando el error absoluto aproximado var Ea = math . abs ( Romb [ i ] [ k ] - Romb [ i ] [ k - 1 ])// Calculando el error aproximado relativo absoluto var epsa = math . abs ( Ea / Romb [ i ] [ k ]) * 100,0// Asignar el valor más reciente a la variable de retorno integ_val = Romb [ i ] [ k ]// devolviendo el valor si se cumple alguna de las tolerancias if (( epsa < tol_rae ) || ( Ea < tol_ae )) { return ( integ_val ) } } } // devolviendo el último valor calculado de la integral si se cumple o no la tolerancia return ( integ_val ) }
Referencias
- Richardson, LF (1911), "La solución aritmética aproximada por diferencias finitas de problemas físicos que involucran ecuaciones diferenciales, con una aplicación a las tensiones en una presa de mampostería", Transacciones filosóficas de la 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 , Trondheim, 28 (7): 30–36
- Thacher Jr., Henry C. (julio de 1964), "Observación sobre el algoritmo 60: integración de Romberg" , Comunicaciones del ACM , 7 (7): 420–421, doi : 10.1145 / 364520.364542
- Bauer, FL; 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 , AMS (15): 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
- Presione, 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