Before delving into Julia programming, let's first explore the fundamental concept of numerical integration, which involves using numerical methods to evaluate the integral of a function. Consider the following one-dimensional integral: \[ S=\int_{a}^{b} dx \, f(x). \tag{1} \] Now, envision an \(N\)-division of the integral interval (\(a \le x \le b\)). The width of each division, \(\Delta x\), is defined as \(\Delta x = (b - a)/N\). We can express the integral \(S\) as the sum of integrals over each divided interval: \begin{align} S &= \int_{a}^{b} dx \, f(x) \\ &= \int_{a+\Delta x}^{a} dx \, f(x) + \int_{a+2\Delta x}^{a+\Delta x} dx \, f(x) + \int_{a+3\Delta x}^{a+2\Delta x} dx \, f(x) + \cdots + \int_{b}^{a+(N-1)\Delta x} dx \, f(x) \\ &= \sum_{j=0}^{N-1} \int_{a + j\Delta x}^{a + (j+1)\Delta x} dx \, f(x) \tag{2} \end{align}

Assuming that the division width \(\Delta x\) is sufficiently small, the integrand \(f(x)\) can be approximated as the average of the values at the ends of each small integral interval. In other words, Eq. (2) can be approximated as follows: \begin{align} S &= \sum_{j=0}^{N-1} \int_{a + j\Delta x}^{a + (j+1)\Delta x} f(x) \, dx \nonumber \\ &\approx \sum_{j=0}^{N-1} \int_{a + j\Delta x}^{a + (j+1)\Delta x} \frac{f \left ( a + (j+1)\Delta x \right )+f\left (a + j\Delta x \right )}{2} \, dx \nonumber \\ &= \sum_{j=0}^{N-1}\frac{f \left ( a + (j+1)\Delta x \right )+f\left (a + j\Delta x \right )}{2} \int_{a + j\Delta x}^{a + (j+1)\Delta x} \, dx \nonumber \\ &= \sum_{j=0}^{N-1} \frac{f \left ( a + (j+1)\Delta x \right )+f\left (a + j\Delta x \right )}{2} \, \Delta x \nonumber \\ &= \frac{f(a)}{2}\Delta x + \left [ \sum_{j=1}^{N-1} f(a + j\Delta x) \, \Delta x \right ] + \frac{f(b)}{2}\Delta x \tag{3} \end{align}

The formula presented in Eq. (3) is commonly known as the Trapezoidal Rule, named so because it approximates the integral value by summing the areas of trapezoids. In the limit where the number of divisions is large (\(N\rightarrow \infty\), \(\Delta x \rightarrow 0\)), the trapezoidal formula in Eq. (3) converges to the exact integral value \(S\).

As an exercise, let's compute the following integral using the trapezoidal formula: \begin{align} S = \int_{0}^{1} dx \, \frac{4}{1+x^2} = \pi . \end{align} Below is a Julia code that numerically evaluates the aforementioned integral. To gain familiarity with numerical integration and Julia programming, try creating your own code based on the provided sample code.

# define your function
function my_func(x)
    return 4.0 / (1.0 + x^2)
end

# integrate the function from a to b with n points
function integrate(f, a, b, n)
    h = (b - a) / n

    s = (f(a) + f(b)) / 2.0
    for i in 1:n-1
        x = a + i * h
        s += f(x)
    end

    s *= h
    return s
end

a = 0.0
b = 1.0
n = 64

num_integral = integrate(my_func, a, b, n)

println("num. integral = ", num_integral)
println("pi            = ", 1.0*π)

The provided Julia code defines two functions and performs a numerical approximation of the integral for a specific mathematical function.

1. The first function, my_func(x), takes a numerical input x and calculates the value of \(\frac{4.0}{1.0+x^2} \).
2. The second function, integrate(f, a, b, n), computes the numerical integral of a given function f over the interval \([a, b]\) using the trapezoidal rule with n points. It initializes the variable h as the width of each trapezoid, computes the average height of the function at the interval endpoints, iterates through the interior points, adding their contributions to the total area. Finally, it multiplies the result by the width of each trapezoid to obtain the approximate integral value.

After defining these functions, the code sets specific values for the interval \([a, b]\) (0.0 to 1.0) and the number of points n (64). Subsequently, it calls the integrate function with the function my_func and the specified interval and number of points, storing the result in the variable num_integral.

Finally, the code prints the numerical integral result (num_integral) and compares it to the value of π, displaying both values.