[Math] Trapezoidal integration in matlab – code explanation

MATLAB

I have a simple matlab script for calculating an Integral using trapezoidal rule.

The idea is calculate Fourier series, and it used for the Fourier coefficients.

function [area] = recArea(func, xl, xr, N, n)
    if nargin < 5
        n = 1;
    end
    w = ((xl - xr) / N);
    y = eval(func);
    area = sum((y(1:end-1) + y(2:end)) / 2*w);
end

Call example:

$a_n$ = recArea(func, xl, xr, N)

in the function:

function [y] = fourier(m, N, func, xl, xr)

Where:

$func$ is the function to calculate it's fourier sereis

$xl, xr$ the integral bounds: $\int_{xl}^{ar}$.

Can you please explain me how this code works?

Thanks in advance.

Best Answer

The domain $[x_l, x_r]$ is uniformly partitioned with $N$ sub-intervals, meaning that the width of any sub-interval is $\frac{x_r-x_l}{N}$.

For trapezoid rule integration, you compute the area of the trapezoid formed by the line segment connecting $(x_i,f(x_i))$ to $(x_{i+1},f(x_{i+1}))$ (draw it). The area under this line is basically the area of the rectangle formed by drawing a line segment parallel to the x-axis from $(x_i,f(x_i))$ to $(x_{i+1},f(x_i))$, which can be expressed by:

$$A_1 = \textrm{width}\times\textrm{height} = \frac{x_r-x_l}{N}\times f(x_i).$$

Then, you add the area of the triangle whose base is the line segment just described, and whose hypotenuse is the line from $(x_i,f(x_i))$ to $(x_{i+1},f(x_{i+1}))$. This has area:

$$A_2 = \frac12\textrm{base}\times\textrm{height} = \frac12 \frac{x_r-x_l}{N} \times \left(f(x_{i+1})-f(x_i)\right).$$

Add $A_1$ and $A_2$, and you get the area of a sub-interval as:

$$A_i = A_1+A_2 = \frac{x_r-x_l}{2N}\left(f(x_i)+f(x_{i+1})\right).$$

Sum these over all intervals, and that exactly describes what's happening by the line of code that says area = sum(....

Related Question