A Mathematical Smörgåsbord

We have spent most of this book discussing the facilities MATLAB provides; now let’s take a look at a few hands-on examples. After all, scripts and plots and OOP are all very nice, but sometimes you just want to calculate something, and it is great to know in advance what function to call.

Elementary Math

No, really, I’m not going to reiterate primary school material. What I will talk about are some of the most basic facilities that are implemented via functions.

Let’s start with a simple fact: every operator, be it * or .^ or the matrix \ operator, has a correspondingly named function. For example, + can be invoked as

These notations are of little use to us (who would want to be using functions instead of operators?) except in cases where we want to use them as function handles, or when we want to implement (overload) it for our classes. For example, say you have a class called Size and you want to support two sizes being added together. In this case, you define your class like this:

The plus function takes care of the operator calls, and here’s how you would use it:

Next up, we have various cumulative calculation functions. For example, sum and prod calculate the sum or product of all the elements in an array (these functions generally do not care about array dimensions—they treat all arrays as flattened). In addition to cumsum (cumulative sum) that we’ve already seen, there is a corresponding cumprod (cumulative product) function.

There is also a function called diff that calculates the differences between adjacent elements in an array. One of its uses is to calculate the “delta” value of an array initialized with linspace:

Lastly, we have functions for calculating the floor and ceiling of a number—these are floor and ceil, respectively. Also, we have the round function that rounds to the nearest integer, and fix, which rounds towards zero:

Finally, mod does modulo division (MATLAB hijacked the % operator for comments, unfortunately) and rem calculates the remainder of division. These are identical for positive values but yield different results for negative ones:

Trigonometry, Exponents, and Logs

Life wouldn’t be fun without trigonometry, would it? As one would expect, MATLAB comes with support for all the common trig functions and their hyperbolic variations. Bizarrely, each function exists in two variations: an ordinary function that takes its arguments in radians, and the same function postfixed with –d (e.g., cosd) that takes degrees. Also, inverse functions are prefixed with a rather than arc, so the inverse cotangent would be acot.

The only function of note is hypot, which calculates the square root of the sum of squares (or the hypotenuse of a triangle).

Moving on, we have the fairly obvious exp (calculates ⅇx, note there is no “e” constant in MATLAB), log (the natural logarithm), and sqrt (x). log10 calculates the base 10 logarithm. There are also more obscure functions to calculate logarithms and exponents for small values of x as well as versions of exp/log/sqrt prefixed with real- that throw exceptions on negative input:

Complex Numbers

The support for complex numbers is built into MATLAB, and you get the i/j constants (functions, really, but you get the idea) for defining and manipulating complex values. In addition to ordinary arithmetic, you can calculate the real and imag (imaginary) parts of complex numbers, the complex conj (conjugate) and a few other, less exciting things.

Notice how the result of real is still a complex number.

Note that there is a difference between calling x' and transpose(x) on a matrix that has complex values: the first performs a conjugate transpose, the second a non-conjugate one:

Special Functions

While MATLAB doesn’t give you every single special function under the sun, it does give you the popular ones such as Beta, Gamma, a variety of Bessel functions, the Error function, and its complementary functions such as the Legendre function.

Also, MATLAB isn’t shy in using special functions in its own results when it cannot derive something. For example, using the Symbolic Toolbox to calculate ∫xn ex we get:

Of course, it would be very surprising if we actually got an analytic solution for the above integral!

Discrete Math

As we’ve mentioned in our discussion on basic syntax, MATLAB has hijacked the ! operator to send calls to the underlying operating system. This unfortunately prevents us from being able to use ! for factorials—instead MATLAB gives us the factorial function. Oh, and there is no double factorial or subfactorial support, which can be a problem if you want to solve that famous “2008” problem.

Discrete math functions let you calculate lots of independent things, like the greatest common divisor and least common multiple (gcd and lcm respectively), prime values such as the lowest prime less or equal to a value (primes), or all the prime factors of a number (factor).

Incidentally, you’ll sometimes see toolboxes reuse (overload) existing functions for their own, nefarious ends. The factor function is a good example:

Test Matrices

MATLAB has lots of different functions for generating various matrices useful for testing or demonstrations. We’ve already used the magic square matrix to show sorting, but there are lots of other matrices such as the Hilbert, Hadamard, Toeplitz, Hankel, and other matrices.

In addition to these, there is also a gallery function that has even more sample matrices (use help gallery to get a listing). For example, to generate a 5 × 5 matrix filled with zeros or ones, you can write:

Interpolation

MATLAB supports several types of interpolation. Let’s look at 1D interpolation first—this requires juts a single list of values. The spline function, for instance, lets us perform cubic interpolation, which lets us ”fill in the dots” on an incomplete data set:

And here’s the plot with x marks indicating the original values, and the line indicating the interpolated spline:

Original values and the interpolated spline.

In fact, the spline function is a specialization of the interp1 function with a ‘spline’ parameter. For 2D and 3D data we simply have interp2 and interp3 respectively. To illustrate 3D interpolation, let’s use the flow function to generate a coarse approximation of fluid flow and then interpolate over a finer mesh:

To render the image shown below, we used the slice function to effectively take and show slices out of a volumetric plot:

Slice plot of an interpolated flow model.

Incidentally, it’s worth mentioning that in addition to meshgrid, which works in 2D and 3D space, there is also a corresponding ngrid function that works in N-D space. Its operation is very similar to meshgrid: you simply provide the coordinate values for however many dimensions you have and get the corresponding set of arrays.

Optimization

MATLAB supports typical optimization problems: minimizing or maximizing single- or multivariate functions and solving constraint problems. Optimization is a very big topic, so we’ll just look at a few examples, and you can consult the MATLAB documentation and relevant texts on all the scary details of various optimization problems.

Consider the function fx=sinxex

A plot of the function fx=sinxex.

To locate the local minimum in the range (-4, 0) we need, first of all, to define the function to search over—MATLAB’s support for function literals comes in handy here. Then, we use the fminbnd function to find the minimum within the specified bounds:

There is no corresponding maximizing function: if you need the maximum of fx, simply calculate the mimimum of –fx.

In addition to supporting functions of one variable, MATLAB also supports searching for minimums of functions of many variables. For instance, let us suppose that you want to find the maximum of fx1,x2=x1ⅇ-x12-x22, which looks like this:

Surface plot of the function fx1,x2=x1ⅇ-x12-x22. A particularly good-looking plot that MATLAB documentation really likes to show off.

We can instruct MATLAB to search for a minimum starting somewhere in the (10, 20) range:

Note that, unlike the bounded fminbnd, fminsearch performs unconstrained optimization.

Let’s move on. The fzero function attempts to locate a root of one equation with one variable. For example, given the following function:

A plot of the function fx=x2+10x+16. This quadratic equation actually has two roots,
but we constrained the range in this chart, so you only see one.

We can locate the root by making the assumption that it lies somewhere in the (-6, 6) range:

We have deliberately turned on verbose output options for the above call so as to show you the iterative process. Without the extra setting, the fzero call just gives you the final result.

MATLAB supports a few more optimization algorithms, but a lot more fun can be had with the Optimization Toolbox (for local optimizations) and the Global Optimization Toolbox.

Numerical Integration and Differential Equations

First and foremost, if you want symbolic differentiation or integration (i.e., analytical solutions rather than numeric ones), you need the Symbolic Math Toolbox, because MATLAB doesn’t include this functionality by default. What MATLAB does give you is the ability to perform numerical integration and differentiation as well as methods for solving ordinary (ODE) and partial (PDE) differential equations.

Let’s start with numeric integration. This is done using the integral function, and requires the function itself as well as the integration limits:

Note: Can you see anything unusual about this code sample? That’s right—when working with numerical methods, your x value might be an array—ergo, we use the elementwise operator (x.^2). Had we used an ordinary exponentiation, we would get a rather cryptic error.

Double and triple integrals are evaluated with integral2 and integral3, respectively. MATLAB also comes with support for a variety of numerical integration methods such as Trapezoidal (trapz) and Simpson (quad).

Now, let’s move on to the fun stuff: solving ODEs (ordinary differential equations). MATLAB comes with several functions all prefixed with ode- followed by numbers and letters, such as ode15i and ode23tb. These are used to distinguish whether the equation to be solved is stiff or nonstiff, and what method (low order, medium order, etc.) is used to solve the equation.

For example, let’s solve the equation:

y'=-ty2-y2

with initial conditions y(0)=1,t∈[0,5]. For this, we first define our function:

We then use ode45 (nonstiff, medium order method) to solve the equation and plot the results:

A plot of the solution to our differential equation.

The PDE solver that comes with MATLAB supports 1D parabolic-elliptic PDEs with initial-boundary values. The function pdepe is used to solve the equation, and pdeval evaluates the numerical solution using pdepe’s output. Additional PDE support is provided by functions in the Partial Differential Equation Toolbox.

Fourier Analysis and Filtering

In the realm of Fourier analysis the key operation is, of course, the Fast Fourier Transform function fft as well as its 2D and N-D variants (fft2 and fftn) and their inverse transforms (ifft, ifft2 and ifftn respectively).

As is custom, we illustrate FFT by first creating some data that includes both oscillations as well as noise:

We now convert to the frequency domain and apply FFT to get the discrete Fourier transform of the noisy signal:

In the above, nextpow2 gets us the next power of 2 of a variable. Having performed the transform, we can now plot the spectrum:

One-sided spectrum plot of the function yt.

And there you have it—two clear peaks at 50 and 120. Of course, thanks to the noise, we don’t get the exact amplitudes, but then life isn’t perfect either.

So let’s now talk about filtering. MATLAB has functions for 1D and 2D digital filtering, 2D and N‑D convolutions, and a few other things besides. Filtering, which in the 1D case is done with the filter function, is rather interesting. Its signature takes three parameters: the first two are the numerator and denominator coefficient vectors, and last one is the data to filter.

The filter function is useful for lots of things. For example, let’s generate a random path and then plot the simple moving average of a few periods of data:

And here is the end result:

A Brownian motion path (blue) and its 30-period moving average (green).

Note: The author of this manuscript makes no claims whatsoever as to the effectiveness of technical analysis techniques (such as use of moving averages) for purposes of financial analysis

It’s virtually impossible to cover every feature of MATLAB in a single book, so we’ll stop here. It is a Succinctly title, after all.

The “2008” Problem

As promised, here’s the solution to the problem presented earlier. First, a reminder: you are given thirteen zeros (0 0 0 0 0 0 0 0 0 0 0 0 0) and you need to use these zeros with some mathematical operations to arrive at a value of 2008.

The trick here is to use the factorial function and the fact that 0!=1 to make some real numbers. The solution can then be acquired from adding up several 0! terms and taking their factorial so, for example

2008=3!+1!2-232

The use of double factorials allows for even shorter solutions. For example,

2008=3*-3!-13*+1-1

where 3*=3!‼=48. Then there are subfactorials and lots of other tricks to reduce the number of required zeros to a minimum. See if you can find some shorter solutions.