Using the MATLAB Meshgrid Command and Array Operators to Implement One-and Two-Variable Functions

One of the more difficult topics covered in our MATLAB Fundamentals course is the use of array operators on matrices to implement functions of one and two variables. Although this is a powerful feature of MATLAB, the concept can be difficult to understand, since it does not exist in other programming languages.

One-Variable Functions

Let's assume that you have a function that depends on a single variable x,

y=f(x)=x2

You can define and visualize the function over a certain range in MATLAB, as follows:

1. First, create a vector with the desired resolution:

>> x=0:0.1:10;

2. Next, define the function on the vector, using the array operator '.^':

>> y=x.^2;

Since each element in the vector y is a function of the corresponding element in the vector x, you must use the array operator '.^' to act on each element separately. If you use the matrix operator '^' instead, it will act on the whole vector and fail.

3. Finally, visualize the function in the x and y dimensions using the plot command:

>> plot(x,y)

Two-Variable Functions

Now let's assume that you have a function that depends on two variables, x1 and x2:

y=f(x1,x2)=x12+x22

In the previous example, the single free variable is represented by a single axis on a 2-D plot. In this example, the two free variables must be represented as two axes on a 3-D plot. The two axes, when viewed from above, appear as follows:

The function is defined on the plane that is spanned by these two axes.
You can define and visualize this function as follows:

1. First, create vectors with the desired resolution:

>> x1=0:4; x2=0:5;

2. Next, define the function on these two variables. The first instinct is to apply array operators to the vectors:

>> y=x1.^2+x2.^2;

However, this method does not work, for two reasons: First, x1 and x2 are of different sizes and cannot be added together. Second, even if the vectors were the same size, the above command defines the function on the two axes, but not on the whole plane.
To define the function on the whole plane, you must define grid points on the plane as follows:

This grid is a collection of coordinate pairs, arranged in 2-D. Now, let's put the first coordinates in a matrix, x1g:

  x1g=  0    1    2    3    4
        0    1    2    3    4
        0    1    2    3    4
        0    1    2    3    4
        0    1    2    3    4
        0    1    2    3    4

and the second coordinates in another matrix, x2g:

  x2g=  0    0    0    0    0
        1    1    1    1    1
        2    2    2    2    2
        3    3    3    3    3
        4    4    4    4    4
        5    5    5    5    5

Now, if we use array operators on x1g and x2g, we are operating on first and second coordinates of each point on the grid, and that's exactly what we want. Thus,

>> y=x1g.^2+x2g.^2;

You can create x1g and x2g using the MATLAB meshgrid command. Given two free variables, x1 and x2, you can create the first and second coordinate matrices of a grid by entering:

>> [x1g,x2g]=meshgrid(x1,x2);

To summarize, here are the required steps to define and visualize a function of two variables.
First, define the free variables:

>> x1=0:4;x2=0:5;

Second, create the grid:

>> [x1g,x2g]=meshgrid(x1,x2);

Third, define the function on the grid using array operators:

>> y=x1g.^2+x2g.^2;

Finally, use a 3-D plotting command to visualize your data. For example:

>> surf(x1g,x2g,y)

Vectorization in MATLAB

Now let's consider why we do math in this manner in MATLAB. MATLAB is an interpreted language. When you type a command, the interpreter parses the command (checking syntax) and sends instructions to the processor to carry out the operations you have described with your command. This is different from compiled languages, like FORTRAN and C, where your code is first converted into instructions, and then these instructions are run to generate your results.

Every command in standard M-code requires that the MATLAB interpreter parse the code, checking the syntax and creating the instructions that need to be run. This process is generally so fast that you don't notice it. However for looping constructs (FOR and WHILE loops), the interpretation process consumes valuable time, and can produce a noticeable difference in the speed at which your code runs.

Consider the following two examples of array operations:

>> x = 0:0.1:10;
>> y = x.^2;

In C and FORTRAN, you must write a loop to perform this operation. An equivalent M-code, non-vectorized computation is:

>> for i = 1:length(x)
y(i) = x(i)^2;
end

The first example executes faster than the second, and is the way MATLAB is meant to be used.

You can try another example with larger vectors, to see the difference:

>> x = -10:.001:10;
>> tic; y = x.^2; toc;
elapsed_time =
0.1000
>> clear
>> x = -10:.001:10;
>> tic; for i = 1:length(x)
y(i) = x(i)^2;
end; toc;
elapsed_time =
21.8350

This example clearly demonstrates the benefits of vectorization in MATLAB. Looping constructs can usually be rewritten using array operations to generate faster, simpler, vectorized code.

Even in cases where your code cannot be vectorized, you may be able to use the following tips to improve code speed:

Pre-allocation

Dynamically sizing vectors requires system time to find a new piece of memory to hold the current matrix plus the elements being added. For example:

>> len = 1000;
>> for i = 1:len
x(i) = rand(1);
end

Variables in MATLAB require a contiguous space in memory for all elements of the array. Growing a vector by just one element requires MATLAB to find a new location in memory large enough for the entire array, copy the array to the new location, and add the new element(s) to the end of the array. Pre-allocating for variables whose size you already know can save on time spent accessing memory. The above example can be rewritten as:

>> len = 1000;
>> x = zeros(len,1);
>> for i = 1:len
x(i) = rand(1);
end

The ZEROS command in the above example allocates a single block in memory for the array. Then the variable is indexed and changed; it is not increased in size.

Reverse Indexing

Another way to improve speed is to reverse-index your looping:

>> len = 1000;
>> for i = len:-1:1
x(i) = rand(1);
end

In this example, the first run through the loop preallocates the memory for x.