I recently had to implement multi-application trapezoidal rule for some work I was doing. I decided to implement the algorithm in MATLAB, and figured I’d post a tutorial for it, since I know many people struggle with implementing algorithms in code.
This implementation takes advantage of MATLAB’s built in functions, and as such, requires no loop. This technique is referred to as ‘vectorization’ and in many cases, can speed up your code significantly, as it lets MATLAB take advantage of its built in matrix features, which rely on LAPACK and LINPACK - matrix math libraries that have been hand-optimized for various processor architectures. Pretty cool!
We’ll build this up in a script so it’s easy to see all of the variables. I often play around with scripts as I’m developing in MATLAB, and encapsulate the code in a function later.
Step 1: Define Function
Trapezoidal rule is a method of numerical integration. It’s not a great one, but it is easy to understand. I’ll assume the reader already has an understanding of how to do it.
Let’s use an anonymous function to contain the specific function we’d like to integrate. This lets us use standard math notation F(x) to fun the function in our code.
F = @(x) x^3 + 2*x +5;
Step 2: Setup Integration Variables
Now we need some variables to contain the start and end points of the integral, and the number of sections we’d like to use in our trap rule.
start = 0; endpt = 5; n = 3;
We’ll also generate the points to evaluate the function at using the
linspace function, which generates a set of n evenly spaced points between start and endpt.
Also, we’ll generate
h, our step width using the points we just calculated.
points = linspace(start, endpt, n+1); h = abs(points(2) - points(1));
Step 3: Implement Algorithm
The algorithm for multi application trapezoidal rule is the following:
Multiple application trapezoidal rule
We can implement this in MATLAB, using the
sum() function as well as the
arrayfun() function, which will apply each element of the points array to our anonymous function
F. One important note is that the algorithm above has the points zero-indexed. MATLAB arrays use 1-indexing, so we need to add 1 to each subscript of
x in the image above.
int = (h/2) * (F(points(1)) + 2*sum(arrayfun(F, points(2:n))) + F(points(n+1)))
Step 4: Vectorize
Right now, we have
arrayfun(), which is (likely) one step better than using a loop. Even better would be using full vectorization. To do this, we need to reformat the f function to take vectors. So, we’ll use the much-misunderstood dot operator
. ahead of our caret
^ to make the power element-wise. This means the power will be applied to each element of the input vector, which is what we want. It’s not needed in this application, but using the . before multiplication will result in a more robust function as well, so let’s do it to demonstrate.
f = @(x) x.^3 + 2*x +5;
The trap expression will now look like this:
int = (h/2) * (f(points(1)) + 2*sum(f(points(2:n))) + f(points(n+1)));
Let’s plot the original function at a higher resolution, and plot the trap points we evaluated. I’ll take advantage of linspace again, to plot our “base” function with 30 times more points than n. For small values of
n, this will clearly show the points we use in trap rule. For higher n values, it’ll all smush together.
% plot a high resolution curve of the function and superimpose the trap % points plot(linspace(start,endpt,n*30),f(linspace(start,endpt,n*30))) hold on plot(points, f(points)) title(sprintf('Trap Rule n = %i', n))
Notice that I used
sprintf to display
n in the title of the plot. Sprintf is used for formatting data into strings, which is what the
title() function takes.
The Entire Thing
Here’s the entire script.
%% Multi application trap rule % Richard Arthurs, March 20, 2017 % Richarthurs.com f = @(x) x.^3 + 2*x +5; % Function to integrate start = 0; endpt = 5; n = 3; % number of sections of trap rule points = linspace(start, endpt, n+1); h = abs(points(2) - points(1)); int = (h/2) * (f(points(1)) + 2*sum(f(points(2:n))) + f(points(n+1))) % the integral approximation % plot a high resolution curve of the function and superimpose the trap points plot(linspace(start,endpt,n*30),f(linspace(start,endpt,n*30))) hold on plot(points, f(points)) title(sprintf('Trap Rule n = %i', n))