Complex Numbers and Plotting in Matlab
When plotting in Matlab, whether it be in two or three dimensions, a number of issues involving complex numbers will arise.
If you take the square root of a negative number, the result is a complex number.
>> sqrt(-9)
ans =
0 + 3.0000i
That is, `sqrt(-9)=0+3i`. In general, complex numbers have the form `a+bi`, where `a` and `b` are arbitrary real numbers and `i=sqrt(-1)`.
Two points regarding complex numbers bear particular importance:
- The number `a` is called the real part of the complex number `a+bi`.
- If a number `r` is real, then the real part of `r` is itself.
Both of these facts are extremely useful when plotting in Matlab. Let's demonstrate with an example.
Sketch the graph of `y=sqrt(4-x^2)` on the interval `[-3,3]`.
Solution: The reader should note that the domain of the function `y=sqrt(4-x^2)` is `[-2,2]`, so on the interval `[-3,-2) U (2,3]`, we will be taking the square root of a negative number. Normally, one would pay attention to the domain of the function and not attempt to sketch the plot at points outside the domain, but it's instructive to plunge ahead anyway and see how Matlab handles the situation.
x=linspace(-3,3); y=sqrt(4-x.^2); plot(x,y)
Let's add a grid, equalize the axes, then adjust the view.
grid on axis equal axis([-3,3,-3,3])
The above sequence of commands will produce the plot shown in Figure 1.
![A plot of y=sqrt(4-x^2) over [-3,3].](figure1.png)
A plot of `y=sqrt(4-x^2)` over the interval `[-3,3]`.
If you square both sides of `y=sqrt(4-x^2)`, you can then manipulate the result into the form `x^2+y^2=4`, a circle of radius 2 centered at the origin. Because a nonnegative square root is called for in the equation `y=sqrt(4-x^2)`, it is not surprising that the graph is the upper half of the circle, as seen in Figure 1. However, note the strange horizontal line segments at each end of the plot in Figure 1. The one on the left is the graph of `y=0` on the interval `[-3,-2)`, and the one on the right is the graph of `y=0` on the interval `(2,3]`. How can these segments be explained?
The answer to our question lies in how Matlab's plot command treats data with complex entries. First, use Matlab's isreal command to check if all of the entries in the variable y are real numbers.
>> isreal(y)
ans =
0
If all of the entries of y were real, Matlab's isreal command would have responded with an answer of 1 (true). The zero response implies that the variable y has come entries that are complex numbers. You can see this for yourself by entering the variable y at the Matlab command prompt and hitting Enter.
>> y
y =
Columns 1 through 3
0 + 2.2361i 0 + 2.1541i 0 + 2.0706i
Columns 4 through 6
0 + 1.9855i 0 + 1.8985i 0 + 1.8093i
Columns 7 through 9
0 + 1.7177i 0 + 1.6231i 0 + 1.5251i
Columns 10 through 12
0 + 1.4230i 0 + 1.3157i 0 + 1.2019i
Columns 13 through 15
0 + 1.0795i 0 + 0.9452i 0 + 0.7931i
Columns 16 through 18
0 + 0.6098i 0 + 0.3495i 0.3468
Columns 19 through 21
0.5961 0.7636 0.8964
Columns 22 through 24
1.0082 1.1055 1.1919
Columns 25 through 27
1.2695 1.3399 1.4041
Columns 28 through 30
1.4630 1.5173 1.5673
Columns 31 through 33
1.6135 1.6562 1.6956
Columns 34 through 36
1.7321 1.7657 1.7966
Columns 37 through 39
1.8250 1.8510 1.8746
Columns 40 through 42
1.8961 1.9153 1.9325
Columns 43 through 45
1.9477 1.9608 1.9720
Columns 46 through 48
1.9813 1.9887 1.9943
Columns 49 through 51
1.9979 1.9998 1.9998
Columns 52 through 54
1.9979 1.9943 1.9887
Columns 55 through 57
1.9813 1.9720 1.9608
Columns 58 through 60
1.9477 1.9325 1.9153
Columns 61 through 63
1.8961 1.8746 1.8510
Columns 64 through 66
1.8250 1.7966 1.7657
Columns 67 through 69
1.7321 1.6956 1.6562
Columns 70 through 72
1.6135 1.5673 1.5173
Columns 73 through 75
1.4630 1.4041 1.3399
Columns 76 through 78
1.2695 1.1919 1.1055
Columns 79 through 81
1.0082 0.8964 0.7636
Columns 82 through 84
0.5961 0.3468 0 + 0.3495i
Columns 85 through 87
0 + 0.6098i 0 + 0.7931i 0 + 0.9452i
Columns 88 through 90
0 + 1.0795i 0 + 1.2019i 0 + 1.3157i
Columns 91 through 93
0 + 1.4230i 0 + 1.5251i 0 + 1.6231i
Columns 94 through 96
0 + 1.7177i 0 + 1.8093i 0 + 1.8985i
Columns 97 through 99
0 + 1.9855i 0 + 2.0706i 0 + 2.1541i
Column 100
0 + 2.2361i
Note that there are complex entries at the beginning and end of the vector y. Note especially that the real part of each of these complex numbers is zero. This can be emphasized with the following command.
>> real(y)
ans =
Columns 1 through 7
0 0 0 0 0 0 0
Columns 8 through 14
0 0 0 0 0 0 0
Columns 15 through 21
0 0 0 0.3468 0.5961 0.7636 0.8964
Columns 22 through 28
1.0082 1.1055 1.1919 1.2695 1.3399 1.4041 1.4630
Columns 29 through 35
1.5173 1.5673 1.6135 1.6562 1.6956 1.7321 1.7657
Columns 36 through 42
1.7966 1.8250 1.8510 1.8746 1.8961 1.9153 1.9325
Columns 43 through 49
1.9477 1.9608 1.9720 1.9813 1.9887 1.9943 1.9979
Columns 50 through 56
1.9998 1.9998 1.9979 1.9943 1.9887 1.9813 1.9720
Columns 57 through 63
1.9608 1.9477 1.9325 1.9153 1.8961 1.8746 1.8510
Columns 64 through 70
1.8250 1.7966 1.7657 1.7321 1.6956 1.6562 1.6135
Columns 71 through 77
1.5673 1.5173 1.4630 1.4041 1.3399 1.2695 1.1919
Columns 78 through 84
1.1055 1.0082 0.8964 0.7636 0.5961 0.3468 0
Columns 85 through 91
0 0 0 0 0 0 0
Columns 92 through 98
0 0 0 0 0 0 0
Columns 99 through 100
0 0
This output explains the strange horizontal segments at each end of the plot in Figure 1. Matlab's plot command ignores the imaginary part of any complex data and plots only the real part, which in the case of this example, is always zero. Hence the horizontal segments `y=0` at each end of the plot in Figure 1.
We can replace each complex entry in the vector y with a NaN. First, use Matlab's find command to find all entries that are not real. If a number is real, then its real part equals itself. Thus, we want to find the entries in the vector y whose real part does not equal itself.
>> k=find(real(y)~=y)
k =
Columns 1 through 12
1 2 3 4 5 6 7 8 9 10 11 12
Columns 13 through 24
13 14 15 16 17 84 85 86 87 88 89 90
Columns 25 through 34
91 92 93 94 95 96 97 98 99 100
Note that k now contains the indices of the entries in vector y that are complex. We will now replace the entries in y with these indices with NaN's.
>> y(k)=NaN
y =
Columns 1 through 7
NaN NaN NaN NaN NaN NaN NaN
Columns 8 through 14
NaN NaN NaN NaN NaN NaN NaN
Columns 15 through 21
NaN NaN NaN 0.3468 0.5961 0.7636 0.8964
Columns 22 through 28
1.0082 1.1055 1.1919 1.2695 1.3399 1.4041 1.4630
Columns 29 through 35
1.5173 1.5673 1.6135 1.6562 1.6956 1.7321 1.7657
Columns 36 through 42
1.7966 1.8250 1.8510 1.8746 1.8961 1.9153 1.9325
Columns 43 through 49
1.9477 1.9608 1.9720 1.9813 1.9887 1.9943 1.9979
Columns 50 through 56
1.9998 1.9998 1.9979 1.9943 1.9887 1.9813 1.9720
Columns 57 through 63
1.9608 1.9477 1.9325 1.9153 1.8961 1.8746 1.8510
Columns 64 through 70
1.8250 1.7966 1.7657 1.7321 1.6956 1.6562 1.6135
Columns 71 through 77
1.5673 1.5173 1.4630 1.4041 1.3399 1.2695 1.1919
Columns 78 through 84
1.1055 1.0082 0.8964 0.7636 0.5961 0.3468 NaN
Columns 85 through 91
NaN NaN NaN NaN NaN NaN NaN
Columns 92 through 98
NaN NaN NaN NaN NaN NaN NaN
Columns 99 through 100
NaN NaN
The beauty of this approach is the fact that Matlab's plot command ignores NaN's when plotting.
plot(x,y) grid on axis equal axis([-3,3,-3,3])
The above sequence of commands produces the plot in Figure 2.
The plot with NaN's replacing complex entries.
Note that the horizontal line segments are gone at each end because the plot command ignores the data in y containing NaN's. However, now we've introduced a new problem. Our linspace command is not fine enough to accurate portray the semicircle near its endpoints. Readers are encouraged to repeat the experiment above, increasing the number of entries in x. For example, try x=linspace(-3,3,500) or x=linspace(-3,3,1000) to see what differences occur.
Of course, some advance preparation would have eliminated the need for this discussion. By simply limiting the plot to the domain of the function, namely `[-2,2]`, we could have secured data with no complex entries. However, the discussion of complex number and NaN's will prove invaluable when plotting surfaces that introduce complex entries. Let's look at an example.
Sketch the graph of `z=sqrt(4-x^2-y^2)` on the domain `{(x,y): -3 le x,y le 3}`.
Solution: This problem bears a striking similarity to the function of Example 1, so once again, expect difficulties with complex numbers. We begin by initailizing the grid. (Note: Readers who have not yet worked with surfaces should consider reading Planes in Matlab or Cylinders in Matlab before attempting the following.)
x=linspace(-3,3,40); y=linspace(-3,3,40); [x,y]=meshgrid(x,y);
Now compute `z`.
z=sqrt(4-x.^2-y.^2);
Now, here lies the difficulty. To find the domain of `z=sqrt(4-x^2-y^2)`, set
`4-x^2-y^2 ge 0`,
then rearrange.
`x^2+y^2 le 4`,
Hence, the domain of `z=sqrt(4-x^2-y^2)` is `{(x,y): x^2+y^2 le 4}`. The domain is all of the pairs `(x,y)` that line on or in the interior of `x^2+y^2=4`, a circle of radius 2 centered at the origin. However, we are computing `z=sqrt(4-x^2-y^2)` on the square grid `{(x,y): -3 le x, y le 3}`, many of whose points lie outside the boundary of the circle `x^2+y^2=4`. Thus, our vector z will contain complex entries.
>> isreal(z)
ans =
0
The false response (0 indicates 'false' in Matlab) is an indication that not all of the entries of z are real; i.e., some entries are complex.
You can view the entries of z, much as we did the entries in y in Example 1, by typing z at the command prompt, followed by the Enter key. If you scroll through the entries in z you'll find the complex entries.
How Matlab's mesh command will perform depends on your version of Matlab. During the development stages of version 6, Matlab's mesh command would fail if it detected complex numbers in the data. Then, in the early development of Matlab 7, an adjustment was made to Matlab's mesh command to allow it to work with complex data; specifically, an adjustment in the code allowed it to deal with complex data in the same way as the the plot command dealt with complex data: it plotted the real part of any complex data entry. Users must have complained, because later versions of Matlab 7 have returned to the days of Matlab 6 and the mesh command will not work if entries have data that is complex.
By taking the real part of `z`, complex numbers are avoided. The following command should work on all systems.
mesh(x,y,real(z))
grid on
box on
axis equal
xlabel('x-axis')
ylabel('y-axis')
zlabel('z-axis')
title('The graph of z = sqrt(4 - x^2 - y^2) in three-space.')
The result of the above sequnce of commands is the surface shown in Figure 3.
Plotting the real values of `z`.
If you square both sides of `z=sqrt(4-x^2-y^2)`, then it's simple to produce the result `x^2+y^2+z^2=4`, a sphere of radius 2 centered at the origin. Because `z=sqrt(4-x^2-y^2)` calls for a nonnegative square root, it is not surprising that we get the upper hemisphere shown in Figure 3.
As you can see, all entries outside the domain `{(x,y): x^2+y^2 le 4}` have their `z`-value set to zero (because the real parts of these complex entries all equal zero). However, these are points that should not be plotted. Again, we can eliminate these points as we did in Example 1, by finding the complex entries of z and setting those entries to NaN's.
k=find(real(z)~=z); z(k)=NaN;
If you leave the semicolon off the end of the find command, you will be able to examine the entries in the vector k, which are indices into entries in z that contain complex numbers.
Now, because z has no complex entries, we no longer have to use real(z) in the mesh command. The mesh command will ignore data containing NaN's. Similar comments are in order for the vector z; i.e., remove the semicolon at the end of the command that assigns NaN's to those values of z with indices from the vector k.
mesh(x,y,z)
Now repeat the grid, turn "on" the box, equalize the axes, and add the annotations.
grid on
box on
axis equal
xlabel('x-axis')
ylabel('y-axis')
zlabel('z-axis')
title('The graph of z = sqrt(4 - x^2 - y^2) in three-space.')
The result of the above sequence of commands is the surface shown in Figure 4.
Replacing complex entries in z witn NaN's.
As we saw in Example 1, we have some erratic behavior on the boundary (`x^2+y^2=4`) of the hemisphere. You can try to eliminate this behavior by increasing the number of points in the mesh, but in upcoming activities we'll demonstrate how a change of coordinates is a better solution.
Matlab Files
Although the following file features advanced use of Matlab, we include it here for those interested in discovering how we generated the images for this activity. You can download the Matlab file at the following link. Download the file to a directory or folder on your system.
The file complex.m is designed to be run in "cell mode." Open the file complex.m in the Matlab editor, then enable cell mode from the Cell Menu. After that, use the entries on the Cell Menu or the icons on the toolbar to execute the code in the cells provided in the file. There are options for executing both single and multiple cells. After executing a cell, examine the contents of your folder and note that a PNG file was generated by executing the cell.
Exercises
- Consider the function `y=sqrt(9-x^2)`.
- Use Matlab's plot command (with no adjustments) to draw the graph of `y=sqrt(9-x^2)` on the interval `[-4,4]`. Use writng to explain any mysterious horizontal segments in your plot.
- Find and replace all complex entries with NaN's, then replot the function `y=sqrt(9-x^2)`.
- Submit plots for parts (a) and (b) above and the M-code that produced them.
- Consider the function `z=sqrt(1-x^2/4-y^2/9)`.
- Use Matlab's mesh command to draw the graph of `z=sqrt(1-x^2//4-y^2//9)` on the domain `{(x,y): -3 le x le 3, -4 le y le 4}`. Use writng to explain any mysterious horizontal planes in your plot.
- Find and replace all complex entries with NaN's, then replot the function `z=sqrt(1-x^2//4-y^2//9)`.
- Submit plots for parts (a) and (b) above and the M-code that produced them.