<Text-field style="_cstyle301" layout="Heading 1"><Font bold="true">Means and Moments and Exploring New Plotting Techniques</Font></Text-field>

<Text-field style="_cstyle301" layout="Heading 1"><Font bold="true">Means and Moments and Exploring New Plotting Techniques</Font></Text-field>

Note: You may notice differences between this Maple worksheet and the equivalent Mathematica notebook. These differences were introduced to preserve the content of these modules and were necessary because of major functional differences between Maple and Mathematica.

<Text-field style="Heading 1" layout="Heading 1">Introduction</Text-field>

OBJECTIVE: Extend the concept of the moments of a density function of a solid to applications in probability and engineering.

What do means and multiple integrals have to do with probabilities? How can the method of moments be used to help determine whether or not an object will float in an upright position? These questions will be answered as you explore this module. You will also get to practice new and interesting ways to plot functions.

<Text-field style="_cstyle256" layout="Heading 1"><Font size="18">Technology Guidelines</Font></Text-field>

NOTE: If you have just finished a worksheet, restart Maple before executing a new worksheet. TO OPEN SECTIONS, Click on the Arrow sign at the left hand side of the screen or select Expand All Sections from the View drop down menu.

TO STOP AN EXECUTION Click on STOP button from the toolbar.

ORDER OF EXECUTION Execute commands in the order given. Do not skip any Maple Input lines within a given worksheet

Alternatively, you can execute the entire worksheet by selecting the Execute Worksheet command from the Edit drop down menu.

SAVING WORKSHEETS.

You can save anytime to any directory you choose, and it is wise to save often. EXPERIENCING MAJOR PROBLEMS Save if appropriate, and then sut down Maple and start it up again.

<Text-field style="Heading 1" layout="Heading 1">Part I: Masses and Moments Applied to Probability Functions</Text-field>

In multivariable probability, probability density functions are defined over specified regions (domains) in the same way that as density functions are defined over regions of a solid. These probability functions are designed so that the total "mass" is always 1. Probabilities are then defined to be the integrals of the density functions over a particular subregion of the domain. Consider the following example.

Suppose that a national fast-food outlet is interested in the joint behavior of the random variables x, defined as the total time between a customer's arrival at the store and departure from the service window, and y, the time that a customer waits in line before reaching the service window. Since x contains the time a customer waits in line, x must be greater than y. The relative frequency distribution of observed values of x and y can be modeled by the probability density function: LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2KEYrLUYjNiYtRiw2JFEiZkYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnL0Y3USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQS8lKXN0cmV0Y2h5R0ZBLyUqc3ltbWV0cmljR0ZBLyUobGFyZ2VvcEdGQS8lLm1vdmFibGVsaW1pdHNHRkEvJSdhY2NlbnRHRkEvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZQLUkobWZlbmNlZEdGJDYlLUYjNiYtRiw2JFEieEYnRjYtRjo2LlEiLEYnLyUnaXRhbGljR0ZBRj1GPy9GQ1EldHJ1ZUYnRkRGRkZIRkpGTEZOL0ZSUSwwLjMzMzMzMzNlbUYnLUYsNiRRInlGJ0Y2RjZGaG5GPUY2LUY6Ni5RIj1GJ0ZobkY9Rj9GQkZERkZGSEZKRkwvRk9RLDAuMjc3Nzc3OGVtRicvRlJGZW8tRiM2Ji1JI21uR0YkNiVRJTAuMDRGJ0ZobkY9LUY6Ni1RMSZJbnZpc2libGVUaW1lcztGJ0Y9Rj9GQkZERkZGSEZKRkxGTkZRLUklbXN1cEdGJDYlLUY6Ni5RLyZFeHBvbmVudGlhbEU7RidGaG5GPUY/RkJGREZGRkhGSkZMRk4vRlJRLDAuMTExMTExMWVtRictRiM2KC1GOjYuUSomdW1pbnVzMDtGJ0ZobkY9Rj9GQkZERkZGSEZKRkwvRk9RLDAuMjIyMjIyMmVtRicvRlJGXnEtRiM2Ji1Gam82JVEkMC4xRidGaG5GPUZdcEZYRjYtRjo2LlEoJm1pbnVzO0YnRmhuRj1GP0ZCRkRGRkZIRkpGTEZdcUZfcS1GIzYmLUZqbzYlUSQwLjNGJ0ZobkY9Rl1wRl5vRjZGK0Y2LyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJ0Y2RitGNkYrRjY= for LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2Ji1JKG1mZW5jZWRHRiQ2JS1GIzYmLUkjbW5HRiQ2JVEiMEYnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYuUSUmbGU7RidGOkY9LyUmZmVuY2VHRjwvJSpzZXBhcmF0b3JHRjwvJSlzdHJldGNoeUdGPC8lKnN5bW1ldHJpY0dGPC8lKGxhcmdlb3BHRjwvJS5tb3ZhYmxlbGltaXRzR0Y8LyUnYWNjZW50R0Y8LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGVC1GLDYkUSJ5RicvRj5RJ2l0YWxpY0YnRlpGOkY9RkAtRiw2JFEieEYnRlpGWkYrRlo= <LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JFEoJmluZmluO0YnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJ0Yv and 0 elsewhere, where x and y are measured in minutes.

restart; with(plots):

f:=0.04*exp(-.1*x-.3*y);

We examine the domain of this function. First read in the plots package to assist us in graphing.

plot(x, x=0..5, labels=[`time in system`, `time in line`], filled=true, color=grey);

We begin by verifying that the integral of the density function over the entire domain is 1.

round(int(int(f, y=0..x), x=0..infinity));

To compute the average or expected values of x and y over the domain, we find the first moments of the density function with respect to the x and y axes. The second integral is difficult to evaluate as x approaches LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JFEoJmluZmluO0YnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJ0Yv, so we substitute a large value of y for the upper bound.

time1:=int(int(x*f, y=0..x), x=0..infinity): print( `Expected time in system = `, time1, `minutes`); time2:=int(int(y*f, y=0..x), x=0..1000000): print( `Expected time waiting in line = `, time2, `minutes`);

Now we compute the probability that y is at least one-half as large as x. First we can draw the region over which we will integrate (the part for x < 5).

polygonplot([[0,0],[5,5],[5,5/2],[0,0]],color=COLOR(RGB,0,1,1));

The probability of this event is the integral of the density function over this region.

prob1:=evalf(int(int(f,y=.5*x..x),x=0..infinity)): print(`The probability that the time in line is at least half as long as the time in the system is `,prob1);

If we want to know the probability that a customer spends no more than six minutes at the service window itself, we represent that quantity by x - y, and we integrate over the portion of the domain where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2J0YrLUYjNiYtRiw2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RKCZtaW51cztGJy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yMjIyMjIyZW1GJy8lJ3JzcGFjZUdGUy1GLDYlUSJ5RidGNkY5RkAtRj02LVElJmxlO0YnRkBGQkZFRkdGSUZLRk1GTy9GUlEsMC4yNzc3Nzc4ZW1GJy9GVUZnbi1JI21uR0YkNiRRIzYuRidGQEZARitGQA==

polygonplot([[0,0],[20,20],[20,14],[6,0]],color=grey);

Looking at the domain, you can see that it is easier here to integrate with respect to x first.

prob2:=int(int(f,x=y..y+6),y=0..infinity):print(`The probability that the time spent at the service desk is no more than 6 minutes is`, prob2);

<Text-field style="Heading 1" layout="Heading 1">You Try It: Part I</Text-field>

Try your hand at computing the probability that the time at the service desk will be at least 8 minutes. Change the entries to the appropriate terms.

prob3:=int(int(f,x=y..y+6),y=0..infinity): print(`The probability that the time spent at the service desk is no more than 6 minutes is`, prob3);

Compute the probability that the time in line is no more than 30% of the time in the system. Change the appropriate terms.

prob4:=int(int(f,y=.5*x..x),x=0..infinity): print(`The probability that the time in line is half as long as the time in the system is`, prob4);

<Text-field style="Heading 1" layout="Heading 1">Part II: Using Moments to Determine Acceptable Heights</Text-field>

Designing a Buoy

Suppose that you have a buoy in the form of a solid hemisphere of radius 40 centimeters and uniform density that is topped with a solid cylinder of the same material and radius 40 centimeters. This buoy floats with the center of the hemisphere above the water surface. Determine the maximum height that the cylinder may attain that permits the buoy to continue floating upright. To picture this, load two packages, and plot the shape of the body.

hemisplot:=sphereplot(40, r=0..2*Pi, theta=0..Pi, labels=["x","y","z"], view = [-1..1,-1..1,-1..0], axes=boxed): pcyl:=plottools[cylinder]([0,0,1],40,40): flagpole:=plot3d([0,0,z], z=35..80,x=-10..10,axes=boxed, thickness=3): flag:=plot3d([0,u,v], u=0..20, v=60..80, axes=boxed): display({hemisplot, pcyl, flagpole, flag}, view=[-50..50, -50..50, -50..80], scaling=CONSTRAINED);

You may assume that the flag on top is weightless.

The method of moments can help us determine that the buoyancy force exerted by the body of water due to the amount of water displaced will be in a direction through the point that would be the midpoint of the sphere, as long as the object is upright or tilted so that no part of the cylinder is immersed. Hence, if the center of gravity of the solid is at that point or closer to to the bottom of the hemisphere, the object will stay upright. If the center of gravity is in the region of the cylinder, the object will tilt. We first determine the z-coordinate of the center of gravity of the hemisphere. Recall that the volume of a hemisphere with radius of r is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbWZyYWNHRiQ2KC1GIzYoLUkjbW5HRiQ2JFEiMkYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RMSZJbnZpc2libGVUaW1lcztGJ0Y0LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRj0vJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRkwtSSNtaUdGJDYlUSUmcGk7RicvJSdpdGFsaWNHRj1GNEY3LUklbXN1cEdGJDYlLUZQNiVRInJGJy9GVFEldHJ1ZUYnL0Y1USdpdGFsaWNGJy1GMTYkUSIzRidGNC8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGNC1GIzYkRmluRjQvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmZvLyUpYmV2ZWxsZWRHRj1GNA==.

volhemis:=2/3*Pi*40^3: hemis:=-sqrt(1600-r^2): unassign('r,z,theta'); zbar:=int(int(int(r*z,z=hemis..0), r=0..40),theta=0..2*Pi)/volhemis;

The z-coordinate of the center of mass of the cylinder is h/2, so if we equate the first moment of the hemisphere to the first moment of the cylinder, we can solve for h. Since the radius is 40 here, the volume of the cylinder is 1600*Pi* h.

height:=solve(-zbar*volhemis=h/2*(1600*Pi*h),h);

Therefore, the height can be at most 20 times the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbXNxcnRHRiQ2Iy1JI21uR0YkNiRRIjJGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGMg==. In general, if the radius were R, the maximum height of the cylinder could be at most R times the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkmbXNxcnRHRiQ2Iy1JI21uR0YkNiRRIjJGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGMg== / 2.

<Text-field style="Heading 1" layout="Heading 1">You Try It: Part II</Text-field>

Experiment with some of the 3-D plot functions presented in this section. First just practice with a cylindrical plot. In the following example, z is a function of r and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJnRoZXRhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGMg==. The first argument of cylinderplot( ) is a list of the three cylindrical coordinates of points on the surface, and we express the z-coordinate as a function of the first two, r and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJnRoZXRhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRidGMg==. Try some of your own functions. Change the required terms.

unassign('theta','r'); z:=r*(cos(theta))^2; cylinderplot([r,theta,z], r=0..10, theta=0..2*Pi,labels=["x","y","z"],axes=boxed);

Try it out with some other functions.

You can create spherical plots with the same package you used to create cylindrical plots. In these plots, it is assumed that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEmJnJobztGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnL0YzUSdpdGFsaWNGJw== is a function of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEmJnBoaTtGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnL0YzUSdpdGFsaWNGJw== and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEoJnRoZXRhO0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRicvRjNRJ2l0YWxpY0Yn. Try some of your own by changing the terms in the function LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUkjbWlHRiQ2JVEmJnJobztGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnL0YzUSdpdGFsaWNGJw==. Not everything will give satisfying results.

rho:=theta*cos(phi);

sphereplot(rho, theta=0..2*Pi, phi=0..Pi, labels=["x","y","z"],axes=boxed);

You may have noticed that to draw the flag we had to first generate points and then plot the points. The difference between the two was that the pole could be thought of as a line, whereas, the flag was more like a surface. Experiment with the following commands and draw something else by changing the terms in the plot commands. Remember to leave one parameter in the first expression and two parameters in the second experession.

unassign('z'); newflagpole:=plot3d([0,0,z], z=0..30,x=-10..10,axes=boxed, thickness=3): print(newflagpole); newflag:=plot3d([0,u,v], u=0..20, v=15..30, axes=boxed): print(newflag); print(display({newflag, newflagpole}));