Integration and Differentiation in MATLAB: Video Tutorial

This video tutorial reviews and implements integration and differentiation in MATLAB. It is part of the Numerical Computations course. To watch the full course, visit the course page.

Two kinds of problems

The four videos cover two situations:

• You have a function. The integral $\int_a^b f(x)\,dx$ is computed by adaptive quadrature: the solver evaluates $f$ at chosen points and refines the intervals until the error is small enough. MATLAB offers integral and quadgk for one variable, integral2 and quad2d for double integrals, and integral3 for triple integrals.
• You have samples. When only values $y_i$ at points $x_i$ are known, sums replace integrals and differences replace derivatives. MATLAB offers sum, cumsum, trapz and cumtrapz for integration, and diff, gradient and del2 for differentiation.

Watch online

The four sections of this tutorial are available on YouTube and are embedded here as a playlist.

Play playlistLoads from YouTube

What each video covers

Section 1: Numerical integration using integral and quadgk

Watch this part on YouTube (opens in a new tab).

The basic call is integral(f, a, b), where f is a function handle. The video starts with $\int_1^2 x\,dx = 1.5$, then integrates $x e^{-x^2}$ from 0 to 1 (about $0.3161$) and from 0 to infinity (exactly $1/2$). The integrand uses element-wise operators (.*, ./, .^), because integral evaluates it on many points at once.

f = @(x) x.*exp(-x.^2);
q = integral(f, 0, Inf)

Next comes $\int_0^1 \ln x\,dx = -1$, which has a singularity at the lower limit. The video uses this example to explain the error tolerances. If $Q$ is the exact value and $q$ the computed one, the solver tries to satisfy

Setting 'RelTol' to 0 leaves only the absolute tolerance in control. A smaller 'AbsTol' requests a tighter error estimate, but it does not remove floating point limits or guarantee that every difficult integrand reaches that accuracy. Check any warning from the solver.

The video then computes contour integrals in the complex plane. The integral of $1/z$ around the unit circle is $2\pi i$ (Cauchy’s theorem says it would be zero if the path enclosed no pole). It is first computed with the parametrization $\gamma(t) = \cos t + i \sin t$:

A simpler way uses the 'Waypoints' option: the path starts and ends at 1 and passes through $i$, $-1$ and $-i$. A triangular path through $2 + i$ and $2 - i$ does not enclose the pole at the origin, and the result is almost zero.

Finally, quadgk (adaptive Gauss-Kronrod quadrature) is called in the same way and gives the same results. It accepts the same tolerance and waypoint options, plus 'MaxIntervalCount'.

Section 2: Double and triple integrals

Watch this part on YouTube (opens in a new tab).

integral2(f, xmin, xmax, ymin, ymax) and quad2d compute double integrals, and integral3 computes triple integrals. The limits of $x$ are numbers, while the limits of $y$ (and $z$) can be numbers or function handles of the outer variables. The video solves three examples:

In the second integral, the upper limit of $y$ depends on $x$, so it is passed as a function handle:

g = @(x, y) 1./(sqrt(x + y).*(1 + x + y).^2);
ymax = @(x) 1 - x;
q = integral2(g, 0, 1, 0, ymax)

Section 3: Discrete integrals and summations

Watch this part on YouTube (opens in a new tab).

In discrete systems, a sum plays the role of an integral. The video shows sum and cumsum on vectors and on a matrix along either dimension, and sum(Y(:)) for the total of all elements.

It then approximates $\int_0^{\pi/2} \sin x\,dx = 1$ with 100 samples. The plain sum of the samples times $\Delta x$ is close to 1. trapz uses the trapezoidal rule and gives a better result. With equal spacing you can pass the spacing $\Delta x$ or multiply by it. With unequal spacing you must pass the vector of $x$ values. cumtrapz returns the running integral, which the video plots against $1 - \cos x$. Nested calls of trapz give double and triple integrals.

x = linspace(0, pi/2, 100);
y = sin(x);
q = trapz(x, y)

Section 4: Numerical differentiation

Watch this part on YouTube (opens in a new tab).

diff is the opposite of cumsum: for a vector it returns the differences $x_{k+1} - x_k$, so the result has one element less. The derivative of $\sin x$ on 100 points from 0 to $2\pi$ is then approximated by the forward difference

and compared with the exact derivative, $\cos x$. Reducing the spacing improves the approximation in this smooth example until rounding error becomes important. Differencing noisy measurements can amplify the noise.

gradient returns a result of the same length as the input and gives a smaller error than diff in this example. It also works in several dimensions. For $f(x, y) = x\,e^{-(x^2 + y^2)}$ on a grid with spacing $h$, [gx, gy] = gradient(Z, h) returns both partial derivatives, and quiver draws the gradient vectors over a contour plot.

The last function is del2, which is related to the Laplacian

del2 returns $\nabla^2 f / (2N)$, where $N$ is the number of dimensions. For a function of two variables, the video multiplies the result by 4 to get the Laplacian.

Applications

• Areas, volumes and other quantities defined by integrals that have no closed form.
• Integrating measured data, such as velocity samples, to get a total or a running total.
• Estimating rates of change and gradients from sampled or grid data.

Check units and spacing

Pass the actual sample coordinates when measurements are not evenly spaced. A raw sum or difference does not include the sample spacing, so it does not by itself have the units of an integral or derivative. For a new integrand, compare tolerances and check against an analytic result when one is available. For sampled data, repeat the calculation on a finer grid when possible and check whether the result stabilizes.