Romberg Extrapolation by James Pate Williams, Jr.

/*
Author:  Pate Williams (c) January 22, 1995
The following program is a translation of the FORTRAN
subprogram found in Elementary Numerical Analysis by
S. D. Conte and Carl de Boor pages 343-344. The program
uses Romberg extrapolation to find the integral of a
function.
*/

#include "stdafx.h"
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>

typedef long double real;

real t[10][10];

real f(real x)
{
	return(expl(-x * x));
}

real romberg(real a, real b, int start, int row)
{
	int i, k, m;
	real h, ratio, sum;

	m = start;
	h = (b - a) / m;
	sum = 0.5 * (f(a) + f(b));
	if (m > 1)
		for (i = 1; i < m; i++)
			sum += f(a + i * h);
	t[1][1] = h * sum;
	if (row < 2) return(t[1][1]);
	for (k = 2; k <= row; k++)
	{
		h = 0.5 * h;
		m *= 2;
		sum = 0.0;
		for (i = 1; i <= m; i += 2)
			sum += f(a + h * i);
		t[k][1] = 0.5 * t[k - 1][1] + sum * h;
		for (i = 1; i < k; i++)
		{
			t[k - 1][i] = t[k][i] - t[k - 1][i];
			t[k][i + 1] = t[k][i] - t[k - 1][i] /
				(powl(4.0, (real)i) - 1.0);
		}
	}
	if (row < 3) return(t[2][2]);
	for (k = 1; k <= row - 2; k++)
		for (i = 1; i <= k; i++)
		{
			if (t[k + 1][i] == 0.0)
				ratio = 0.0;
			else
				ratio = t[k][i] / t[k + 1][i];
			t[k][i] = ratio;
		}
	return(t[row][row - 1]);
}

int main(void)
{
	long double experimental = 0.7468241328124271;
	int row = 8;

	printf("Romberg integration of f(x) = exp(- x * x)\n");
	printf("Internet value = 0.7468241328124271\n");
	printf("from website = https://www.integral-calculator.com/\n");
	printf("Approximation\t\tPercent Difference\tSteps\tRuntime (s)\n");
	for (long steps = 100000; steps <= 900000; steps += 100000)
	{
		clock_t start = clock();
		long double trial = romberg(0.0, 1.0, steps, row);
		clock_t endTm = clock();
		long double runtime = ((long double)endTm - start) /
			CLOCKS_PER_SEC;
		long double pd = 100.0 * fabsl(experimental - trial) /
			(0.5 *(experimental + trial));
		printf("%16.15lf\t%16.15le\t%6d\t%4.3lf\n", 
			trial, pd, steps, runtime);
	}
	return(0);
}