Chapter 3: MATLAB® graphics – Matlab® in Bioscience and Biotechnology

3

MATLAB® graphics

Plots are one of the best and most visually accessible means for presenting information in general and results of calculations in particular, especially in science and engineering. MATLAB® has a wide variety of commands for graphics. Available commands allow us to design various types of two- (in other words, XY or 2D) and three-dimensional (XYZ or 3D) plots.

2D graphics include plots with linear, semilogarithmic or logarithmic axes, bars, histograms, pies, polar and multiple plots, and many more. A plot can comprise one or several curves or surfaces, and several graphs can be plotted in the same Figure Window. Formatting can regulate the desired line style or marker type, its thickness and color, and add a grid, text or legend.

Plots having three axes are frequently used to represent data involving more than two variables. MATLAB® provides a variety of features for visualization of 3D data. These include spatial line plots, mesh and surface plots, geometrical figures and images. Generated plots can be formatted by commands or interactively from the Figure Window.

The most important commands for 2D and 3D plotting are described below. It is assumed that the reader has mastered the preceding chapter, and therefore lesser commands have inline explanations, and are written not within special frames (as in the second chapter), but as inline MATLAB® comments, next to the percentage sign (%).

3.1 Generation of XY plots

The basic command used for XY plotting is plot, which in its simplest form can be written as:

plot(x,y) or plot(y)

where x and y are two vectors of equal length, the first being used for horizontal and the second for vertical axes.

In the second plot-form the y values are plotted versus their indices.

After inputting the plot command with given values of x and y, the curve y(x) is created in the MATLAB® Figure Window with a linear scale by default. For example, we have mass (g/L) – time (min) data for an aerobic biomass process which are presented in Table 2.2. Let us plot a graph in which the x-axis has time values and the y-axis has biomass values. To create the plot we have to type in the Command Window:

> > t = [0:25:150 200];

> > b_mass = [5.15 5.21 5.5 6.55 7.15 7.75 7.59 7.45];

> > plot(t,b_mass)

Entering these commands in the Figure Window, the biomass–time plot resembles that shown in Figure 3.1.

Figure 3.1 Biomass data plotted in the Figure Window with default settings.

The line style and marker type, its thickness and color can be regulated using the plot command with additional optional arguments:

plot (x,y, 'Line Specifiers', 'Property Name', 'Prperty Value')

where Line Specifiers determines the line type, the marker symbol and the color of the plotted lines (see Table 3.1); Property Name assigns properties to the specified Property Values and can be taken from Table 3.2.

Table 3.1

Line style, color and marker type specifiers*

*Incomplete.

Table 3.2

Property names and property values*

Property name Purpose Property value
LineWidth or
linewidth
Specifies the width of the line A number in points (1 point = 1/72 inch). The default width is 0.5
MarkerSize or
markersize
Specifies the size of the marker (by the ‘.’
symbol)
A number in points
MarkerEdgeColor or markeredgecolor Specifies the color of the marker or the edge color for filled markers Color in accordance with specifiers in Table 3.1
MarkerFaceColor or markerfacecolor The fill color for markers that are closed shapes Color in accordance with specifiers in Table 3.1

*Incomplete.

Line specifiers, property names and property values are typed in the plot commands as strings in quotes (inverted commas). The specifiers and property names with their values can be written in any order, and one or more of them can be omitted. The omitted properties would be taken by default.

Some examples of application of specifiers and properties:

plot(y,'–c') plots the cyan solid line with x-equidistant y points.
plot(x,y,'g') plots the greed solid (default) line that connects the points.
plot(x,y,'--p') plots the blue (default) dashed line, points marked with the five-pointed asterisks.
plot(x,y,'k:h') plots the black dotted line, points marked with the six-pointed asterisks.
plot(t,b mass, … '–mo', 'linewidth',3, 'markersize',10, 'MarkerEdgeColor', … ' k ', 'MarkerFaceColor', 'y')  
  plots the magenta three-point solid line (1 point has size 1/72 inch), and x,y-values marked with the 10-point black-edged yellow circles.

Entering the last command with previously used biomass–time data, we get the plot shown in Figure 3.2.

Figure 3.2 Biomass data generated with specifiers and property settings in the plot command.

To close one Figure Window type and enter the close command in the Command Window; to close more than one Figure Window, use the close all command.

3.1.1 Two or more curves on a single 2D plot

To plot two or more curves in the same graph, at least two options can be used in MATLAB®: the first by typing the pairs of x,y-vectors into the plot command and the second by using the hold on … hold off commands.

3.1.1.1 The plot command option

The form of the plot command to create two or three curves in a single plot is:

plot(x1,y1,x2,y2) or plot(x1,y1,x2,y2,x3,y3)

These commands create graphs with two and three curves, respectively, where x1 and y1, x2 and y2, x3 and y3 are the pairs of equal-length vectors containing the x,y data. For more than three curves, new x,y-pairs should be added in the plot command. As an example, to plot two trigonometric functions, sine and cosine, in the same plot, we enter in the Command Window the following commands (without comments):

> > x = 0:pi/100:3; % create vector x with values between 0 and 3 with step % pi/100
> > y1 = sin(x);y2 = cos(x); % create vectors y1 and y2
> > plot(x,y1,x,y2, '--k') % create two lines in the same plot

The resulting plot is shown in Figure 3.3 on p. 54.

Figure 3.3 Two curves (sine and cosine) in a single plot.

3.1.1.2 The hold command option

This command is suitable when one plot already exists and we wish to add a new curve to it. In this case the hold on command must be typed and the new curve created by entering a new plot command. To complete the hold on process, we enter the hold off command that stops the hold process and shows the next graph in the new Figure Window.

As an example, the function sin(x2) can be added to an existing graph on Figure 3.3 by entering the additional commands:

> > y = sin(x.^2); % create new vector y
> > hold on  
> > plot(x,y, ':r') % add vector y to the same plot
> > hold off  

The resulting plot is shown in Figure 3.4.

Figure 3.4 sin x, cos x and sin x2 in a single plot.

3.1.1.3 Generation of several graphs on the same page

It is often necessary to place several plots on one page or, in other words, to multiply plots on the same page. For this, the subplot command has to be used, in the form:

subplot(m,n,p)

where m and n are the rows and columns of the panes into which the page is divided; p is the plot number that this command makes current (see Figure 3.5).

Figure 3.5 Possible arrangements of one page in four panes.

For instance:

subplot(2,2,4) creates 4 panes ordered in 2 rows and 2 columns and makes subplot 4 current
subplot(2,3,2) creates 6 panes ordered in 2 rows and 3 columns and makes subplot 2 current
subplot(2,1,2) creates 2 panes in the same column and makes the last subplot current
subplot(1,2,1) creates 2 panes in the same row and makes the first subplot current.

As an example, to generate four plots on one page: a one-point plot, sine plot, cosine plot and a circle in parametric form x = sin(t) and y = cos(t). The MATLAB® command to arrange the plots in this way are:

> > subplot(2,2,1) % makes pane 1 current
> > plot(0.1, 'p', 'MarkerSize',10) % plots a five-point asterisk with a size of 10
> > x = 0:pi/100:2*pi;
> > subplot(2,2,2) % makes pane 2 current
> > plot(x,sin(x)) % plots the sine curve
> > subplot(2,2,3) % makes pane 3 current
> > plot(x,cos(x)) % plots the cosine curve
> > t = x;
> > subplot(2,2,4) % makes pane 4 current
> > plot(sin(t),cos(t)); % plots the circle

The resulting plot is shown in Figure 3.6.

Figure 3.6 Four plots on the same page.

The plots can be arranged asymmetrically so that one can be placed into two or more columns or rows; for example, subplot(2, 2, [3, 4]) creates four panes on the page and spans the third and fourth panes by the third plot (the bottom of the current page).

3.1.2 Formatting of 2D plots

All plotting commands described above produce bare plots. In practice, a figure must have a title, grid, axis labels, suitable axes ranges and some text. The plot can be formatted by including the specifying commands in the program created or interactively by the plot editor in the Figure Window. The first method is preferable when we intend to use the written program repeatedly, and the second can be used when the figure created is intended to be saved, e.g. for demonstration.

3.1.2.1 Formatting of 2D plots with specify commands

The formatting commands must be entered after the plot command. Some of these commands are described below.

The grid command

The grid or grid on commands add a grid to the created plot. The grid off command removes the grid lines from the latticed plot. For example, typing grid in the Command Window immediately after the commands that produced Figure 3.3 will add the grid to the figure.

The axis command

Some possible forms are:

axis([xmin xmax ymin ymax])

axis equal

axis square

axis tight

axis off

The first command adjusts the axes to the limits written in the square brackets; the second sets the same scale for both axes, the third sets the axes region to be square, the fourth sets the axis limits to the range of the data to be performed on the plot, and the last removes the axis from the plot.

As an example, enter the sequence of commands:

> > x = 0:pi/100:pi; % pi is maximal value of x
> > y1 = sin(x);y2 = cos(x);
> > plot(x,y1,x,y2, '--k') % creates the plot with maximal x-axis limits 3.5
> > axis tight % sets xmax and ymax to pi and 1, respectively

Figure 3.7 will be produced by these commands.

Figure 3.7 Sine and cosine plot constructed with the axis tight command.

The xlabel, ylabel and title commands

These commands provide text to each of the axes and at the top of the plot. The text must be written in string form. The commands have the form:

xlabel('text string')

ylabel('text string')

title('text string')

There are provisions within the text string command to use Greek letters, font size and style, and also property options which can be written after the text for text angle, font name and color – see ‘Formatting of text strings’.

The text and gtext commands

Existing forms are:

text(x,y,'text string')

gtext('text string')

The first command directs the text at the point with coordinates x and y. The second command places it at the location chosen by the user. When this command is entered the Figure Window appears with two crossed lines and the user moves these lines by shifting the mouse to the point required, after which the text is entered by clicking on the mouse.

As an example, add the title, xlabel, ylabel, text and grid commands to the commands that constructed the plot in Figure 3.2:

> > title('Biomass vs Time')

> > xlabel('Time, min'),ylabel('Biomass,g/L')

> > text(75,6.3, 'Biomass') % display word Biomass at point with x = 75 and y = 6.3

> > grid

The resulting plot is shown in Figure 3.8 on p. 60.

Figure 3.8 Biomass data plot formatted with the xlabel, ylabel, title, text and grid commands.

The legend command

This be written as:

legend('text string1','text string2',..,'Location',location)

The command shows the line type for each plotted curve and prints explanations written in 'text string1', 'text string2', … The 'Location' property is optional and specifies the location for the legend. Thus, for instance, location = 'NorthEastOutside' places the legend outside the axes on the right; location = 'Best' places the legend inside the axes at a better location (least conflict with the data in the plot). The default legend location is in the upper right-hand corner of the plot.

So, for Figure 3.7 inputting the command legend('Sine', 'Cosine') generates Figure 3.9 (see p. 61).

Figure 3.9 Plot of the sine and cosine functions with legend.

Formatting of text strings

The text strings in the described commands can be formatted by writing modifiers (special characters) inside the string or by including the options PropertyName with PropertyValue in the command after the text string.

Typing modifiers can be font name, style, size and color, Greek letters, or sub- and superscripts. Some useful modifiers are:

\b \it \rm sets the type in bold, italic or normal, respectively.
– ^ sets subscript and superscript, respectively; for example, ′^oC′ sets superscript for ‘o’ and the resulting text will display as °C.
\name of the Greek letter sets a Greek letter; for example, \sigma sets σ and \Sigma sets Σ.
\fontsize{number} specifies the size of the type; for example, \fontsize{12} sets type at 12 point.
\fontname{name} specifies the name of the font used; for example, \fontname{Arial} sets Arial font.

By including in the command the property (the name and the value), we can regulate the color of the text, or its background. Some property names and their values are:

’Color’, ’color specifiers from Table 3.1 sets the text color; for example, Color, b sets blue color for the text string.
’BackgroundColor’, ’color specifiers from Table 3.1 sets the background color (rectangular area); for example, ’BackgroundColor’, ’y’ sets yellow color for background area.

A detailed explanation is available in MATLAB® Help and in the Text Properties section of the MATLAB® documentation.

3.1.2.2 Formatting of 2D plots with the plot editor

The Figure Window contains an assortment of formatting buttons and menu items. To start the plot edit mode the Edit Plot button, , must be clicked on the bar under the menu. The head of the Figure Window, with the menu and bar with most frequently used buttons, is shown in Figure 3.10. The properties of the axes and lines, and the whole figure, can be changed by means of the pop-up menu accessed by clicking the Edit Menu button. The title, axis labels, text and legend can be archived by means of the pop-up menu accessed by clicking the Insert button.

Figure 3.10 Plot Editor buttons in the Figure Window.

The text, legend and objects in the plot can be shifted by clicking on each of them. Double clicking of the curves, plotted points or axes activates the Property Editor, which allows many of the characteristics of the object selected to be changed or edited. Detailed information is available in the Help Window under the Editing Plot section.

3.2 Generation of XYZ plots

In MATLAB® there are three main groups of commands for lines, meshes and surfaces. These groups of commands together with various formatting commands are described below.

3.2.1 Generating lines in 3D plots

A line in 3D space connects points that are each described by three coordinates. Similar to 2D line plotting with the plot command, the plot3 command is used to plot a 3D line. In its simplest form this command can be written

plot3(x,y,z),

or in a more complicated form

plot3(x,y,z, 'Line Specifiers', 'Property Name', 'Property Value')

Here x, y and z are the equivalent vectors with point coordinates, and the Line Specifiers, Property Name and Property Value have the same sense as in the 2D case.

In 3D plotting the grid, xlabel and ylabel commands and, by analogy, the zlabel commands are also used.

For example, a 3D plot can be produced as follows:

> > t = 0:pi/100:6*pi;

> > x = t.*cos(2*t);

> > y = t.*sin(2*t);

> > z = t;

> > plot3(x,y,z, 'k', 'LineWidth',4)

> > grid

> > xlabel('x'),ylabel('y'),zlabel('z')

This algorithm computes the parametrically given coordinates x = t · cos(2 t), y = t · sin(2 t) and z = t with t changed from 0 up to 2π in steps of π/100.

With the procedure completed the plot appears in the Figure Window and the line has the form shown in Figure 3.11 (see p. 64).

Figure 3.11 Line in 3D coordinates.

3.2.2 Meshes in 3D plots

To understand the main commands of the 3D plot, namely mesh and surf, it is useful to understand the mesh construction in MATLAB®. As every point in 3D space has three coordinates x, y and z, these must be given in reconstructing a surface. In other words, we must have two 2D matrices for the x- and y-coordinates, respectively, and calculate the matrix of z-coordinates for every pair. The area of x and y coordinates for which the z-coordinates must be obtained is called the domain. Figure 3.12 overleaf shows an example of point representation in 3D space. The grid in the x–y plane is the domain given by the vectors x = − 2…2 and y = − 2…2. Each node in the x–y plane has a pair of x, y values. Writing all x-values ordered by rows, we obtain the X-matrix, and the same procedure yields the Y-matrix. For the case in Figure 3.12 matrices X and Y read:

Figure 3.12 Points in 3D interpretation and their x, y-plane projection.

With X and Y matrices available the z-coordinates must be obtained for every grid point using element-by-element calculation. The whole surface can then be plotted.

MATLAB® has a special command called meshgrid which allows us to calculate X and Y from the given vectors of x and y. The command has the form:

Here X and Y are the matrices of grid coordinates calculated in this function on the basis of the given x and y vectors that determine the domain dividing. When the x and y vectors are equal, this command may be simplified to

For example, the X and Y matrices can be defined for the case in Figure 3.12 as follows:

> > x = − 2:2;

> > [X,Y] = meshgrid(x)

The Z matrix can be calculated for instance via the expression

which corresponds to the Z-coordinates of the points in Figure 3.12.

We can now generate a graph with the aid of the mesh command. This command has the simplest form:

where X, Y and Z are the same matrices as in the example. This command also enables the mesh lines to be colored in the plot.

Thus, the program that plots the mesh surface reads:

> > x = - 2:0.1:2;

> > [X,Y] = meshgrid(x);

> > Z = X + Y + Y.*X;

> > mesh(X,Y,Z)

> > xlabel('x'),ylabel('y'),zlabel('z')

> > grid on

The resulting plot is shown in Figure 3.13.

Figure 3.13 Mesh plot.

3.2.3 Surfaces in 3D plots

If the plot needs colored surfaces between the mesh lines, the surf command can be used. The form of this command is:

X, Y and Z being the same matrices as in the mesh command.

Using this command for the function in the previous example we can enter commands as follows:

> > x = - 2:0.1:2;

> > [X,Y] = meshgrid(x);

> > Z = X + Y + Y.*X;

> > surf(X,Y,Z)

> > xlabel('x'),ylabel('y'),zlabel('z')

> > grid on

The resulting plot is shown in Figure 3.14.

Figure 3.14 Surface plot.

The surf command and the mesh command can be used in the form surf(Z) or mesh(Z). In this case the Z values are plotted versus the indices.

3.2.4 Formatting and rotation of 3D plots

Many of the 2D commands described in section 3.1.3 are suitable for 3D plot formatting, such as grid, title, xlabel, ylabel and axis. Many additional commands exist to format a 3D plot. Some of these are described below.

3.2.4.1 The box command

The box on command draws a box around the plot, and on entering the box off command the drawn box disappears. Figure 3.15 illustrates use of the command for the previous example:

Figure 3.15 Boxed surface plot.

> > box on

3.2.4.2 The colormap command

Color plays a more important role in 3D plotting than in 2D plotting. Default color is automatically generated with mesh or surf according to the z-values. With the colormap command the colors can be set at a constant value. The command has the form

colormap(c)

where c is a three-element vector in which the first element specifies the red color intensity, the second the green color intensity and the third the blue intensity (RGB); intensities are graded from 0 to 1, for example:

c = [0 0 0] black
c = [1 1 1] white
c = [1 0 0] red
c = [0 1 0] green
c = [0 0 1] blue
c = [1 1 0] yellow
c = [1 0 1] magenta
c = [0.5 0.5 0.5] gray.

If colormap([0 0 1]) is added to the commands that produced the mesh plot in Figure 3.14, we obtain this figure with the mesh lines changed in color to blue.

MATLAB® has some built-in color regulation commands in the form

colormap name

where the name can be jet, cool, winter, spring and a few others. For example, the colormap spring changes colors to shades of magenta and yellow.

3.2.4.3 The view command

The plot orientation relative to the viewer look can be regulated by the view command, having the form

view(az,el)

where az and el are azimuth and elevation angles: the first is the horizontal (x,y-plane) angle relative to the negative direction of the x-axes; the second is the vertical angle defined as the elevation from the x,y-plane. A positive az value is defined in the counter-clockwise direction. A positive el value corresponds to opening the angle in the direction of the z-axes. Both angles must be in degrees; default values are az = − 37.5° and el = 30°.

The azimuth and elevation angles are shown on Figure 3.16.

Figure 3.16 Viewpoint, azimuth and elevation.

The various planes can be viewed at the chosen appropriate angle, as follows:

• the x,y-projection of the 3D plot can be obtained with az = 0 and el = 90 – top view;

• the x,z-projection of the 3D plot can be obtained with az = el = 0 – side view [can be entered simply as view(2)];

• the y,z-projection of the 3D plot can be obtained with az = 90 and el = 0 – side view.

For example, the function can be plotted for default viewing angles, angles az = 20° and el = 35°, and for the top and side views by the commands:

> > x = − 3:.25:3;

> > [X,Y] = meshgrid(x);

> > Z = X.*exp(− X.^2-Y.^2);

> > subplot(2,2,1),surf(X,Y,Z)  % az = − 37.5 and el = 30°

> > subplot(2,2,2),surf(X,Y,Z)

> > view(20,35)  % az = 20° and el = 35°

> > subplot(2,2,3),surf(X,Y,Z)

> > view(2)  % az = 0° and el = 90° – top view

> > subplot(2,2,4),surf(X,Y,Z)

> > view(0,0)  % az = 0° and el = 0° – side view

which gives the result shown in Figure 3.17.

Figure 3.17 The function with viewing angles: (a) default, az = − 37.5° and el = 30°; (b) az = 20° and el = 35°; (c) az = 0° and el = 90° (top view); (d) az = el = 0° (side view).

3.2.4.4 The rotate3d command

This command has the forms

rotate3d on

rotate3d off

Typing the rotate3d on or pressing on the figure toolbar enables us to rotate the plot with the mouth pointer (see Figure 3.18), and to view the azimuth and elevation angles appearing in the bottom left-hand corner of the figure. For this, after entering the command in the Command Window, we should:

Figure 3.18 Plot in rotation regime with the rotate cursor and values of azimuth and elevation angles.

• go to the Figure Window,

• and by holding down the mouse button, rotate the plot and simultaneously view the az and el values that appear.

A view of the Figure Window in this regime is shown in Figure 3.18.

Plot rotation ends when the rotate3d off command is entered in the Command Window.

3.3 Specialized 2D and 3D plots

There are accessories to the commands of 2D and 3D graphics with the aid of which various specialized plots can be constructed. Some of them, such as errorbar, hist, semilogx and semilogy, are described below.

3.3.1 The errorbar command

Usually y = f(x) data have some uncertainty or an error in the y-value and it is desirable to set limits to it in each data point. The errorbar command makes it possible to plot the points with error limits. The simplest form of this command is

errorbar(x,y,e)

where x and y are the data vectors and e is that of the symmetric (two sides equal) error.

For example, if the biomass–time data used in section 3.1 have a symmetric error of ± 0.15 g/l at each point, the commands

yield, on entering, the plot shown in Figure 3.19.

Figure 3.19 Plot of biomass–time data with error bars.

A line style specifier may be included into the errorbar command for setting marker and/or line color and style; for example, in the previous plot the errorbar(t,b_mass,'--o') command changes the line to a dashed one and signs the data points with circles.

3.3.2 The hist command

The histogram is one of the most popular forms of data representation. Data values in a histogram are divided into intervals (bins) and plotted in the form of vertical bars, the height of which represents the number of data in each of them. Such a plot can be plotted with the aid of the hist command, whose simplest form is

hist(y) or n = hist(y)

where y is the vector containing the data points and n is the vector containing the number of data points in each of the 10 (default) equally spaced bins.

The first command yields the histogram plot, while the second yields only numerical output.

For example, the following data points are the weights (in mg) of 33 vials from a lot: 65, 80, 95, 93, 67, 81, 90, 93, 92, 83, 86, 83, 90, 94, 93, 96, 96, 100, 96, 50, 75, 81, 65, 88, 81, 60, 57, 61, 68, 77, 75, 76 and 70. A histogram can be plotted with the commands

> > y = [65 80 95 93 67 81 90 93 92 83 86 83 90 94 93 96 96 100 96 50 …

75 81 65 88 81 60 57 61 68 77 75 76 70];

> > hist(y)

> > xlabel('Weight,mg'),ylabel('Number of weights per one bin')

The resulting plot is shown in Figure 3.20.

Figure 3.20 Histogram of the weight data.

By using the n = hist(y) command in this example, the following frequencies will be displayed in the Command Window:

> > n = hist(y)

n =

1 2 3 3 2 3 5 4 6 4

The hist command has additional forms that allow us to set the required number of bins or output the values of the bin locations; detailed information is obtainable by entering help hist in the Command Window.

3.3.3 Plots with semi-logarithmic axis

In bioscience and engineering, coordinates with one axis on a logarithmic scale frequently have to be used. Log scales allow us to present values in a wider range than can be done with linear axes and can be used to plot the various exponential-like relationships that are widespread in biotechnology, such as expressions for first-, second- or higher-order reactions, population growth rates, radioactive decay and sterilization. For this, the semilogy or semilogx commands should to be used. In their simplest form:

semilogy(x,y, 'LineSpec'), semilogx(x,y, 'LineSpec')

The semilogy(x,y) command generates a plot with a log-scaled (base 10) y-axis and a linear x-axis; semilogx(x,y) generates a plot with a log-scaled x-axis and a linear y-axis; LineSpec (optionally) specifies the line type, plot symbol and color for the lines drawn in the semi-log plot.

For example, data for the decay of a radioactive substance are: 400, 200, 100, 50 and 25 disintegrations per minute at 0, 1, 2, 3 and 4 hours, respectively. A semi-logarithmic plot can be plotted with the commands:

> > N = [400 200 100 50 25]; t = 0:4;

% plot N on a semi-log y-axis and a linear x-axis

> > semilogy(t,N, '-o')

> > xlabel('Time, h'), ylabel('Disintegrations/minute')

> > grid

The data in the semi-log graph of Figure 3.21 are linear; these data on a normal graph produced by the plot command generate an exponential curve.

Figure 3.21 Semi-logarithmic graph for the relationship between time elapsed and residual radioactivity.

3.3.4 Additional commands for 2D and 3D graphics

Table 3.3 lists additional commands for 2D and 3D plotting; it gives the format of the corresponding basic command with short explanations, examples and the resulting plots. A complete list of 2D, 3D and specialized plotting functions can be obtained by entering the following commands: help graph2d, help graph3d or help specgraph.

Table 3.3

Additional commands and plots for 2D and 3D graphics*

*The commands are described in their simplest form.

3.4 Application examples

3.4.1 Two species concentration change in an opposed reaction

Consider that for two species AB with initial concentrations [A]0 and [B]0 equal 0.25 and 0 mol/L, respectively, the concentrations at time t = 0, 0.1, …, 1 can be determined from the following equations:

The reaction rate constants kf and kb equal 2 and 1 min−1, respectively.

Problem: Calculate and plot the results, and present two forms of a concentration–time graph: (a) [A] – t and [B] – t curves on the same graph, and (b) each of these curves on different graphs but on the same page.

The commands for presentation (a) are

> > % determine the initial concentration and reaction rate constants

> > Ao = .25;kf = 2;kb = 1;

> > t = 0:.1:1;  % create time vector

> > A = Ao*1/(kf + kb)*  % calculate concentrations A

(kb + kf*exp(−(kf + kb)*t));

> > B = Ao* kf/(kf + kb)  % calculate concentrations B

*(1-exp(−(kf + kb)*t));

> > % plot A and B vs. T; A-solid line (default), B-dashed

> > plot(t,A,t,B, '--')

> > % this and the next commands are used to format the current plot

> > xlabel('Time, min'),

> > ylabel('Concentration,g/l')

> > grid

> > title('Two species concentration vs time')

> > legend('The species A', 'The species B')

The plot created by these commands is shown at the top of p. 81.

For the second representation, the plotting command in the preceding text should be changed to:

> > subplot(2,1,1)

> > plot(t,A)

> > title('Concentration of the A-spacies vs. Time')

> > xlabel('Time, min'),ylabel('Concentration,g/l')

> > grid

> > subplot(2,1,2)

> > plot(t,B)

> > title('Concentration of the B-spacies vs. Time')

> > xlabel('Time, min'),ylabel('Concentration,g/l')

> > grid

After inputting, these commands produce two graphs on the same page:

3.4.2 Microorganism growth curve

One type of bacterium has a doubling time of 30 minutes. The initial number N0 of bacteria in a flask was 100. The general equation for bacterial population growth is

where n is the number of generations that have elapsed; in our case after 3 hours the number of elapsed generations n is 3 × 60/30 = 6.

Problem: Calculate and plot bacterial growth over 3 hours in two plots on one page: (a) N versus time t, and (b) log N versus t; mark calculation points with a circle (the letter ‘o’).

The commands used to solve this problem are:

> > No = 100;t_double = .5;  
> > n = 0:6; % Number of generations elapsed
> > N = No*2.^n;  
> > t = n* t_double; % Time elapsed
> > subplot(2,1,1)  
> > plot(t,N, '-o')  
> > title('Bacterial growth during 3 hours')  
> > xlabel('Time elapsed, hours'), ylabel('Number of bacteria, cell')  
> > grid  
> > subplot(2,1,2)  
> > semilogy(t, N, '-o') % the command for a semi-logarithmic plot
> > title('log plot of bacterial growth during 3 hours')  
> > xlabel('Time elapsed, hours'), ylabel('Log of number of cells')  
> > grid  

The figure produced is shown on p. 83.

3.4.3 Air density at atmospheric pressure

Under the ideal-gas approach, the density, ρ, of air at atmospheric pressure, p, and different temperatures, T, can be calculated with the ideal gas equation of state as:

with p = 101.3 kN/m3, R = 0.286 kJ/(kg°K) and T in degrees Kelvin.

The experimental air density values measured with mean two-sided error δρ = 0.1% at 0–100 °C (or 273.15–373.15 °K) are given in Table 3.4.

Table 3.4

Air temperature-density data

Problem: Calculate and plot the air density as determined theoretically (via the equation of state) and experimentally (from measured data with data errors in each point).

The commands to calculate and plot air densities at atmospheric pressure and temperatures of 273.15–373.15 °K are:

> > T = 273.15:20:373.15; % in °K
> > ro_exp = [1.293 1.205 1.127 1.067 1 0.946]; % in kg/m3
> > p = 101.3; % atmospheric pressure in kN/m2
> > R = .286; % in kJ/(kg°K)
> > ro_theor = p./(R*T);  
> > error = 0.001* ro_exp;  
> > errorbar(T, ro_exp, error,'.b'),hold on  
> > % plots experimental points with error bars  
> > plot(T, ro_theor), grid  
> > % plots theoretically calculated curve  
> > xlabel('Temperature, ^oK'),ylabel('Density, kg/m^3')  
> > title('Air Density'),  
> > legend('experimental points', 'theoretical curve'), hold off  

The resulting plot is:

3.4.4 Elliptic Gaussian function

The 2D elliptic Gaussian distribution function for the uncorrelated variables x (e.g. species weight) and y (e.g. species volume) having a bivariate normal distribution is given by

where μx and μy are the mean values of x and y, and σx and σy are the standard deviations.

Problem: Calculate and plot the distribution function for μx = 90, μy = 130, σx = 10, σy = 20 in x-axis limits μx ± 3σx and y-axis limits μy ± 3σy.

The steps used to solve this problem are:

• define 40-point vectors x and y with the linspace command;

• create an X,Y grid in the ranges of the x, y vectors by using the meshgrid command;

• calculate f via the above expression;

• generate an X,Y,f plot with the surf command.

The commands are:

> > sigma1 = 10; sigma2 = 20;

> > mu1 = 90; mu2 = 130;

> > setting limits for x, y, sigma_x, sigma_y

> > xlim = 3*sigma1; ylim = 3*sigma2;

> > xmin = mu1-xlim;xmax = mu1 + xlim;

> > ymin = mu2-ylim;ymax = mu2 + ylim;

> > x = linspace(xmin,xmax,40);  % defining vector x

> > y = linspace(ymin,ymax,40);  % defining vector y

> > create X,Y grid from the x,y vectors

> > [X,Y] = meshgrid(x,y);

> > f = 1/(2*pi*sigma1*sigma2)*exp(−((X-mu1).^2/(2*sigma1^2) + ( Y-mu2).^2/(2*sigma2.^2)));  % calculate element-wise the f

> > surf(X,Y,f)  % plot surface graph

> > xlabel('x'),ylabel('x'),zlabel('f(x,y)')

The configuration generated by these commands is shown on p. 86:

3.4.5 One-dimensional transient diffusion

The analytical solution of the 1D diffusion equation has the form

where c is a species concentration that changes with coordinate x and time t, and k is the diffusion coefficient.

Problem: Calculate and generate a 3D surface plot in which x and t are in the horizontal plane and the z-axis is the concentration; if D = 1 cm2/s, x = 0, 0.1, …, 2 cm, and t = 0.1, 0.2, …, 2.1 s, present the plot with azimuth and elevation angles of 123° and 26°, respectively.

The program follows these steps:

• create the vectors of x and t;

• create an X,Y grid in the ranges of the x and t vectors, respectively, by using the meshgrid command;

• calculate c for each pair of X and Y values by the above expression;

• generate a 3D plot by the X,Y,c-values determined;

• set required view point for the plot generated.

The commands to solve this problem are:

> > x = 0:.1:2;t = 0.1:0.2:2.1;  %D = 1 and is not used here

> > % creates the X,Y grid from the x,t vectors

> > [X,Y] = meshgrid(x,t);

> > % calculate element-wise the value of c

> > c = 1./sqrt(4*pi*Y).*exp(−(X).^2./(4*Y));

> > surf(X,Y,c);  % plot surface graph

> > xlabel('Distance');ylabel('Time');zlabel('Concentration')

> > view(123,26)  % az = 123° el = 26°

3.4.6 First-order reaction: concentration changes with time and temperature

In the first-order reaction A → Product, the concentration [A] changes with time t according to

where [A]0 is the initial concentration of A.

The rate constant k changes with temperature according to the Arrhenius equation:

in which a is a frequency factor, Ea is the activation energy, T is temperature and R is the gas constant.

The values of the parameters in these equations are a = 1014 s−1, Ea = 77 000 J/mol, T = 293, 294, …, 303 K, R = 8.314 J/(K mol), [A]0 = 0.25 mol/L, and t = 1, 1.2, …, 4 s.

Problem: Calculate and generate a 3D surface plot of the concentration (z-axis) as a function of time (x-axis) and temperature (y-axis). Present the plot with azimuth and elevation angles of 67° and 22°, respectively. Change colors to the autumn combination of colors.

The following steps should be used to solve this problem:

• determine a, Ea, R and [A]0 and create the vectors for t and T;

• create an X,Y grid in the ranges of the t and T vectors, respectively, by using the meshgrid command;

• calculate k and [A]0 for each pair of X and Y values by the above expression;

• generate a 3D plot using the defined values of X, Y and [A];

• set the required view point and color combination.

The commands are:

> > ea = 77e3;a = 1e14;R = 8.314;a0 = .25;

> > t = 1:.2:4; T = 293:303;

> > [X,Y] = meshgrid(t,T);

> > k = a*exp(- ea./(R*Y));

> > A = a0*exp(- k.*X);

> > surf(X,Y,A)

> > xlabel('Time, sec'), ylabel('Temperature, K'),

> > zlabel ('Concentration, mol/L'),

> > view(67,22)

> > colormap spring

The resulting plot is shown at the top of p. 89.

3.5 Questions for self-checking and exercises

1. The plot3 command generates on a 3D graph: (a) a surface, (b) a surface mesh, (c) a line?

2. Which sign or word specifies a circled point on a plot: (a) +, (b) the word ‘circle’, (c) o, (d) t?

3. The command subplot is intended to: (a) place a new plot into a previously generated plot, (b) place several plots on the same page, (c) change colors on the graph?

4. A biotechnology company is using a bacterial broth to produce an antibiotic; the pH of the broth was measured for 5 hours every 30 minutes, giving: 6.12, 5.13, 5.84, 6.53, 6.12, 6.3, 6.04, 5.79, 5.94, 6.03, 6.12. Plot the data with a solid line between data points marked with a diamond; add axis labels, grid and a caption.

5. Generate three plots on one page for the following geometrical figures.

Astroid:

with 0 ≤ φ ≤ 6π.

Archimedes’ spiral

with 0 ≤ φ ≤ 6 π.

Leminscate

with 0 ≤ φ ≤ 2 π.

Place the astroid and Archimedes’ spiral in the first graph lines and the lemniscate in the second line. Make the size of the current axis square for the astroid and Archimedes’ spiral, and add a caption to each plot.

6. Plot the function

for 0 ≤ x ≤ 4.

Add to the graph the new function

with 1 ≤ x ≤ 5.

Add axis labe–ls, grid and legend to the graph.

7. The kinematic viscosity η (in cP) of a bio-oil as a function of temperature T (in K) is given by the Andrade-type expression

where A = 0.0005806 and B = 3382.

Measured viscosities at each temperature are given in the table below

The uncertainty of the experimental data is ± 0.5%.

Plot the calculated and measured viscosity values. Give the error bars at each measured point, add axis labels, a caption and legend.

8. Plot the surfaces of the following bodies on one page. Make the current axis equal in size and add a caption to each figure:

(a) Ellipsoid

where a = 1.2, b = 1.7, c = 2.1, u = 0, …, 2π, v = 0, …, π.

(b) Cross-cap

where u and v range between 0 and π.

9. The turbulent-flow Nusselt number (Nu) for a plate of length l is a function of the Reynolds, Re, and the Prandtl, Pr, numbers:

where 5 × 105 ≤ Re ≤ 107 and 0.6 ≤ Pr ≤ 2000; use 100 values of Re and Nu each.

Plot log(Nu) as a function of x = log(Re) and y = log(Pr). Add axis labels and a caption.

10. The height of common prairie aster (in centimeters) plants was measured in the field, giving: 151, 174, 146, 155, 105, 122, 144, 154, 130, 151, 151, 154, 125, 140, 182, 160, 127, 141, 122, 140, 160, 157, 126, 129, 138, 156, 141, 182, 131, 156, 180, 180, 133, 162, 141, 189, 147. Plot a histogram for these data and add a title to the graph. Calculate the mean and standard deviation and write the determined values into the graph (with the text commands; take the initial text coordinates as follows: x = 150 and y = 7 for the mean, and x = 150 and y = 6.5 for the standard deviation).

11. Weight–height data for students in an academic group are (see Subsection 2.2.5.4): h = 155, 175, 173, 175, 173, 162, 173, 188, 190, 173, 173, 185, 178, 168, 162, 185, 170, 180, 175, 180, 175, 175, 180, 165 cm; w = 54, 66, 66, 71, 68, 53, 61, 86, 92, 57, 59, 80, 70, 59, 50, 145, 68, 78, 67, 90, 75, 74, 84, 53 kg. The data were fitted by the polynomial expression wf = 906.14 – 11.39 h + 0.03780 h2. Plot the wf(h) curve calculated by the fit expression and the data points. Take the heights for the fit expression in the interval between the minimal and maximal h-data values. Add axis labels, grid, a caption and a legend to the graph.

12. The M2 molar concentration of a solute in water is calculated by the expression M2 = M1V1/V2 (see Subsection 2.2.5.4), where V2 is the final solution volume, and M1 and V1 are the initial solute concentration and solution volume, respectively. Given the initial solution volume V1 = 0.1 L, plot the 3D graph M2(M1,V2) in which M1 changes in the range 0.5–2 mol/L and V2 in the range 0.1–0.9 L. Show the 3D graph with azimuth and elevation angles of − 135° and 35°, respectively. Add axis labels, a grid and a caption to the graph.

13. The partial pressure, PA, of the gaseous A-component in the reaction A → B + C can be calculated by the expression PA = P0ekt where the initial partial pressure P0 = 55 Torr, the reaction rate constant k = 3 min−1 and the time t changes from 1 to 5 min in steps of 0.25 min. Plot two graphs on one page: ln(PA) as a function of t, and PA(t). Add axis labels, a grid and a caption to each graph.

14. The molecular velocity distribution of gaseous helium is a function of the velocity v and temperature T and can be described by the Maxwell–Boltzmann equation:

where the molecular weight M is 0.004 kg/mol, the gas constant R = 8.314 J/(mol K), and v and T are in the ranges 0–1300 m/s and 50–500 K, respectively. Plot the graph fv(v,T) with axis labels, a grid and a caption.

15. In bio-fluid mechanics the motion of a rigid particle is described by time-dependent coordinate equations. Plot the trajectory of a particle with coordinates:

Use the LineWidth property equal to 5. The time t is given in 0–2 dimensionless units. Add axis labels, a grid and a caption.

16. The relationship between the surface coverage θ of an adsorbed gaseous substance on the substance pressure P above the surface at constant temperature (Langmuir isotherm) is:

where b is constant for the given substance at given temperature. Plot the graph of θ(P,b) when P and b are in the range 0–1 and 0.1–100, respectively. Add graph axis labels, a grid and a caption.

17. The surface roughness height h of bio-films used in ophthalmology (lenses) and other medical applications can be modeled by the equation

where x and y are dimensionless surface size and change in the range 0–1 at steps of 0.01; k is the wave number and is equal to 5.

Plot a rough surface graph h(x,y) with azimuth and elevation angles of − 24° and 82°, respectively; add axis labels, grid and a caption to the graph.

3.6 Answers to selected exercises

3. (b) placing several plots on the same page. 5.

8.

12.

16.