Mètode de Romberg

En càlcul numèric, el mètode de Romberg (Romberg 1955) genera una taula triangular que consisteix en estimacions numèriques de la integral definida

${\displaystyle \int _{a}^{b}f(x)\,dx}$

A base d'utilitzar l'extrapolació de Richardson (Richardson 1910) repetidament sobre el mètode trapezial. El mètode de Romberg avalua l'integrand a punts equidistants. L'integrand ha de tenir derivades contínues tot i que es poden obtenir força bons resultats encara que només existeixin unes quantes derivades. Si és possible avaluar l'integrand en punts desigualment espaiats, llavors altres mètodes com la quadratura de Gauss i la quadratura de Clenshaw-Curtis són, en general, més exactes. El mètode va ser descrit i desenvolupat per Werner Romberg en la década de 1950.

Mètode

El mètode es pot definir per inducció així:

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

on

${\displaystyle n\geq 1}$

${\displaystyle m\geq 1}$

${\displaystyle h_{n}={\frac {b-a}{2^{n}}}.}$

En notació de Landau, l'error de R(n,m) és:

${\displaystyle O\left(h_{n}^{2^{m+1}}\right).}$

La primera extrapolació, ${\displaystyle R(n,1)}$, és equivalent al mètode de Simpson amb ${\displaystyle n+2}$ punts.

Quan les avaluacions de la funció són costoses en termes computacionals, pot ser preferible substituir la interpolació polinòmica de Richardson per la interpolació racional proposada per Bulirsch & Stoer (1967).

Implementació del mètode de Romberg en Python

Aquesta és una implementació del mètode de Romberg en Python.

from math import *

def imprimeix_fila(lst):
 print ' '.join('%11.8f' % x for x in lst)

def romberg(f, a, b, eps = 1E-8):
 """Aproxima la integral definida de f des d'«a» fins a «b» pel mètode de Romberg.
 eps és la precisió desitjada."""
 R = [[0.5 * (b - a) * (f(a) + f(b))]] # R[0][0]
 imprimeix_fila(R[0])
 n = 1
 while True:
 h = float(b-a)/2**n
 R.append((n+1)*[None]) # Afegeix una fila buida.
 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)) # per límits adequats
 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)
 imprimeix_fila(R[n])
 if abs(R[n][n-1] - R[n][n]) < eps:
 return R[n][n]
 n += 1

# En aquest exemple s'avalua la funció error erf(1).
print romberg(lambda t: 2/sqrt(pi)*exp(-t*t), 0, 1)

Exemple

Com a exemple s'integra la funció de Gauss des de 0 fins a 1, és a dir la funció error ${\displaystyle {\rm {erf}}(1)\doteq 0.842700792949715}$. La taula triangular es calcula fila per fila i el càlcul s'acaba si els dos útims nombres de l'última fila defereixen menys de 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 resultat de la cantonada de baix a la dreta de la taula triangular és exacte en tots els dígits que es presenten. És notable que aquest resultat s'ha obtingut a partir de les menys exactes aproximacions obtingudes pel mètode trapezial de la primera columna de la taula triangular.

Referències

• Richardson, L. F. «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, 1911, pp. 307–357.
• Romberg, W. «Vereinfachte numerische Integration». Norske Videnskabers Selskab Forhandlinger [Trondheim], 28, 1955, pp. 30–36.
• Thacher, Jr., Henry C. «Remark on Algorithm 60: Romberg integration». Communications of the ACM, 7, 1964, p. 420-421.
• Bauer, F.L.; Rutishauser; Stiefel, E. «New aspects in numerical quadrature». Experimental Arithmetic, high-speed computing and mathematics, Proceedings of Symposia in Applied Mathematics. AMS, 1963, pp. 199–218.
• Bulirsch, Roland; Stoer, Josef «Handbook Series Numerical Integration. Numerical quadrature by extrapolation». Numerische Mathematik, 9, 1967, p. 271–278.
• Mysovskikh, I.P.. Encyclopaedia of Mathematics. Springer-Verlag, 2002. ISBN 1-4020-0609-8. «Romberg method» 

Enllaços externs

• ROMBINT Arxivat 2008-06-09 a Wayback Machine. -- code for MATLAB (author: Martin Kacenak)
• Romberg's method is implemented in Maxima CAS
• ROMBERG^{[Enllaç no actiu]} -- c++ code for romberg integration
• Module for Romberg Integration