CS310 Team Lab #2
Functions and Plotting in MATLAB
Numeric Computing

OBJECTIVES

INTRODUCTION

MATLAB is a "computing environment" that allows the user to type mathematical expressions into a Command Window and immediately see the resulting answers. But MATLAB is much more than a graphing calculator for your computer. MATLAB also allows the user to define their own functions that can be used the same way that built-in functions like sin, cos, and sqrt are used.  In this team lab, we will begin to explore the powerful features of MATLAB by using some of the vast number of internal functions of MATLAB, defining our own specialized functions, and plotting the results of invoking or calling our newly defined functions.

PROBLEMS

1. Use MATLAB's Command Window to compute each of these results.

  1. What is the decimal approximation of sine of pi over three? How do you enter the value π in MATLAB?
  2. What is the exponential of the sine of pi over three?  Compute the exponential of the result from part a.  Use the default variable name that MATLAB assigned or repeat the previous command and assign the result to a variable. Repeat the calculation on a calculator to check the results.
  3. Compute the exponential of the sine of a range of angles from 0 to pi over three in (pi over three/100) increments.

    This is a tedious calculation to do with a hand-calculator, unless you program it. But with MATLAB, this is calculation is completed in one (or three) easy steps. This calculation can be done in one command. We split it up here to help you see the results of each part of the computation.

    1. Use the colon operator ( : ) to create a vector of values that ranges from 0 to pi over three in increments of (pi over three/100) and assign the vector result to the variable named vec.
    2. Compute the sine of that vector and save it as sin_vec.
    3. Compute the exponential of the vector result and save it as exp_vec.

    The sin and exp functions evaluate all of the values that are contained in the variable that is passed to the function.  Many predefined functions in MATLAB can accept either scalar variables, vector variables as shown here, or even a matrix of values. 

    View the contents of each variable you have defined by double-clicking on the variable's name in the Workspace window.  The Variables window shows the contents of each variable you select.

  4. Enter the following commands in the Command Window and compare each result with the contents of vectors vec, sin_vec, and exp_vec, respectively.  Determine and describe to each other what the parenthesis notation in each statement means.

    >> val1 = vec(2)
    >> val2 = sin_vec(1:5)
    >> val3 = exp_vec(2:3:20)

  5. Enter the following statements which clear all current variables and recompute values for the range of angles from 0 to two pi.  Use the semi-colon to suppress the output of these commands

    >> clear
    >> vec = 0 : pi/200 : 2*pi;
    >> sin_vec = sin(vec);
    >> exp_vec = exp(sin_vec);

    When you type a semicolon at the end of the command, MATLAB does not echo the answer to the screen. What would be the value of this, if you can’t see your result?

  6. Enter the following command to plot sin_vec vs. vec as red circles.

    >> plot( vec, sin_vec, 'ro' )

    You can now see a plot of your results. Notice how MATLAB plots pairs of points defined by the vectors vec and sin_vec (you may need to enlarge the Figure window to see the individual points). It is much more common to have MATLAB plot the points as a smooth line (by drawing lines between the points). Modify your previous plot command so that it plots sin_vec vs. vec as a red line:

    >> plot( vec, sin_vec, 'r-' )

    Enter the following command to plot exp_vec vs. vec as a blue dotted line.

    >> plot( vec, exp_vec, 'b:' )

    Enter the following command to plot both curves on the same plot figure.

    >> plot( vec, sin_vec, 'r-'  ,  vec, exp_vec, 'b:' )

    The characters in quotes specify the format for each set of points.  Try help plot for more format options.

The Command Window is useful for short calculations like this, but it is more convenient to define a new function if you are doing a complex multi-step calculation.

SHOW YOUR PLOT WITH BOTH CURVES TO A LAB LEADER

2.  Defining your own MATLAB function.

Defining your own functions shows more of the real power of MATLAB. You will now create your own unique function that computes the exponential of the sine of a value that is passed to your function as a parameter.

  1. Within your MATLAB folder in the My Documents folder on your computer's L: drive, create a new folder named TeamLab2. You will save the new function definitions that you define for this team lab here.

  2. Create a new MATLAB function file by selecting the menu option New→Function (on the HOME tab). Dock the Editor window, if it is not already docked.  A quick and easy way to do this is to type Ctrl+Shift+D (alternatively, you can find the Dock editor option by selecting the icon that looks like a point-down triangle in the upper-right corner of the Editor window immediately to the right of the ? icon). Edit or enter the following lines into the Editor window.

    function results = exp_of_sin(value)
    s_val = sin(value);
    e_val = exp(s_val);
    results = e_val;
    end

    Save this m-file with the name exp_of_sin.m. This will be the name MATLAB wants to give it when you click on the Save icon (on the HOME tab).  Be sure to navigate to the MATLAB folder you created in part a and save the file in that folder.

  3. In the Command Window, try your new function by entering

    >> cd TeamLab2
    >> exp_of_sin(pi/3)

    If you get an error message when you run the cd command, try the following instead: cd L:\My' Documents'\MATLAB\TeamLab2

    The command cd tells MATLAB to change the Current Folder; cd stands for change directory and a directory is just another name for a folder. Note that when the second command is executed, it returns only the final result of 2.3774. Note also that you simply enter the function name with an input value, just as you do when calling built-in functions like sin or exp.

  4. Modify (and save) the function definition by changing the last line to

     results = [ s_val , e_val ];

    In the Command Window, use the up arrow button to repeat the previous command.  Note that both results are returned in a single row vector.  The sine value is listed first and then the exponential of the sine.  The comma indicates that each value is one element in a row.  This type of return value is rarely useful. (Think about why this might be the case.)

  5. Modify (and save) the function definition by changing the last line to

              results = [ s_val ; e_val ];

    In the Command Window, use the up arrow button to repeat the previous command.
    The results are now in a matrix with a single column.  The sine is listed in the first row and the exponential of the sine is in the second row.  The semi-colon indicates that each value is in its own row of the result matrix. This type of return value (a column vector of values) is more useful.

  6. Assign the value pi/3 to a variable, such as angle, and then use that as the argument of a call to your exp_of_sin function.

    >> angle = pi/3; 
    >> soln = exp_of_sin(angle)

    Your own function works the same way as built-in functions. You can assign its results to a variable that you can later use for additional calculations.

  7. Use your function to compute the sine and exponential of a vector (a list) of values. Enter the following in the Command Window:

    >> angle = [ 0 : pi/100 : 2*pi ]; 
    >> soln = exp_of_sin(angle)

    Note that the results are returned as a matrix. The first row contains the results for the sine value and the second row contains the results for the exponential value. This is the matrix equivalent to the result that you obtained for a single input value, but now with multiple input values, you are returned multiple output results.

  8. Separate the sine and exponential rows into individual variables using the commands

    >> s_values = soln(1, :);
    >> e_values = soln(2, :);

    What does the notation (1, :) mean?  What would soln(:, 1) return?  Make a guess and then try it from the Command Window to see.

    Plotting results in MATLAB always requires that you evaluate the function that you want to plot at closely spaced points and then plot the point pairs. We will do this many times as we learn to use MATLAB. Try these commands to plot. Be sure that you understand what each command line does.

    >> plot(angle, s_values, 'g', angle, e_values, 'k')
    >> title('sin and exponential of sin for angles 0 to 2\pi')
    >> xlabel('angle (radians)')
    >> ylabel('sin(angle) and exp(sin(angle))')

  9. Define a new function exp_of_sin_2 that returns the sine and exponential of sine as two separate variables instead of two rows in a single variable. The function header changes to:

    function [s_val , e_val] = exp_of_sin_2(value)

    and the function body will be:

    s_val = sin(value);
    e_val = exp(s_val);

    Save this m-file with the name exp_of_sin_2.m and then call it from the Command Window with this command:

    >> [sineVals, expVals] = exp_of_sin_2(angle)

    Review the contents of the Workspace to see that you now have two separate variables, sineVals with the sine of each value and expVals with the exponential of the sine of each value.

3. Grade Percentage Function

  1. Define a function named calcPercent that computes a weighted percentage based on several vectors of scores passed as arguments to the function. The function must be named calcPercent; must use sum() and length() and not mean(); and the function header must be as follows:
    function final_percent = calcPercent( exams, quizzes, hws, labs )
    % CALCPERCENT Computes a weighted percentage for the grade
    % components entered.
    %
    % AUTHOR: Beck Hasti
    % 
    % exams is a vector of exam scores with a max of 60 pts per exam
    % quizzes is a vector of quiz scores with a max of 10 pts per quiz
    % hws is a vector of homework scores with a max of 5 pts per hw
    % labs is a vector of lab scores with a max of 2 pts per lab
    % 
    % final_percent = exams*60% + quizzes*13% + homeworks*15% + labs*12%

    Note that all you've entered so far is the function header and the comments describing what the function is supposed to do. You will now need to create the function body by adding the code that actually does the computation described in the comments.

  2. Create a script to test your calcPercent function: create a new (empty) script by selecting New Script on the HOME tab, or by selecting New→Script on the HOME tab or the EDITOR tab. Put a clear command at the top of this script and save the script as TeamLab2.m

  3. Next add code to your script to test your function:
    1. Test your function with these values.
      • calcPercent( 0, 0, 0, 0 )
      • calcPercent( 60, 10, 5, 2 )
      • calcPercent( 48, [10 9.5 6.8 9.8 10], [5 5 4.5 3.5 3 5], [2 2 2 2 2] )
      Make sure to run your script after each line of code you enter to check what value your calcPercent is returning.
    2. Calculate the percent for a student with the scores given below. In your TeamLab2 script:
      1. create a vector for each type of grade component:
        • exam scores: 54, 42, 48
        • quiz scores: 10, 9.5, 10, 10, 7.8, 6.5, 8.5, 7, 10, 9.2, 4.5, 8.75
        • homework scores: 5, 5, 4, 5, 5, 3, 4, 5, 4.5, 4, 4, 5, 5, 4.5, 4
        • lab scores: 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2
      2. add code to call your newly defined function and save the results in a new variable.
      3. add code to display your result value in a statement like: The final grade percentage is 84.4229%

It is possible to publish functions to HTML (or PDF) just like you did with the script from Team Lab 1. In this course, you will not be required to publish your functions. Instead, we will require you to turn in your function m-files so that we can run our test scripts with your functions.

SHOW YOUR CALCPERCENT FUNCTION AND RESULTS OF PART C TO A LAB LEADER

4.  More complex and realistic functions

The previous examples of sine and exponential are familiar but they have no physical meaning. Following is a formula for the reaction rate of thermonuclear fusion between the hydrogen isotopes of deuterium and tritium as a function of their temperature measured in units of keV. This formula is valid for 1<T<25 keV. Don’t be concerned about the meaning of it, it could as easily be a formula for chemical reactions or combustion reactions.

fusion power formula

  1. Define a new MATLAB function that returns the results of this calculation when value(s) of T are provided as input.  Name the function DT_reaction_rate. The function header is:

    function rr = DT_reaction_rate(T)

  2. Test your new function by entering the following in the Command Window.

              >> DT_reaction_rate(4)

    The result should be         5.1183e-18
    If this is not the result, debug your function's formula before proceeding.

  3. Now evaluate the formula for values of T between 1 and 25 in increments of 0.1. Be sure to use a (;) after entering the commands so you don’t get a lot of output in the Command Window.  Assign this result to a variable like DTrr.

    What happens?  Do you get an error message?  If you do, it is probably because the arithmetic operators in the function DT_reaction_rate are trying to do matrix arithmetic. This is not what you want.  Use element-wise operators so that they do their computation on each element of the input vector and do not try to do matrix arithmetic. Plot the results that you obtained.

  4. We now wish to compute the fusion power density produced by combining deuterium and tritium together and heating them to high temperatures like 1 to 25 keV. Write another function named DT_fusion_power that evaluates the following formula.

    fusion power formula

    where nD and nT are the densities of the deuterium and tritium reactants measured in atoms / cm3.

    The function header for the fusion power function is:

    function power = DT_fusion_power( nD, nT, T )

    As you enter the formula for power, note that DT_reaction_rate(T) replaces the expression reaction rate term in the formula above. So now you have a function calling another function. This is modular programming! Modular programming is where one programming unit (function) calls another programming unit to solve an even bigger problem.  Very complex problems are solved by combining the solutions of many smaller problems.

  5. In the Command Window enter

    >> DT_fusion_power( 2.7e19, 2.7e19, 4 )

    You should get the result               2.0895e+09

    These reactant densities are gas densities at atmospheric pressure and the power density is 2 Megawatts per cubic centimeter!!! Notice that you can write functions that have several input quantities. This function has three inputs and the reaction rate function has only one input.

TRY THESE PLOTS WHILE YOU WAIT TO SHOW US YOUR RESULT

5. Cool Plots

Plot the following graphs using the commands and functions provided.  Between each exercise use the clf command to clear the previous plot.

  1. Fast Fourier Transform function in MATLAB.  This is a very powerful function used to fit periodic data.  The eye function forms the identity matrix, all ones on the diagonal.  Enter in the Command Window:

    >> plot(fft(eye(17))), axis equal, axis off

  2. Epicycloid.  Cut and past the following function into an empty function definition document and save it with the name cool_plot1.m

    function xx = cool_plot1(a,b)
    t=0:0.05:10*pi;
    x=(a+b)*cos(t)-b*cos((a/b+1)*t);
    y=(a+b)*sin(t)-b*sin((a/b+1)*t);
    plot(x,y)
    axis equal
    axis([-25 25 -25 25])
    grid on
    title('Epicycloid: a=12, b=5')
    xlabel('x(t)'), ylabel('y(t)')
    xx = 1;


    Call the function from the Command Window as:

    >> cool_plot1(12, 5)

  3. Open a new document and cut and paste the following function into it and save it.

    function xx = cool_plot2(r)
    t=-5:0.005:5;
    x = (1+t.^2).*sin(r*t);
    y = (1+t.^2).*cos(r*t);
    z = t;
    plot3(x,y,z)
    grid on
    xlabel('x(t)'), ylabel('y(t)'), zlabel('z(t)')
    title('\it{cool plot 2 example}','FontSize',14)
    xx = 1;

    Call the function from the Command Window as:

    >> cool_plot2(20)

ADDITIONAL PROBLEMS or FURTHER READING

Challenge problem: Modify the calcPercent function so that it correctly calculates the grade percentage for this semester of CS 310. This means you will need to use the appropriate weighting of the grade components, take into account the fact that your lowest quiz and team lab scores get dropped, and figure out how to handle the computation of the homeowork portion of the grade. To make things a bit simpler, you may assume that there are always at least two team lab scores and two quiz scores.

Concluding remarks: Being able to define your own functions is at the core of all functional programming. A common technique for solving large and complex problems is to break them into smaller problems that can be solved independent of the larger problem. Each smaller problem is likewise solved and code written to accomplish that task. Once all of the smaller problems are solved (and their solutions are tested), we can then solve the larger problem by calling each functional part of the larger solution in the correct sequence.

Good modularization of our solutions often means that our functions are useful for other problems and can thus be reused in different ways to solve different problems. Can you think of any large problems that you were able to solve by breaking into smaller problems? What are the pros of defining a function instead of simply typing that code into the "main" solution script? What are the cons of defining functions?