1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#include <math.h>
#include <stdio.h>

double
f(double x)
{
    return x * x;
}

double
F(double x)
{
    return x * x * x / 3;
}


// some simple rules for numerical integration

double
midpoint(double a, double b, size_t n, double f(double))
{
    double h = (b - a) / n;

    double res = 0;
    for (size_t i = 0; i < n; ++i) {
        res += f(a + (i + 0.5) * h);
    }
    return res * h;
}

double
trapezoidal(double a, double b, size_t n, double f(double))
{
    double h = (b - a) / n;

    double res = f(a) + f(b);
    for (size_t i = 1; i < n; ++i) {
        res += 2 * f(a + i * h);
    }
    return res * h / 2;
}

// apply an integration rule to a given function

void
applyRule(const char *ruleName, double a, double b, size_t n, double f(double),
          double rule(double, double, size_t, double(double)))
{
    printf("%s(%lf, %lf, %zu) = %lf\n", ruleName, a, b, n, rule(a, b, n, f));
}

// print a report for comparing numerical results

void
compare(double a, double b, size_t n, double f(double))
{
    printf("compare:\n");
    applyRule("  midpoint", a, b, n, f, midpoint);
    applyRule("  trapezoidal", a, b, n, f, trapezoidal);
}

int
main(void)
{
    printf("f(x) = x^2\n");
    compare(1, 2, 5, f);
    compare(1, 2, 10, f);
    printf("exact result: %lf\n\n", F(2) - F(1));

    printf("f(x) = sin(x)\n");
    compare(1, 2, 5, sin);
    compare(1, 2, 10, sin);
    printf("exact result: %lf\n\n", -cos(2) + cos(1));
}