## 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.

- You can select
**File->New->M-file**from Matlab's menu. - You can click the
**New M-file**icon on Matlab's toobar. - 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:

- The keyword
**function**dicatates that this file is a function file, not a script. - The function name is
**f**. - The input variable is
**t**. - 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:

- First, the input
**t**is a vector. - 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**. - 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);

Return to the Matlab command window and try again.

>> 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.

### Using Matlab's QUAD Command

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 adaptive Simpson quadrature. FUN is a function handle. The function 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:

- The name of the function is
**skunk**. Therefore, the file should be saved as**skunk.m**in the ArcLength directory. - The input variable is named
**rattled**instead of**t**in this function M-file. - 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.

- 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. - 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.
- Set up the integral on the printout of your plot for determining the length of the arc.
- Write a function M-file for the integrand and obtain a printout of the file.
- 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) - Turn in the printout of your plot and the printout of your function M-file.