<Previous Page: Numerical Differentiation|Next Page: Numerical Solutions to First-Order ODEs>
In this page, we will learn about numerical integration, evaluating integrals numerically using the C programming language.
For users of other programming languages, please refer to the following pages:
Numerical Integration with Python
Numerical Integration with Julia
Numerical Integration with Fortran
To evaluate integrals numerically, let's first review integration. Here, we consider the integration of a function \(f(x)\) as follows: \[ S = \int^b_a dx f(x). \tag{1} \] The value of the integral \(S\) in equation (1) is known to equal the area of the region shaded with red diagonal lines in Figure (a). Based on this fact, the integral in equation (1) can be approximated as the sum of the areas of multiple trapezoids. For instance, the area shaded with red diagonal lines in Figure (a) can be approximated as the sum of the areas of three trapezoids shaded with blue diagonal lines in Figure (b). Assuming the integration interval \((a \le x \le b)\) is divided into three equal-width regions, the approximate value of the integral using the three trapezoids in Figure (b) can be calculated as follows: \[ h = \frac{b-a}{3}, \tag{2} \] \[ S_3 = \frac{f(a)+f(a+h)}{2}h + \frac{f(a+h)+f(a+2h)}{2}h + \frac{f(a+2h)+f(b)}{2}h. \tag{3} \] Here, \(h\) represents the width of each divided interval, and \(S_3\) represents the area of the region shaded with blue diagonal lines in Figure (b). This equation can be further rewritten as: \[ S_3 = h \times \left [ \frac{f(a)}{2} + f(a+h) + f(a+2h) + \frac{f(b)}{2} \right ]. \tag{4} \] This expression is useful for constructing approximation formulas for integration using more trapezoids, as we will see next.
In the above discussion, we considered evaluating the integral in equation (\(1\)) using the approximation formula in equation (\(4\)) with three trapezoids. As expected from Figure (b), increasing the number of divided intervals and using more trapezoids to approximate the area allows for a more accurate approximation of the integral. For example, if the integration interval \((a \le x \le b)\) is divided into \(N\) subintervals, and the integral value for each subinterval is approximated using the area of a trapezoid, the approximate value of the integral as the sum of \(N\) trapezoidal areas can be expressed as: \[ h = \frac{b-a}{N}, \tag{5} \] \[ x_n = a + n \times h ~~~\rm{for}~~~ (0 \le n \le N), \tag{6} \] \[ S_n = h \times \left [\frac{f(x_0)}{2} + f(x_1) + \cdots + f(x_{N-1}) + \frac{f(x_N)}{2} \right ]. \tag{7} \] Here, \(h\) represents the width of each subinterval, \(x_n\) represents the coordinate of the \(n\)-th grid point, and \(n\) is an integer between \(0\) and \(N\). In this equation (\(7\)), the value of the integral in equation (\(1\)) is approximated using the areas of \(N\) trapezoids.
The approximation formula for integration in equation (\(7\)) is known as the trapezoidal rule. In the limit where the number of divisions \(N\) becomes large (\(N \rightarrow \infty\)), the formula gives the exact integral value.
In this section, we numerically evaluate the following integral using the trapezoidal rule (Equation (\(7\))): \[ \int_0^1 dx \frac{4}{1+x^2}=\pi \tag{8}. \] The code below divides the integration range (\(0 \leq x \leq 1\)) into \(N=20\) trapezoids and evaluates the integral. Additionally, an explanation of the code is provided below the code for your reference.
#include <stdio.h>
#include <math.h>
int main() {
int i, n;
double a, b;
double x, h, f, s;
n = 20;
a = 0.0;
b = 1.0;
h = (b - a) / n;
x = a;
f = 4.0 / (1.0 + x * x);
s = 0.5 * h * f;
for (i = 1; i < n; i++) {
x = a + h * i;
f = 4.0 / (1.0 + x * x);
s = s + h * f;
}
x = b;
f = 4.0 / (1.0 + x * x);
s = s + 0.5 * h * f;
printf("Result:\n");
printf("%26.16e\n", s);
return 0;
}
The above C code performs numerical integration using the trapezoidal rule. Here is an explanation of the code:
f
is defined to calculate the value of the function \( \frac{4.0}{1.0+x^2} \) for a given numeric value x
.trapezoidal_rule
is defined, which takes the integration range \([a, b]\) and the number of subdivisions n
as inputs and evaluates the integral using the trapezoidal rule.trapezoidal_rule
function, the width of each trapezoid is calculated as h
. The function computes the average value at the endpoints and then adds the contributions from the interior points to the sum. Finally, the sum is multiplied by h
to obtain the approximate value of the integral.main
function, the integration range \([a, b]\) and the number of subdivisions n
are set, and the trapezoidal_rule
function is called to display the result.Does the result approximate \(\pi = 3.14159...\)? In the above code, increasing n
will yield more accurate results. Try changing the value of n
to observe how the error changes.
<Previous page: Numerical Differentiation|Next page: Numerical Solutions of First-Order Ordinary Differential Equations>