Romberg's method

<h2 id="method">Method</h2>
<p>Using 
  
    
      
        
          h
          
            n
          
        
        =
        
          
            
              (
              b
              −
              a
              )
            
            
              2
              
                n
              
            
          
        
      
    
    {\textstyle h_{n}={\frac {(b-a)}{2^{n}}}}
  
, the method can be inductively defined by

R
                (
                0
                ,
                0
                )
              
              
                
                =
                
                  h
                  
                    0
                  
                
                (
                f
                (
                a
                )
                +
                f
                (
                b
                )
                )
              
            
            
              
                R
                (
                n
                ,
                0
                )
              
              
                
                =
                
                  
                    
                      1
                      2
                    
                  
                
                R
                (
                n
                −
                1
                ,
                0
                )
                +
                2
                
                  h
                  
                    n
                  
                
                
                  ∑
                  
                    k
                    =
                    1
                  
                  
                    
                      2
                      
                        n
                        −
                        1
                      
                    
                  
                
                f
                (
                a
                +
                (
                2
                k
                −
                1
                )
                
                  h
                  
                    n
                    −
                    1
                  
                
                )
              
            
            
              
                R
                (
                n
                ,
                m
                )
              
              
                
                =
                R
                (
                n
                ,
                m
                −
                1
                )
                +
                
                  
                    
                      1
                      
                        
                          4
                          
                            m
                          
                        
                        −
                        1
                      
                    
                  
                
                (
                R
                (
                n
                ,
                m
                −
                1
                )
                −
                R
                (
                n
                −
                1
                ,
                m
                −
                1
                )
                )
              
            
            
              
              
                
                =
                
                  
                    1
                    
                      
                        4
                        
                          m
                        
                      
                      −
                      1
                    
                  
                
                (
                
                  4
                  
                    m
                  
                
                R
                (
                n
                ,
                m
                −
                1
                )
                −
                R
                (
                n
                −
                1
                ,
                m
                −
                1
                )
                )
              
            
          
        
      
    
    {\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}}}

where 
  
    
      
        n
        ≥
        m
      
    
    {\displaystyle n\geq m}
  
 and 
  
    
      
        m
        ≥
        1
        
      
    
    {\displaystyle m\geq 1\,}
  
. 
In <a href="/facts/Big_O_notation/weFFjSWg">big O notation</a>, the error for <i>R</i>(<i>n</i>, <i>m</i>) is:<a class="footnote-ref" id="fnref:3" href="#fn:3"><sup>3</sup></a> 
  
    
      
        O
        
          (
          
            h
            
              n
            
            
              2
              m
              +
              2
            
          
          )
        
        .
      
    
    {\displaystyle O\left(h_{n}^{2m+2}\right).}

</p><p>The zeroeth extrapolation, <i>R</i>(<i>n</i>, 0), is equivalent to the <a href="/facts/Trapezoidal_rule/KEg8cLGa">trapezoidal rule</a> with 2<i>n</i> + 1 points; the first extrapolation, <i>R</i>(<i>n</i>, 1), is equivalent to <a href="/facts/Simpson%2527s_rule/DpCal7T8">Simpson's rule</a> with 2<i>n</i> + 1 points. The second extrapolation, <i>R</i>(<i>n</i>, 2), is equivalent to <a href="/facts/Boole%2527s_rule/faxq9c1B">Boole's rule</a> with 2<i>n</i> + 1 points. The further extrapolations differ from Newton-Cotes formulas. In particular further Romberg extrapolations expand on Boole's rule in very slight ways, modifying weights into ratios similar as in Boole's rule. In contrast, further Newton-Cotes methods produce increasingly differing weights, eventually leading to large positive and negative weights. This is indicative of how large degree interpolating polynomial Newton-Cotes methods fail to converge for many integrals, while Romberg integration is more stable.
</p><p>By labelling our 
  
    
      
        O
        (
        
          h
          
            2
          
        
        )
      
    
    {\textstyle O(h^{2})}
  
 approximations as 
  
    
      
        
          A
          
            0
          
        
        
          
            (
          
        
        
          
            h
            
              2
              
                n
              
            
          
        
        
          
            )
          
        
      
    
    {\textstyle A_{0}{\big (}{\frac {h}{2^{n}}}{\big )}}
  
 instead of 
  
    
      
        R
        (
        n
        ,
        0
        )
      
    
    {\textstyle R(n,0)}
  
, we can perform Richardson extrapolation with the error formula defined below:

∫
          
            a
          
          
            b
          
        
        f
        (
        x
        )
        
        d
        x
        =
        
          A
          
            0
          
        
        
          
            (
          
        
        
          
            h
            
              2
              
                n
              
            
          
        
        
          
            )
          
        
        +
        
          a
          
            0
          
        
        
          
            (
          
        
        
          
            h
            
              2
              
                n
              
            
          
        
        
          
            
              )
            
          
          
            2
          
        
        +
        
          a
          
            1
          
        
        
          
            (
          
        
        
          
            h
            
              2
              
                n
              
            
          
        
        
          
            
              )
            
          
          
            4
          
        
        +
        
          a
          
            2
          
        
        
          
            (
          
        
        
          
            h
            
              2
              
                n
              
            
          
        
        
          
            
              )
            
          
          
            6
          
        
        +
        ⋯
      
    
    {\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 }

Once we have obtained our 
  
    
      
        O
        (
        
          h
          
            2
            (
            m
            +
            1
            )
          
        
        )
      
    
    {\textstyle O(h^{2(m+1)})}
  
 approximations 
  
    
      
        
          A
          
            m
          
        
        
          
            (
          
        
        
          
            h
            
              2
              
                n
              
            
          
        
        
          
            )
          
        
      
    
    {\textstyle A_{m}{\big (}{\frac {h}{2^{n}}}{\big )}}
  
, we can label them as 
  
    
      
        R
        (
        n
        ,
        m
        )
      
    
    {\textstyle R(n,m)}
  
.
</p><p>When function evaluations are expensive, it may be preferable to replace the polynomial interpolation of Richardson with the rational interpolation proposed by Bulirsch & Stoer (1967).
</p>
<h2 id="a-geometric-example">A geometric example</h2>
<p>To estimate the area under a curve the trapezoid rule is applied first to one-piece, then two, then four, and so on.
</p>

<p>After trapezoid rule estimates are obtained, <a href="/facts/Richardson_extrapolation/46hdzKDA">Richardson extrapolation</a> is applied.
</p>
<ul><li>For the first iteration the two piece and one piece estimates are used in the formula ⁠4 × (more accurate) − (less accurate)/3⁠. The same formula is then used to compare the four piece and the two piece estimate, and likewise for the higher estimates</li>
<li>For the second iteration the values of the first iteration are used in the formula ⁠16 × (more accurate) − (less accurate)/15⁠</li>
<li>The third iteration uses the next power of 4: ⁠64 × (more accurate) − (less accurate)/63⁠ on the values derived by the second iteration.</li>
<li>The pattern is continued until there is one estimate.</li></ul>

<table><tbody><tr><th>Number of pieces</th><th>Trapezoid estimates</th><th>First iteration</th><th>Second iteration</th><th>Third iteration</th></tr><tr><td></td><td></td><td>⁠4 MA − LA/3⁠</td><td>⁠16 MA − LA/15⁠</td><td>⁠64 MA − LA/63⁠</td></tr><tr><td>1</td><td>0</td><td>⁠4×16 − 0/3⁠ = 21.333...</td><td>⁠16×34.667 − 21.333/15⁠ = 35.556...</td><td>⁠64×42.489 − 35.556/63⁠ = 42.599...</td></tr><tr><td>2</td><td>16</td><td>⁠4×30 − 16/3⁠ = 34.666...</td><td>⁠16×42 − 34.667/15⁠ = 42.489...</td><td></td></tr><tr><td>4</td><td>30</td><td>⁠4×39 − 30/3⁠ = 42</td><td></td><td></td></tr><tr><td>8</td><td>39</td><td></td><td></td><td></td></tr></tbody></table>
<h2 id="example">Example</h2>
<p>As an example, the <a href="/facts/Gaussian_function/5bG8Kpp9">Gaussian function</a> is integrated from 0 to 1, i.e. the <a href="/facts/Error_function/Ca2H7pGr">error function</a> erf(1) ≈ 0.842700792949715.  The triangular array is calculated row by row and calculation is terminated if the two last entries in the last row differ less than 10−8.
</p>
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

<p>The result in the lower right corner of the triangular array is accurate to the digits shown.
It is remarkable that this result is derived from the less accurate approximations
obtained by the trapezium rule in the first column of the triangular array.
</p>
<h2 id="implementation">Implementation</h2>
<p>Here is an example of a computer implementation of the Romberg method (in the <a href="/facts/C_(programming_language)/Ky2No763">C programming language</a>):
</p>
#include <stdio.h>
#include <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");
}

/*
INPUT:
(*f) : pointer to the function to be integrated
a    : lower limit
b    : upper limit
max_steps: maximum steps of the procedure
acc  : desired accuracy

OUTPUT:
Rp[max_steps-1]: approximate value of the integral of the function f for x in [a,b] with accuracy 'acc' and steps '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 is previous row, Rc is current row
  double h = b-a; //step size
  Rp[0] = (f(a) + f(b))*h*0.5; // first trapezoidal step

print_row(0, Rp);

for (size_t i = 1; i < max_steps; ++i) {
    h /= 2.;
    double c = 0;
    size_t ep = 1 << (i-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); // compute R(i,j)
    }

// Print ith row of R, R[i,i] is the best estimate so far
    print_row(i, Rc);

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

// swap Rn and Rc as we only need the last row
    double *rt = Rp;
    Rp = Rc;
    Rc = rt;
  }
  return Rp[max_steps-1]; // return our best guess
}

<p>Here is an implementation of the Romberg method (in the <a href="/facts/Python_(programming_language)/YbuGqofa">Python programming language</a>):
</p>
def print_row(i, R):
    """Prints a row of the Romberg table."""
    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):
    """
    Calculates the integral of a function using Romberg integration.

Args:
        f: The function to integrate.
        a: Lower limit of integration.
        b: Upper limit of integration.
        max_steps: Maximum number of steps.
        acc: Desired accuracy.

Returns:
        The approximate value of the integral.
    """
    R1, R2 = [0] * max_steps, [0] * max_steps  # Buffers for storing rows
    Rp, Rc = R1, R2  # Pointers to previous and current rows

h = b - a  # Step size
    Rp[0] = 0.5 * h * (f(a) + f(b))  # First trapezoidal step

print_row(0, Rp)

for i in range(1, max_steps):
        h /= 2.0
        c = 0
        ep = 2 ** (i - 1)
        for j in range(1, ep + 1):
            c += f(a + (2 * j - 1) * h)
        Rc[0] = h * c + 0.5 * Rp[0]  # R(i,0)

for j in range(1, i + 1):
            n_k = 4**j
            Rc[j] = (n_k * Rc[j - 1] - Rp[j - 1]) / (n_k - 1)  # Compute R(i,j)

# Print ith row of R, R[i,i] is the best estimate so far
        print_row(i, Rc)

if i > 1 and abs(Rp[i - 1] - Rc[i]) < acc:
            return Rc[i]

# Swap Rn and Rc for next iteration
        Rp, Rc = Rc, Rp

return Rp[max_steps - 1]  # Return our best guess

<h3>Citations</h3>

<h3>Bibliography</h3>

<ul><li><a href="/facts/Lewis_Fry_Richardson/kKQR6lFu">Richardson, L. F.</a> (1911), "The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differential Equations, with an Application to the Stresses in a Masonry Dam", <i>Philosophical Transactions of the Royal Society A</i>, 210 (459–470): 307–357, <a href="/facts/Bibcode_(identifier)/9HtdQSGB">Bibcode</a>:<a href="https://ui.adsabs.harvard.edu/abs/1911RSPTA.210..307R">1911RSPTA.210..307R</a>, <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1098%2Frsta.1911.0009">10.1098/rsta.1911.0009</a>, <a href="/facts/JSTOR_(identifier)/YTeVmaJ7">JSTOR</a> <a href="https://www.jstor.org/stable/90994">90994</a></li>
<li>Romberg, W. (1955), "Vereinfachte numerische Integration", <i>Det Kongelige Norske Videnskabers Selskab Forhandlinger</i>, 28 (7), Trondheim: 30–36</li>
<li>Thacher Jr., Henry C. (July 1964), "Remark on Algorithm 60: Romberg integration", <i>Communications of the ACM</i>, 7 (7): 420–421, <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1145%2F364520.364542">10.1145/364520.364542</a></li>
<li>Bauer, F.L.; Rutishauser, H.; Stiefel, E. (1963), Metropolis, N. C.; et al. (eds.), "New aspects in numerical quadrature", <i>Experimental Arithmetic, High-speed Computing and Mathematics, Proceedings of Symposia in Applied Mathematics</i> (15), <a href="/facts/American_Mathematical_Society/fr5gotCt">AMS</a>: 199–218</li>
<li>Bulirsch, Roland; Stoer, Josef (1967), <a href="http://www-gdz.sub.uni-goettingen.de/cgi-bin/digbib.cgi?PPN362160546_0009">"Handbook Series Numerical Integration. Numerical quadrature by extrapolation"</a>, <i>Numerische Mathematik</i>, 9: 271–278, <a href="/facts/Doi_(identifier)/muM9Etpq">doi</a>:<a href="https://doi.org/10.1007%2Fbf02162420">10.1007/bf02162420</a></li>
<li>Mysovskikh, I.P. (2002) [1994], <a href="https://www.encyclopediaofmath.org/index.php?title=Romberg_method">"Romberg method"</a>, in Hazewinkel, Michiel (ed.), <i><a href="/facts/Encyclopedia_of_Mathematics/WC6mGtPm">Encyclopedia of Mathematics</a></i>, <a href="/facts/Springer-Verlag/nAesf6nT">Springer-Verlag</a>, <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 1-4020-0609-8</li>
<li>Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), <a href="http://apps.nrbook.com/empanel/index.html?pg=166">"Section 4.3. Romberg Integration"</a>, <i>Numerical Recipes: The Art of Scientific Computing</i> (3rd ed.), New York: Cambridge University Press, <a href="/facts/ISBN_(identifier)/15AdSPa9">ISBN</a> 978-0-521-88068-8</li></ul>

<h2 id="external-links">External links</h2>
<ul><li><a href="http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=34&objectType=file">ROMBINT</a> – code for <a href="/facts/MATLAB/qPjLISCk">MATLAB</a> (author: Martin Kacenak)</li>
<li><a href="http://www.hvks.com/Numerical/webintegration.html">Free online integration tool using Romberg, Fox–Romberg, Gauss–Legendre and other numerical methods</a></li>
<li><a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romberg.html">SciPy implementation of Romberg's method</a></li>
<li><a href="https://github.com/fgasdia/Romberg.jl">Romberg.jl</a> — <a href="/facts/Julia_(programming_language)/AoB0PJ9C">Julia</a> implementation (supporting arbitrary factorizations, <i>not</i> just 
  
    
      
        
          2
          
            n
          
        
        +
        1
      
    
    {\displaystyle 2^{n}+1}
  
 points)</li></ul>

<h2 id="references">References</h2>

<ol>
<li id="fn:1"><p>Romberg 1955 - Romberg, W. (1955), "Vereinfachte numerische Integration", Det Kongelige Norske Videnskabers Selskab Forhandlinger, 28 (7), Trondheim: 30–36 <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li>
<li id="fn:2"><p>Richardson 1911 - Richardson, L. F. (1911), "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 A, 210 (459–470): 307–357, Bibcode:1911RSPTA.210..307R, doi:10.1098/rsta.1911.0009, JSTOR 90994 <a href="https://ui.adsabs.harvard.edu/abs/1911RSPTA.210..307R" target="_blank">https://ui.adsabs.harvard.edu/abs/1911RSPTA.210..307R</a> <a href="#fnref:2" class="footnote-back-ref">↩</a></p></li>
<li id="fn:3"><p>Mysovskikh 2002 - Mysovskikh, I.P. (2002) [1994], "Romberg method", in Hazewinkel, Michiel (ed.), Encyclopedia of Mathematics, Springer-Verlag, ISBN 1-4020-0609-8 <a href="https://www.encyclopediaofmath.org/index.php?title=Romberg_method" target="_blank">https://www.encyclopediaofmath.org/index.php?title=Romberg_method</a> <a href="#fnref:3" class="footnote-back-ref">↩</a></p></li>
</ol>

Romberg's method open-in-new

Romberg's method