{VERSION 6 0 "IBM INTEL NT" "6.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 "" 1 14 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{PSTYLE "Normal " -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 31 "Moment of inertia calcula tions." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 161 "Now we extend the calculations from page2p3.mws to perform not ju st area (volume) integrals, but compute an integral over a function ev aluated inside the volume." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 187 "As an example to illustrate the idea we calculate the inertia of a cone defined by rad ius r0 and height h. We begin with the volume calculation, and then ex tend to integrate the function." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "Digits:= 15;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'DigitsG\"#:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "r0:=2;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#r0G\"\"#" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "h:=10; " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"hG\"#5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 96 "We work in Cartesian coordinates (for simplicity), i e. we don't use the symmetry of the problem." }}{PARA 0 "" 0 "" {TEXT -1 107 "The reference volume in this case is conveniently chosen as a \+ block of height h and base area 2r0 times 2r0" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "V_S:=(2*r0)^2*h;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$V_SG\"$g\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 174 "We start \+ with random numbers that are uniform in [0,1], but we will need to sca le x to cover -r0 <= x <= r0, and likewise y, while the z variable is \+ to be scaled from 0 to h." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "xymap:=xi->(2*xi-1)*r0;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%&xyma pGj+6#%#xiG6\"6$%)operatorG%&arrowGF(*&,&*&\"\"#\"\"\"9$F0F0F0!\"\"F0% #r0GF0F(F(F(6#\"\"!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "xymap (0);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#!\"#" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "xymap(1);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\" #" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "xymap(0.5);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"\"!F$" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "zmap:=xi->h*xi;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>% %zmapGj+6#%#xiG6\"6$%)operatorG%&arrowGF(*&%\"hG\"\"\"9$F.F(F(F(6#\"\" !" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 30 "Generate the lists of random s:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 107 "S1:=[stats[random, u niform](1000)]: S2:=[stats[random, uniform](1000)]: S3:=[stats[random, uniform](1000)]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "S1[1], S2[1],S3[1];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%$\"-\"3p'>uU!#7$\"-#*e /=CUF%$\",oiabW#F%" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 144 "We need th e criterion which will determine whether the point lies inside the con e or not. This requires some thinking about the geometric shape." }} {PARA 0 "" 0 "" {TEXT -1 158 "For a given height z the requirement is \+ that the distance from the z-axis not exceed the size of a circle of r adius r(z). What is the size of the circle r(z)?" }}{PARA 0 "" 0 "" {TEXT -1 81 "It is given by a line which passes through the points (z =0, r=r0), and (z=h,r=0)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "rmap:=z->r0-r0/h*z;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%rmapGj+6#% \"zG6\"6$%)operatorG%&arrowGF(,&%#r0G\"\"\"*(F-F.%\"hG!\"\"9$F.F1F(F(F (6#\"\"!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "rmap(0);" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"#" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "rmap(h);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"!" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "rmap(h/2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"\"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 24 "We mo dify our procedure:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 91 "MC:= proc(S1,S2,S3) local i,N,x_nu,y_nu,z_nu,r_nu,Nhit,i0,nu,Nhit_a,Nhit_d; global r0,h,V_S;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 122 "N:=nops(S1): \+ if N<>nops(S2) or N<>nops(S3) then RETURN(\"Problem with mismatched po int sequences\",N,nops(S2),nops(S3)) fi: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Nhit:=Vector(10,datatype=float[8]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "for i0 from 1 to 10 do: for i from 1 to N/10 do: nu:=(i0-1)*(N/10)+i;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 79 "x_nu:=xyma p(S1[nu]); y_nu:=xymap(S2[nu]); z_nu:=zmap(S3[nu]); r_nu:=rmap(z_nu); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "if x_nu^2+y_nu^2 <= r_nu^2 then Nhit[i0]:=Nhit[i0]+1; fi: od: od: Nhit:=Nhit/(N/10);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 144 "Nhit_a:=evalf(add(Nhit[i0],i0=1..10)/10); Nhit_ d:=sqrt(add((Nhit[i0]-Nhit_a)^2,i0=1..10)/(9*10)); V_S*[Nhit_a-Nhit_d, Nhit_a,Nhit_a+Nhit_d]; end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 76 "No te: on the last line we multiply the success rate by the reference vol ume!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "Volume:=MC(S1,S2,S3 );" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'VolumeG7%$\"0LLLLLt&R!#8$\"0+ ++++?6%F($\"0nmmmmmE%F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 " V_cone:=evalf(1/3*Pi*r0^2*h);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'V_ coneG$\"0Q'y/-z)=%!#8" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 179 "Now let us modify the procedure to calculate the inertia about the z-axis. Th is is not hard, because all we need to do is to calculate the average \+ of the distance from the z-axis." }}{PARA 0 "" 0 "" {TEXT -1 109 "Inst ead of updating a counter by one if we are inside the cone, we store t he average distance from the z-axis" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "MC_MI:=proc(S1,S2,S3) local i,N,x_nu,y_nu,z_nu,r_nu,N hit,i0,nu,Nhit_a,Nhit_d; global r0,h,V_S;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 122 "N:=nops(S1): if N<>nops(S2) or N<>nops(S3) then RETU RN(\"Problem with mismatched point sequences\",N,nops(S2),nops(S3)) fi : " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Nhit:=Vector(10,datatype=floa t[8]):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "for i0 from 1 to 10 do: f or i from 1 to N/10 do: nu:=(i0-1)*(N/10)+i;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 79 "x_nu:=xymap(S1[nu]); y_nu:=xymap(S2[nu]); z_nu:=zmap( S3[nu]); r_nu:=rmap(z_nu);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "if x_ nu^2+y_nu^2 <= r_nu^2 then Nhit[i0]:=Nhit[i0]+x_nu^2+y_nu^2; fi: od: o d: Nhit:=Nhit/(N/10);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 140 "Nhit_a:=e valf(add(Nhit[i0],i0=1..10)/10); Nhit_d:=sqrt(add((Nhit[i0]-Nhit_a)^2, i0=1..10)/(9*10)); [Nhit_a-Nhit_d,Nhit_a,Nhit_a+Nhit_d]; end:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "r2_average:=MC_MI(S1,S2,S3); " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%+r2_averageG7%$\"0xz^s>k\"G!#:$ \"0iljZH$yHF($\"0Z^vAR-9$F(" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "If we leave out the mass of the body (set it to one) we express the iner tia as" }}{PARA 0 "" 0 "" {TEXT -1 34 "M**V_S/Volume, using M =1:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "r2_average[2]*V_S/Vo lume[2];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"0qLZ1$))e6!#9" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 48 "This is to be compared with the te xtbook result:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 28 "MI_cone:=evalf(M*3/10*r0^2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(MI_coneG,$*&$\"0++++++?\"!#9\"\"\"%\"MGF*F*" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "What are the uncertainty brackets \+ of our calculation?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 262 "We can use the worst-case scenario: combine the upper \+ limit (standard deviation) of the numerator with the lower limit of th e denominator, and vice versa. This yields a conservative standard err or for the estimation of the MI about the symmetry axis of the cone." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "r2_average[3]*V_S/Volume[ 1];" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"0X0<]Q'p7!#9" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "r2_average[1]*V_S/Volume[3];" }} {PARA 11 "" 1 "" {XPPMATH 20 "6#$\"0TUpRdh0\"!#9" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}}{MARK "0 0 0" 31 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 0 1 2 33 1 1 }