Composite Simpson 1/3 algorithm

Composite Simpson's 1/3 rule divides [a, b] into n (even) equal sub-intervals of width h = (b − a)/n and fits a parabola over each pair:

I ≈ (h/3) [ f₀ + 4(f₁ + f₃ + … + fₙ₋₁) + 2(f₂ + f₄ + … + fₙ₋₂) + fₙ ],   fᵢ = f(a + i·h)

Algorithm

1. Start.
2. Define the function f(x).
3. Read the limits a and b and the number of sub-intervals n.
4. If n is odd, add 1 to n (or ask again), because n must be even.
5. Compute h = (b − a)/n.
6. Set sum = f(a) + f(b).
7. For i = 1 to n−1:
   • x = a + i·h
   • if i is odd, sum = sum + 4·f(x); otherwise sum = sum + 2·f(x).
8. Compute I = (h/3)·sum.
9. Print I.
10. Stop.

C program

#include <stdio.h>
#include <math.h>

double f(double x) { return exp(x) + sin(2 * x); }   /* example function */

int main() {
    double a, b, h, x, sum, I;
    int n, i;
    printf("Enter a, b and n (even): ");
    scanf("%lf %lf %d", &a, &b, &n);
    if (n % 2 != 0) n++;                 /* n must be even */
    h = (b - a) / n;
    sum = f(a) + f(b);
    for (i = 1; i < n; i++) {
        x = a + i * h;
        if (i % 2 == 1) sum += 4 * f(x);  /* odd points */
        else            sum += 2 * f(x);  /* even points */
    }
    I = h / 3 * sum;
    printf("Integral = %f\n", I);
    return 0;
}

For a = 0, b = 2, n = 10, the program prints the integral of eˣ + sin 2x as about 7.2160 (exact 7.2159).