{VERSION 5 0 "SUN SPARC SOLARIS" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 257 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Outpu t" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Plot" 0 13 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 256 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 257 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }3 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 11 258 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 259 1 {CSTYLE " " -1 -1 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "" 0 260 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {PARA 11 "" 1 "" {TEXT -1 0 "" }{TEXT 256 0 "" }{TEXT 257 11 " MATH 2270-2" }}{PARA 256 "" 0 "" {TEXT -1 30 "MAPLE PROJECT 2a - Power Laws " }}{PARA 257 "" 0 "" {TEXT -1 16 "October 22, 2001" }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT 259 40 "Using logarithms to discover power laws:" }}{PARA 0 " " 0 "" {TEXT -1 333 " Here's how to use a least squares method to \+ test whether a power law can adequately describe a correlation between two different variables. Suppose we have some data points, \{[xi,yi] \}, and we expect a power law yi=C*(xi)^p to relate the yi's to the x i's. Here p is the power, and we can call C the proportionality cons tant. " }{TEXT 258 206 "Then we expect the logarithms to satisfy ln(y i) = ln(C) + p*ln(xi). In other words, the ln-ln data should satisfy \+ the equation of a line, having slope equal to p and vertical axis inte rcept equal to ln(C)." }{TEXT -1 277 " So we can take the least squar es fit for the ln-ln data, and then use the slope and intercept to ded uce C and p in the power law: take the power p equal to the slope of t he least square fit to the ln-ln data, and take C to equal to the expo nential of its vertical intercept." }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 260 "" 0 "" {TEXT -1 8 "Example:" }}{PARA 0 "" 0 "" {TEXT -1 413 " We work exercise #40 in section 5.4 of the text. This examp le relates to Johannes Kepler's discovery several hundred years ago th at for circular planetary orbits the period of revolution is proportio nal to the radius or the orbit, to the 3/2 power. Later Newton showed that the full set of Kepler's observational Laws was consistent (only ) with the famous \"inverse-square law\" for gravitational attraction. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "with(lin alg):with(plots):\n #we use these two command libraries" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 150 "S:=[[.387,.241],[1,1],[5.20,11.86] ,\n [19.18,84.0],[39.53,248.5]]:\n#the actual data, in table f orm\nA1:=convert(S,matrix); #and as a matrix\n " }}}{PARA 0 "" 0 " " {TEXT -1 188 "Just for fun I have reproduced the table from the text . The units of distance are \"Astronomical Units\", i.e. the distance between the earth and the sun. The time units are (Earth) years." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 213 "planets:=matrix(5,1,[`Mercu ry`,`Earth`,`Jupiter`,\n `Uranus`,`Pluto`]):\nheaders:=matrix(1,3,[ `Planet`,`Mean radius`,`mean period`]):\nbooktable:=stackmatrix(header s,augment(planets,A1));\n #the table in the book" }}{PARA 11 "" 1 " " {XPPMATH 20 "6#>%*booktableG-%'matrixG6#7(7%%'PlanetG%,Mean~radiusG% ,mean~periodG7%%(MercuryG$\"$(Q!\"$$\"$T#F17%%&EarthG\"\"\"F67%%(Jupit erG$\"$?&!\"#$\"%'=\"F;7%%'UranusG$\"%=>F;$\"$S)!\"\"7%%&PlutoG$\"%`RF ;$\"%&[#FD" }}}{PARA 0 "" 0 "" {TEXT -1 106 "Now that we've had our f un, let's get down to the task of discovering the power law, using the A1 matrix." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 209 "A2:=map(ln, \+ A1): \n #the map command will apply a given function,\n #in this \+ case the natural logarithm, to each\n #entry of the specified matrix . So the entries of\n #A2 will be the ln of entries in A1." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 116 "A3:=map(evalf,A2); \n #g et each decimal value (not necessary in\n #this example, but needed for your work on BMI" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 182 "r owdim(A3); \n #rowdim computes the number of rows in\n #a matrix. Of course, for this small matrix,\n #we know there are five rows. \+ Your B.M.I. matrix\n #will be much bigger\n" }{TEXT -1 0 "" }}} {PARA 0 "" 0 "" {TEXT -1 361 "We use A2 to construct the matrix \"A\" \+ and right hand side b, for the least squares line fit. The first col umn of A will be the ln(x)-values, the second column will be a vector \+ of 1's and the vector b will be the corresponding ln(y) values, since \+ we want a least squares line fit (lny)=m(lnx) + b. Maple has commands to extract columns, augment matrices, etc." }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 451 "col2:=vector(rowdim(A3),1); \n #this will be th e second column, a vector of 1's\n #for our least-squares line fit m atrix, see e.g.\n #page 219 of the text. (Out of habit I make this \n #the second column of the matrix rather than the\n #first, and \+ then reverse m and b in the matrix equation.\n #Sorry.) This comman d creates a vector\n #with number of entries equal to the first argu ment\n #and makes each entry equal to the second argument\n" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "A4:=delcols(A3,2..2); \n # remove the second column of A3," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 133 "A:=augment(A4,col2); #This is our matrix for least squares, \n #with the ln(x) terms in the first column, and 1's in\n #the se cond." }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"AG-%'matrixG6#7'7$$!+geI$ \\*!#5\"\"\"7$\"\"!F-7$$\"+E'e'[;!\"*F-7$$\"+p!oQ&HF3F-7$$\"+x)fqn$F3F -" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 97 "b:=delcols(A3,1..1); \+ #This is the right-hand side for least squares,\n #it has the ln(y) \+ values" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"bG-%'matrixG6#7'7#$!+X$e HU\"!\"*7#\"\"!7#$\"+%RrJZ#F,7#$\"+*z;3V%F,7#$\"+YGW:bF," }}}{PARA 0 " " 0 "" {TEXT -1 75 "The next three commands find the least-squares sol ution, as in section 5.4." }}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 110 "ATA:=evalm(transpose(A)&*A);\nATb: =evalm(transpose(A)&*b);\nlinsolve(ATA,ATb);\n #solve the system (AT A)x=(ATb)" }}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{PARA 11 "" 1 "" {TEXT -1 0 "" }}}{PARA 258 "" 1 "" {TEXT -1 73 "Alternately, we could multip ly the right-hand side by the inverse matrix:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "evalm(inverse(ATA)&*ATb);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#-%'matrixG6#7$7#$\"+8k\")*\\\"!\"*7#$\"'!*o[F*" }}} {PARA 0 "" 0 "" {TEXT -1 162 "Actually, least squares for linear syste ms is a standard tool so Maple has the command built right in. You ca n read about it by searching the topic \"leastsqrs\".\n" }}{PARA 0 "" 0 "" {TEXT -1 144 "The first entry above is the least-squares line slo pe (for the ln-ln data), the second point its intercept. See visually how the line-fit went:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 174 "lnlnplot:=pointplot(\{seq([A3[i,1] ,A3[i,2]],i=1..rowdim(A3))\}):\n #those are the points from the ln-l n data in\n #the A3 matrix. the index i ranges over all rows\n #in A3.\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 174 "line:=plot(1.499 816412*t + .4868908465e-3, t=-1..4):\n #I used the mouse to paste in the coefficients\n #from my work above. Here t is standing for the \n #ln(x) variable." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "di splay(\{lnlnplot,line\});" }}{PARA 13 "" 1 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 325 "The line fit should be pretty good, because Ke pler wasn't a moron. (Actually, he didn't get to use Pluto in his dat a set. Pluto wasn't discovered until the 1900's.) The slope of this \+ line will be our experimental power, its intercept will be the ln of o ur proportionality constant. Now work backwards to get these values: " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "C:=exp(.486890e-3); #pr oportionality constant\np:=1.499816413; \n #power. I used the mous e to paste these in." }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"CG$\"+4q[+ 5!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"pG$\"+8k\")*\\\"!\"*" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 157 "Notice, the power came out close \+ to 1.5 and the proportionality constant came out close to 1. Now see \+ how our power law works for the original (xi,yi) data:" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 61 "realplot:=pointplot(\{seq([A1[i,1],A1[i,2]],i= 1..rowdim(S))\}):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "powerp lot:=plot(C*r^p,r=0..40,color=black):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "display(\{realplot,powerplot\});" }}{PARA 13 "" 1 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 259 "" 0 "" {TEXT -1 240 "Your job in the first part of this Maple project will be to carry out the same sort of analysis for the height-weight data whi ch we have accumulated. If you go to our Maple page you may download the template for this part of your project." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{MARK "20 0" 56 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }