221: Mech Computer Lab 1
Hand in the solutions to the three questions in the lab at the end of the lab. Success in many kinds of engineering requires skill at numerical approximation. The computer labs this term build this skill. This lab in particular will build your understanding of and experience with numerical integration techniques. Some MATLAB commands you already know will be reinforced and some new commands introduced. Lab Learning Goals: After completing this lab, you should be able to use MATLAB to: generate vectors and dene functions create plots using either the Command Window or the plot Graphical User Interface write a simple program, save it as an .m le, run it, edit it, and run it again perform numerical integration to a given accurancy Question 1: Dening functions and plotting Let us review how to dene and plot functions in MATLAB. We will also show you how to interface with the plot GUI (Graphical User Interface). You may want to refer to the list of MATLAB commands given in the pre-lab in case you are not familiar with any of the commands used in this lab. Let us start with a few examples. First, we will plot the function f (x) = x2 over the interval [0,1]. Here are some instructions on how to do that: Create a vector x of length 11 containing values between 0 and 1 using the linspace command as follows: >> x = linspace(0,1,11) The vector x corresponds to the end points of 10 equally spaced subintervals of [0,1]. 1
Now we will map these points to the function f (x) = x2 and store the values in the vector y. Enter the following command: >> y = x.^2 Dont forget the dot in the command above, which tells MATLAB to take the square of each entry of x. The command y= x.*x also works, but x*x will give an error message (it tries to use matrix multiplication and the dimensions of the vector x do not allow this). We now have dened our function on the interval. Its time to plot it! We will plot it as a thin red line, with our data points indicated by xs. Type in the following command: >> plot(x,y,r-x) You should now have a window on your screen that contains a plot of y = x2 from x = 0 to x = 1. The next example consists in plotting two dierent functions on the same axis. Specically, we will plot the functions f1 (x) = sin(x) and f2 (x) = sin(x2 ) from x = 0 to x = 2 , using equally spaced points on the x-axis. The distance between successive points on the x-axis is set to 0.01 units. Here are some instructions on how to create these plots: First of all, we clear the workspace and the gure window. To do that, enter the following commands: >> clear all; clf; The clear all command clears all existing variables and other items from the workspace; it also frees up system memory. The clf command clears the current gure window. This practice of clearing the workspace is highly recommended whenever you start a fresh new set of operations in MATLAB. Now create a new vector x containing points between 0 and 2 that increase constantly by 0.01. Since you want to control the step size between vector components, the colon operator : is a better choice to generate the vector x than the linspace command, as shown below: 2
>> x = 0:0.01:2*pi; Now dene the functions by creating two vectors y1 and y2: >> y1 = sin(x); >> y2 = sin(x.^2); You can now create the plots by entering the following: >> >> >> >> plot(x,y1,k-,linewidth,2); hold on; grid on; plot(x,y2,m--,linewidth,2); legend(f1,f2);
Take note of what you see in the plot window. Consider what would have happened if we had not included the hold on command. The example above shows how to create plots by specifying command line options. It is also possible to modify the appearance of plots after they have been generated using GUI commands in the plot window. Now that you know how to dene and plot functions, work through the following questions and hand in your answers. Question 1.1: Using any method you prefer, plot the following function on the given interval. Use 21 points (i.e. 20 subintervals) on the x-axis. Plot it as a solid red line of greater than default thickness. Add a grid and give the plot the title Gaussian. Print out this plot and hand it in as part of the lab. f (x) = ex ,
2
0x1
Note: the MATLAB command for the exponential function, ex , is exp(x). Question 1.2: Plot on the same axis the following functions on the interval 0 x 1, f1 (x) = sin(2x), f2 (x) = cos(2x) Let the points on the x-axis be spaced 0.025 units apart. Plot f1 (x) using black asterisks with no connecting line, and f2 (x) using black circles with no connecting line. Label your x- and y -axis, x and y respectively. Include an appropriate legend. Print out this plot and hand it in as part of the lab. 3
Question 2: Trapezoidal Method and the Gaussian A very important integral which appears in many dierent disciplines is the Gaussian ex integral: dx
2
(1)
In most cases, if you tried to evaluate this quantity analytically you would wind up with a headache. Fortunately we have numerical integration at our disposal. In the following questions you will use MATLAB to estimate a Guassian integral. While the implementation of the numerical method can be straightforward if you use MATLAB built-in functions, it is important that you learn a way to estimate how accurate your numerical calculations are. Question 2.1: Use the trapz command to estimate I = ex dx using 20 0 subintervals. What is the approximate value you obtain? Remember that N subintervals requires N + 1 points. If you are unfamiliar with the trapz command, review the list of MATLAB commands provided in the pre-lab or consult the MATLAB online help by typing help trapz at the command line. Hand in your approximate value of I . Question 2.2: As reviewed in the pre-lab, a good rule of thumb to check if the outcome of your numerical integration is accurate is to double the number of subintervals, N , used to evaluate the integral and observe if the digits of your calculation change. When using a second order or higher method, digits that do not change as N is doubled are usually correct. With this in mind, work through the following problem: Consider the integral I dened in Question 2.1, let TN be the estimate of I computed using the Trapezoidal method. Fill out the following table until you get the value of TN accurate to 3 decimal places. N TN 2 4 8 16 . . . . . . 4
1
2
What value of N is required to reach that accuracy? What is the approximate value for I ? Hand in the answers to these questions along with the table you constructed. Question 2.3: Remember from class that if I is the value of an integral and TN is an approximated value of I computed using the Trapezoidal method with N subintervals, then I TN C/N 2 for N large enough and constant C . This means that I T2N C/4N 2 . We can use this to get that TN T2N 43C2 . Determine an estimate for C and use it N to give an estimate on the accuracy of T20 for the integral dened in Question 2.1. Hand in your estimate of C and the calculations used to get it, and the estimated error in T20 .
Question 3: Introduction to .m les and Simpsons Rule To solve large problems in MATLAB it is impractical for you to continuously input commands in the Command Window. A more ecient way is to create and execute an .m le, that is, to write all the sequence of commands to be executed into a text le, save it with the extension .m in your working directory, and then run it in MATLAB. The sequence of commands in an .m le is known as a script. Once written, the script can be run (and modied if needed) repeatedly with little eort. To show you how to create and run .m les, let us start by creating a le with all the commands you used to answer Question 2.1. You can then modify your .m le as needed to work through the nal part of this lab. Here are some instructions to create and execute an .m le: Create a new m.le (File New M-File). If there is no pull-down menu at New, try rst clicking on the Command Window. Type in all the sequence of commands you used to estimate the integral in Question 2.1. Remember to include a semi-colon ; at the end of any line of code whose output you would like to suppress once you run the le. Save the le in the same directory as your Current Directory with the name traprule.m. Be sure to add the .m extension.
5
Go to the MATLAB command window and enter >> traprule MATLAB will now run the script contained in traprule.m. It is also possible to run the script from the editor window (use the button with the green arrow). If no errors appear, your code has been executed. You can easily verify that running traprule.m produces the same result as executing each line of code separately as you probably did when you worked through Question 2.1. Keep in mind that .m les should be saved in your working directory and the lename should not contain spaces. Now that you know how to create, save, and run .m les, you are ready to work through the last part of this lab. Question 3.1: Create an .m le containing all the necessary commands to compute an estimate of the integral of f (x) = x2 sin(x) on the interval 1 x 2, using Simpsons Rule on N = 20 subintervals. Recall that you wrote a similar code as part of your pre-lab assignment. If you are having trouble with errors when you run your .m le, take note of where the error is occurring in your le (MATLAB will tell you) and attempt to debug the problem. When your code runs without errors, check that it is computing the correct quantities by comparing to the exact value of the integral or approximate values computed with trapz. What value did your Simpsons Rule code give you for the approximate value of the integral? Write down the answer and submit it along with a print-out of your Simpsons Rule code.
6