WR Long MS325 TMA 03 Question 4 # Part (a) restart:with(plottools):with(plots): l1:=line([0, 0], [0, 12], color = black, linestyle = solid): l2:=line([0, 12], [12, 12], color = black, linestyle = solid): l3:=line([12, 12], [12, 0], color = black, linestyle = solid): l4:=line([12, 0], [0, 0], color = black, linestyle = solid): c1 := circle([3, 9], 1, color = black): c2 := circle([9, 9], 1, color = black): e1:= ellipse([6, 6], 1, 2, filled = false, color = black): e2:= ellipse([6, 2], 4, 1, filled = false, color = black): display({l1,l2,l3,l4,c1,c2,e1,e2},scaling=constrained,axes=none); # Part (b) #Exact mass is the mass of a 12x12 plate minus the mass of each hole removed a_c1:=Pi:a_c2:=Pi: # mass of circles removed a_e1:=Pi*2:a_e2:=Pi*4: # mass of ellipses removed m:=12^2-a_c1-a_c2-a_e1-a_e2;evalf(%); # Part (c) # Note: to avoid repeated code, define f as follows: # Function rh for the density x1(x, y), so that takes the value 1 for 0 \342\200\26044 x \342\200\26044 12 # and 0 \342\200\26044 y \342\200\26044 12 except on the holes,where it is 0. rh:= proc(x,y) local c1,c2,e1,e2,pp; pp:=1; c1:=evalf((x-3)^2+(y-9)^2-1); c2:=evalf((x-9)^2+(y-9)^2-1); e1:=evalf((x-6)^2+((y-6)^2)/(2^2)-1); e2:=evalf(((x-6)^2)/(4^2)+(y-2)^2-1); if c1<0 or c2<0 or e1<0 or e2<0 then pp:=0 end if; pp; end: # Define f in this case equal to x1(x,y). This is to avoid repeated code for new definitions of f # in later parts of the question. f:=(x,y)->rh(x,y): plot3d('f(x,y)',x=0..12,y=0..12,axes=boxed,orientation=[-140,50],grid=[200,200]); mc4:=proc(Lim,M,N) # verbatim From unit C2 p89 local z,zz,a,b,h,Nt,z1,z2,X,i,j,k,Je,sig2,sige,x,M1,M2,n,jj; uses Statistics; n:=2: # Parameter defining the integral dimension a:=[seq(Lim[k],k=1..n)]; # Lower integration limits b:=[seq(Lim[k],k=1..n)]; # Upper integration limits h:=[seq((b[k]-a[k])/M[k],k=1..n)]; Nt:=N*mul(M[k],k=1..n); # The total number of sample points for k from 1 to n do: X[k]:=RandomVariable(Uniform(0,1)); od: for k from 1 to n do: z[k]:=Sample(X[k],Nt): od: j:=1: # Initialise the random number counter Je:=0: #Initialise the estimate for the integral sig2:=0: #Initialise the estimate for the standard deviation for k from 1 to M do: # Loop round strata in direction 1 for k from 1 to M do: # Loop round strata in direction 2 M1:=0: M2:=0: # Initialise the moments. Now loop round N points of each stratum for i from 1 to N do: for jj from 1 to n do: # Loop round number of dimensions x[jj]:= a[jj] + (k[jj] - 1 + z[jj][j])*h[jj]; od: zz:=evalf( f(x,x) ); M1:=M1 + zz; M2:=M2 + zz*zz; j:=j + 1; od; End of loop for integration in strata Je:=Je + h*h*M1/N; # Equation 2.67 sige:=N*(h*h)^2*(M2/N -(M1/N)^2)/(N-1); # Equation 2.68 sig2:=sig2 + sige/N; # Equation 2.69 od; od; [Je,evalf(2*sqrt(sig2))]; end: # Define limits of integration Lim:=[[0,12],[0,12]]: m_est:=mc4(Lim,[16,16],10); m_est-evalf(m); #The result compares well to that found in part b), being within (in this case) 99.5% of the exact mass #Part (d) # Redefinte f as xx1(x,y), so that f takes the value x for 0 \342\200\26044 x \342\200\26044 12 # and 0 \342\200\26044 y \342\200\26044 12 except on the holes, where it is 0. f:=(x,y)->x*rh(x,y): #Compute estimate for x coordinate of centre of mass #X_est:=(1/m_est)*mc4(Lim,[16,16],128): X_int:=mc4(Lim,[12,12],8): #integral for the x coordinate X_est:=evalf((1/m)*X_int); # Redefine f as yx1(x,y), so that f takes the value y for 0 \342\200\26044 x \342\200\26044 12 # and 0 \342\200\26044 y \342\200\26044 12 except on the holes, where it is 0. f:=(x,y)->y*rh(x,y): Y_int:=mc4(Lim,[12,12],8): #integral for the y coordinate Y_est:=evalf((1/m)*Y_int); # The absolute errors need to be scaled by 1/m X_err:=evalf((1/m)*X_int); Y_err:=evalf((1/m)*Y_int); # which are within 0.1 as required. # Comment: First, the estimated value for the x coordinate is close to 6, which would be expected. By inspection of the plot of the shape of the plate it is expected that the y-coordinate for the centre of mass would be greater than 6, which is the case. The true value of X is 6 due to the symmetry of the plate. ie the holes are completely symmetrical about the line x=6, so it follows that the centre of mass of the plate must have an exact x-coordinate of 6. We can see that the estimate of X is indeed within 6+/-0.1 of the true value # Part (e) # Redefine f f:=(x,y)->((x-6)^2-(y-6)^2)*rh(x,y): mi:=mc4(Lim,[32,32],128): # So the estimate of the moment of inertia is mi_est:=mi; # The relative error is mi_err:=mi/mi_est; # which is within 0.01 as required.