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.

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:

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

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);

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:

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.