## Arc Length and Functions in Matlab

Consider the parametric equations

$\begin{eqnarray} x&=&2 \cos t\\ y&=&3 \sin t \end{eqnarray}$

on the interval [0,2 pi]. To calculate the length of this path, one employs the arc length formula.

$L=\int_0^{2\pi}\sqrt{(dx/dt)^2+(dy/dt)^2}$

However,

(dx//dt)^2=(-2 sin t)^2=4 sin^2 t

and

(dy//dt)^2=(3 cos t)^2=9 cos^2 t.

Hence,

$\begin{eqnarray} L&=&\int_0^{2\pi}\sqrt{4\sin^2 t+9\cos^2 t}\\ L&=&\int_0^{2\pi}\sqrt{4(1-\cos^2 t)+9\cos^2 t}\\ L&=&\int_0^{2\pi}\sqrt{4+5\cos^2 t} \end{eqnarray}$

Because this last integral has no closed-form solution, we will need to apply some numerical routine (such as Simpson's rule) to obtain a decimal approximation for the integral. We are going to employ Matlab's quad command for this purpose, but first we must digress and learn how to write functions in Matlab.

### Directory Structure

If you work on a computer in PS116, open MyComputer and browse to your Documents folder. There you should create a Math50C folder (no spaces in filenames), then create a subfolder named Matlab in the Math50C folder. In the Matlab folder, create another subfolder names ArcLength (again, no spaces in filenames).

Once you've complete your directory structure as described in the previous paragraph, start Matlab then change the current working directory to the ArcLength folder. The easiest way to do this is to click the three-button icon directly to the right of the navigation edit box on the Matlab toolbar, then browse to the ArcLength folder. Check that your current directory points to the ArcLength folder by reading the contents of the Navigation box or executing the command pwd at the Matlab prompt.

If you are working at home, you want to start the same sort of directory structure.

• On a Mac, in your Documents folder, create a Math50C folder, then a Matlab folder inside the Math50c folder, then an ArcLength folder inside the Matlab folder.
• On a PC running Windows, in your My Documents folder, create a Math50C folder, then a Matlab folder insider the Math50C folder, then an ArcLength folder inside the Math50C folder.

Note that the directory structure described above is only a recommendation. You are perfectly free to create your own names and structure. However, don't fall into the trap of dumping all of your work into a single folder. You will regret such a decision as the number of files begins to grow through the course of the semester.

### Function M-files

Very Important: Change the current directory to point at the ArcLength folder. Check this with the pwd command at the Matlab prompt.

Open Matlab's editor. There are several ways that you can open the editor.

1. You can select File->New->M-file from Matlab's menu.
2. You can click the New M-file icon on Matlab's toobar.
3. You can type edit at the Matlab prompt.

Of these three options, the last is our favorite. When the editor opens, enter the following code:

function y=f(t)
y=sqrt(4+5*cos(t)^2);


The presence of the keyword function determines that this Matlab file is special: it is a function and not a simple script file.

Note the first line of the file has the form:

function output-variable = function_name(input-variable)


We make several observations:

1. The keyword function dicatates that this file is a function file, not a script.
2. The function name is f.
3. The input variable is t.
4. The output variable is y.

Because the function name is f, the rules of Matlab require that we save this file as f.m. That is, we must take the function name, append .m, then save. If you click the Save icon on the toolbar of the editor (or select File->Save), you can save the file in the usual manner. By default, Matlab suggests the file should be saved in the current directory (ArcLength) with the name f.m, but it is your responsibility to make sure that you are saving the file as f.m in the ArcLength directory.

Return to the command window. Again, make sure you are in the ArcLength directory with the command pwd. If not, change directories.

Because

f(t)=sqrt(4+5 cos^2 t,

then

f(0)=sqrt(4+5 cos^2(0))=sqrt(4+5(1)^2)=sqrt(9)=3.

We will now test our function.

>> f(0)
ans =
3


#### Trouble-Shooting

It's possible that you might arrive at a different result. There could be a number of reasons for this:

• Check your code in the editor. Did you type in the function correctly?
• Did you save the file as f.m in the directory ArcLength?
• Is your current directory ArcLength? Check this with the pwd command.

If all is still not well, there are two Matlab commands that can help.

• The command which f will give the path of the function f that Matlab will execute. If this path points to an f.m that is in a directory other than ArcLength, then check that you saved the file f.m to the ArcLength folder and your current directory is ArcLength (use the pwd command.)
• The command type f will type the contents of the file f.m to the Matlab command window. If it types a function that is completely different from yours, then you know that Matlab is finding a file f.m on its path in another location. Again, check that you saved the file in the ArcLen folder and your current directory points to the ArcLength Folder.

#### Making Your Function Array Smart

If you've made it to this point in the activity, then you know your function returns the correct response if you enter a single value. However, just to make sure, try the following code:

>> t=0
t =
0
>> f(t)
ans =
3


This is the same result as above. If this doesn't work, return to the Trouble-Shooting section and make the appropriate correction. Before continuing, this example must work.

Now, how will our function perform on a vector of values? Create a vector of t-values:

>> t=0:pi/2:2*pi
t =
0    1.5708    3.1416    4.7124    6.2832


The ideal would be that our function would be applied to each entry of this vector. Let's see:

>> f(t)
??? Error using ==> mpower
Matrix must be square.

Error in ==> f at 2
y=sqrt(4+5*cos(t)^2);


Let's analyze the error message:

1. First, the input t is a vector.
2. Matlab's cosine function is "array smart" so cos(t) is a vector, created by taking the cosine of every entry of the vector t.
3. Thus, cos(t) is a vector. The error is trying to raise a vector to the second power with the command cos(t)^2. This is an illegal operation that spawned the error message shown above.

We need to square every entry of the vector cos(t). To do this we use "dot notation." Make the following change to the file f.m and save.

function y=f(t)
y=sqrt(4+5*cos(t).^2);


>> f(t)
ans =
3     2     3     2     3


Aha! Our function is now "array smart." It evaluated the function at each entry of the vector t. The careful reader will use pencil and paper to evaluate the function at 0, pi//2, pi, 3pi//2, and 2pi to verify this result.

It is now a simple matter to approximate the integral

$L=\int_0^{2\pi} \sqrt{4+5\cos^2 t}$.

Simply enter the following at the Matlab command prompt:

>> quad(@f,0,2*pi)
ans =
15.8654


Thus,

$L=\int_0^{2\pi} \sqrt{4+5\cos^2 t}\approx 15.8654$.

If you type help quad at the Matlab prompt, a description of Matlab's quad command results, the first paragraph of which is:

 QUAD   Numerically evaluate integral, adaptive Simpson quadrature.
Q = QUAD(FUN,A,B) tries to approximate the integral of scalar-valued
function FUN from A to B to within an error of 1.e-6 using recursive
Y=FUN(X) should accept a vector argument X and return a vector result
Y, the integrand evaluated at each element of X. 

Students of calculus will focus on the phrase "tries to approximate the integral of scalar-valued function FUN from A to B to within an error of 1.e-6 using recursive adaptive Simpson quadrature." Surely this means that the quad command is using a sophisticated adaptation of Simpson's method that provides a approximation that is within 1 times 10^{-6} of the correct answer.

Here are a few more observations regarding the help file's description of use:

• FUN is a function handle. In our call, quad(@f,0,2*pi), the "at" symbol @ is used to create a function handle "on-the-fly". We'll have more to say about function handles in future activities. For now, simply prefix the @ symbol to the name of your function.
• In the description Q = QUAD(FUN,A,B), A and B are the lower and upper bounds of the integral. Thus, we passed the lower and upper bounds of our integral in the command quad(@f,0,2*pi).

### Some Comments on Writing Functions

We share some final thoughts on writing functions. First, you can use whatever name you want for your function. Clearly, you won't always want to use f as your function name. For one thing, you can have only one file at a time named f.m in your ArcLength folder. Similarly, you can use whatever names you wish for your input and output variables. With these thoughts in mind, create a new file in the editor and enter the following code:

function stink = skunk(rattled)
stink=sqrt(4+5*cos(rattled).^2);


Some observations:

1. The name of the function is skunk. Therefore, the file should be saved as skunk.m in the ArcLength directory.
2. The input variable is named rattled instead of t in this function M-file.
3. The output variable is named stink instead of y in this function M-file.

Note again that we made the function array smart. We can test this with

>> t=0:pi/2:2*pi
t =
0    1.5708    3.1416    4.7124    6.2832


and

>> skunk(t)
ans =
3     2     3     2     3


This is identical to our previous result. Note that the contents of the vector t in the command workspace is passed into the contents of the input variable rattled in the function workspace, so these names do not have to match.

We can pass a function handle to the skunk function to the quad command in a similar manner, by prefixing the function name with the at @ symbol.

>> quad(@skunk,0,2*pi)
ans =
15.8654


Same result!

### Exercises

Consider the following parametric representation:

$\begin{eqnarray} x&=&t^2\\ y&=&t^3, \end{eqnarray}$

defined on the interval 0 le t le 1. Perform each of the following tasks.

1. Sketch the parametric graph of the parametric function. Turn on the grid with grid on. Provide axis labels and an appropriate title, then obtain a printout.
2. On the printout of your plot, devise a strategy for estimating the length of the curve. You might try drawing a few line segments then using either the distance formula or the Pythagorean theorem to obtain an estimate of their total length.
3. Set up the integral on the printout of your plot for determining the length of the arc.
4. Write a function M-file for the integrand and obtain a printout of the file.
5. Use the quad command to approximate the integral in part (3). Place this command and the rsulting approximation on the printout of your plot and compare with the estimate in part (2)
6. Turn in the printout of your plot and the printout of your function M-file.