Editing Romberg's method (section)

== Implementation ==
Here is an example of a computer implementation of the Romberg method (in the [[C (programming language)|C programming language]]):

<syntaxhighlight lang="c">
#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
}
</syntaxhighlight>

Here is an implementation of the Romberg method (in the [[Python (programming language)|Python programming language]]):

<syntaxhighlight lang="python">
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
</syntaxhighlight>