Romberg Integration Scheme¶

Trapezoidal Rule¶

The approximation formula of an integral based on Trapezoidal Rule reads

$$ \begin{align} I(m) \approx h_m \frac{1}{2} (f(a) + f(b)) + h_m \sum_{{\rm interior}} f(x_i) \end{align} $$

We employ a scheme of $2^{m-1}$ intervals ($N=2^{m-1}+1$ grids) and thus

$$ \begin{align} h_m = \frac{b-a}{2^{m-1}}. \end{align} $$

This allows us to compute

$$ \begin{align} I(m) = \frac{1}{2} I(m-1) + h_m \sum_{p=1,3,5,...,2^{m-1}-1} f(a + p \, h_m) \end{align} $$

n-point methods¶

  • n=2 Trapezoidal rule
  • n=3 Simpson's rule
  • n=4 Simpson's 3/8 rule

Assume this trend can be generalized...

The estimate of an integral $\mathcal{I}$ by an n-point method with $2^{m-1}$ intervals is written as

$$ \begin{align} \mathcal{I} &= I(n, m) + C \times h^{2n-2}. \end{align} $$

Writing the same formula, but for $2^{m+1}$ intervals, read

$$ \begin{align} \mathcal{I} = I(n, m+1) + C \times h^{2n-2} \times \frac{1}{2^{2n-2}}. \end{align} $$

Solving for $\mathcal{I}$ and $C$ give

$$ \begin{align} \mathcal{I} \approx \frac{1}{4^{n-1}-1} \, [ 4^{n-1} I(n, m+1) - I(n, m) ]. \end{align} $$

Such an estimate can then be associated as $I(n \rightarrow n+1, m+1)$, i.e., the next order in $n$ with $2^{m}$ intervals:

$$ \begin{align} I(n+1, m+1) = \frac{1}{4^{n-1}-1} \, [ 4^{n-1} I(n, m+1) - I(n, m) ]. \end{align} $$

implementation¶

Construct a matrix $R(i_1, i_2)$ to represent approximation formulae for $I(n=i_2 + 1; m=i_1)$:

We can construct $$ \begin{align} R(i_1, i_2) = \frac{ 4^{i_2-1} R(i_1, i_2-1) - R(i_1-1, i_2-1) }{4^{i_2-1}-1}. \end{align} $$ together with

$$ \begin{align} R(i_1, 1) = \frac{1}{2} R(i_1-1, 1) + h_m \sum_{p=1,3,5,...,2^{m-1}-1} f(a + p \, h_m). \end{align} $$

The process can then be repeated via a triangular motion

(1, 1)

(2, 1) (2, 2)

(3, 1) (3, 2) (3, 3)

...

In [1]:
# Romberg Integration

function romberg(f, a, b; Max=25, tol=1E-15)
        
    Rmat = zeros(Max, Max)
    # initial result
    Rmat[1, 1] = 0.5 * (b-a) * (f(a)+f(b))
    estim = Rmat[1, 1]
      
    for mm in 2:Max
        
        # without duplicated calculations
        dx = (b-a)/2^(mm-1)
        inside = dx * sum([f(a + (2*kk-1)*dx) for kk in 1:2^(mm-2)])
        
        Rmat[mm, 1] = 0.5 * Rmat[mm-1, 1] + inside
        
        for nn in 2:mm
            Rmat[mm, nn] = (4^(nn-1)*Rmat[mm, nn-1]-Rmat[mm-1, nn-1])/(4^(nn-1)-1)
        end
 
        if abs(1-Rmat[mm-1, mm-1]/Rmat[mm, mm]) < tol
            # println(["finish at $mm"])
            break
        else
            estim = Rmat[mm, mm]
        end
            
    end
    
    return estim
    
end
Out[1]:
romberg (generic function with 1 method)
In [2]:
romberg(x->sin(x), 0.0, pi), 2.0
Out[2]:
(2.0000000000000004, 2.0)
In [3]:
using BenchmarkTools
@benchmark ret1 = romberg(x->sin(x), 0.0, pi)
Out[3]:
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.012 μs … 845.362 μs  ┊ GC (min … max):  0.00% … 99.62%
 Time  (median):     1.429 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   1.816 μs ±  14.437 μs  ┊ GC (mean ± σ):  20.81% ±  2.80%

                ▁▅▂ ▁▃▃██▅▅▂                                   
  ▁▁▂▃▃▂▂▂▁▁▁▁▂▇███▇████████▆▅▄▃▃▂▂▂▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  1.01 μs         Histogram: frequency by time        2.16 μs <

 Memory estimate: 6.50 KiB, allocs estimate: 17.

further improvement¶

  1. keep only the last two rows
  2. go as large as needed
In [4]:
function romberg_dynamic(f, a, b; tol=1E-15, Nmax=50)

    # keeping the last two rows of Rmat: R1, R2
    
    # Initialize the Romberg table with a single row
    R1 = [0.5 * (b-a) * (f(a) + f(b))]
    estim = R1[1]
    mm = 1

    while true
        mm += 1
        # Expand the table dynamically
        R2 = zeros(mm)
        
        # Compute the next trapezoidal estimate
        dx = (b-a)/2^(mm-1)
        inside = dx * sum(f(a + (2*kk-1)*dx) for kk in 1:2^(mm-2))
        R2[1] = 0.5 * R1[1] + inside

        # Perform Richardson extrapolation
        for nn in 2:mm
            R2[nn] = (4^(nn-1) * R2[nn-1] - R1[nn-1]) / (4^(nn-1) - 1)
        end

        # Check convergence
        if abs(1 - estim/R2[end]) < tol
            break
        end

        # Update previous row and continue
        R1 = R2
        estim = R2[end]

        if mm > Nmax
            println("$Nmax iterations reached")
            break
        end
        
    end

    return estim
end
Out[4]:
romberg_dynamic (generic function with 1 method)
In [5]:
romberg_dynamic(x->sin(x), 0.0, pi), 2.0
Out[5]:
(1.9999999999999996, 2.0)
In [6]:
using BenchmarkTools
@benchmark ret1 = romberg_dynamic(x->sin(x), 0.0, pi)
Out[6]:
BenchmarkTools.Trial: 10000 samples with 37 evaluations.
 Range (min … max):  915.541 ns … 183.775 μs  ┊ GC (min … max): 0.00% … 99.18%
 Time  (median):     935.811 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   967.089 ns ±   1.829 μs  ┊ GC (mean ± σ):  1.88% ±  0.99%

     ▁█▂                                                         
  ▂▃▇████▅▄▃▃▅▆▅▅▂▂▂▂▂▃▄▃▄▂▂▁▁▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  916 ns           Histogram: frequency by time         1.07 μs <

 Memory estimate: 768 bytes, allocs estimate: 16.