201MS Resit Coursework (2016–17) Page 1 of 22 RESIT COURSEWORK This is an individual coursework. Estimated time: 20–25 hours. Percentage of coursework component: 100 Intended Learning Outcomes: 1. Apply mathematical techniques and modern mathematical software to engineering problems; 2. Apply probability distributions and linear regression to engineering problems. 3. use techniques from multivariable calculus; 4. apply mathematical techniques and modern mathematical software to engineering problems; 5. find Fourier series in trigonometric form and apply this to the solution of partial differential equations; Hand Out: 8th May 2017 Hand in: 3rd July 2017, 11:55pm (online) Poisson Distribution: 5 marks Normal Distribution: 8 marks Linear Regression: 6 marks Manufacturing Question: 11 marks The wave equation: 20 marks MATLAB simulation: 10 marks Maxima and minima: 25 marks Total: 85 marks REMEMBER: coursework hand-in now has a hard deadline. Be late, get zero. Plan to hand in at least a day early. Better still, a week. For full marks, all working must be shown. All program files should be submitted on paper as part of your working.201MS Resit Coursework (2016–17) Page 2 of 22 Part 1: Probability Theory 1. Poisson Distribution It is assumed that the number of people telephoning a certain help line is a Poisson process. The help line was open for calls between 0900 and 1700, Monday – Friday and, over a 1 week period, the number of calls per hour were logged. The result was: Mon. 3 0 1 2 0 2 1 0 Tue. 1 1 0 3 4 1 2 1 Wed. 2 1 3 5 1 0 1 3 Thu. 1 2 1 0 1 2 3 0 Fri. 1 0 4 1 0 2 1 2 Use this data to determine the relative frequency for n calls per hour, with n = 0, 1, 2, 3, 4 and 5. Find the mean number of calls per hour and hence estimate the above relative frequencies assuming a Poisson distribution. (5 marks)201MS Resit Coursework (2016–17) Page 3 of 22 Solution: The relative frequencies in the data are: (fill in missing values in row 2 of table) (2 marks) r 0 1 2 3 4 5 f The mean number of calls per hour is µ = (fill in value). (1 mark) If the process is Poisson then P(X = r) = µre−µ r! . For r = 0, 1, 2, 3, 4, 5 this gives (rounding to 4 dp) the following probabilities: (fill in missing values in row 2 of table) (2 marks) r 0 1 2 3 4 5 p 2. Normal Distribution A certain type of battery has a mean shelf life of 2 1 2 years with a standard deviation of 3 months. Assuming a normal distribution, estimate the probability that a battery chosen at random a) will have a shelf life of at least 3 years; b) will break down during the 3rd year; Sketch the corresponding areas under the standard normal curve for both cases. (8 marks)201MS Resit Coursework (2016–17) Page 4 of 22 Solution: Let X represent the random variable measuring the shelf life of a battery, assumed to be continuous. Using time units of months, the mean shelf life is given by µ = (fill in value), with the standard deviation of σ = (fill in value). (1 mark) It follows that Z = will be Standard Normal (fill in expression in X). (1 mark) We then have a) P X ≥ ! = P Z ≥ ! ≈ = % (2 marks) Sketch (1 mark) b) P ≤ X ≤ ! = P ≤ Z ≤ ! ≈ = % (2 marks) Sketch (1 mark) (fill in missing values and sketch the areas under the standard normal curve).201MS Resit Coursework (2016–17) Page 5 of 22 3. Linear Regression The number of miles per gallon in city driving is assumed to depend affinely on its weight, i.e. the fuel efficiency, c, in miles per gallon (m.p.g.), has the form c = a + bw , where w is the weight in kg. The following table lists c and w for 10 cars. (6 marks) w kg. 1955 1225 1008 1323 710 1350 1436 1561 2120 2110 c mpg 12.7 30.1 33.9 21.3 41.6 25.5 21.1 22.6 13.3 13.6 In MATLAB do the following: a) Enter the weights into a vector, w, and the corresponding fuel efficiencies in a vector c. b) Use polyfit to find the (degree one) polynomial p = a + bw that best fits the data in the least squares sense. c) Let x be a linearly spaced vector of weights going from w = 500 to w = 2500 containing 200 values (use the command x = linspace(500,2500,200); ) and calculate the corresponding y-values, where y = a + bx, using the polyval command. d) Produce a plot of y = a + bx together with the data using plot(x,y,w,c,’o’) .201MS Resit Coursework (2016–17) Page 6 of 22 Solution: a) The MATLAB commands used to enter w and c are: (enter commands in box) (1 mark) b) The Matlab command to find best linear fit is: (enter command in box) (1 mark) That is to say that the best linear fit is given by y = with coefficients given to 4 dp (enter expression for best straight line fit). (2 marks) c) The Matlab commands to evaluate this line at the given x values are: (enter the missing polyval command arguments) (1 mark) x = linspace(500,2500,200); y = polyval( , ); d) The Matlab commands: >> plot(x,y,w,c,’o’) >> title(num2str(clock)) produced the following graph: (attach a printout of the required graph) (1 mark)201MS Resit Coursework (2016–17) Page 7 of 22 4. Manufacturing Question A company receives an order to produce 150 customized Widgets of which 2% on average happen to be faulty. a) What is the chance of fulfilling the order by producing only the required number of Widgets ? (1 mark) b) How many Widgets are expected to be good on average in this case ? (1 mark) c) How many Widgets need to be produced to have on average the required number of good Widgets ? (2 marks) d) What is the chance of fulfilling the order by producing this number of Widgets ? (2 marks)201MS Resit Coursework (2016–17) Page 8 of 22 e) How many Widgets should be produced to fulfil the order with less than 1% risk ? (5 marks) Write your answers and all working in the boxes above. Part 1 30 marks201MS Resit Coursework (2016–17) Page 9 of 22 Part 2: The Wave Equation The problem is concerned with a uniform string stretched along the x-axis between x = 0 and x = L, where L is determined by the last digit of student ID, lSID, as L = 1 2 l lSID SID < ≥ 5, 5. Thus, lSID = and L = (fill in boxes). For the rest of the assignment use this value for L rather than variable L. The one-dimensional wave equation is given by 1 2c ∂2u ∂t2 = ∂2u ∂x2 , 0 ≤ x ≤ L, t ≥ 0, (1) where u(x, t) represents the displacement of the string at position x and time t. It is assumed that the string is held fixed at both ends and so we have fill in boxes (1 mark) u( , t) = u( , t) = 0, t ≥ 0. (2) Now assume the string is raised a distance µ at x = b and then released from rest at t = 0. Zero initial velocity implies that fill in boxes (1 mark) ∂u ∂t ( , ) = 0, 0 ≤ x ≤ L. (3) The initial displacement (plotted in figure 1) is given by fill in formula for b ≤ x ≤ L (2 marks) u(x, 0) = u0(x) =  µx/b 0 ≤ x ≤ b b ≤ x ≤ L. (4) We begin the assignment with some hand-calculations. Most of the calculations are provided, you need only fill in the missing details. Space is provided at the end for this. If you wish to word process your solutions, simply create your own template so that it is clear which part of your solution refers to which part of the assignment. Step 1: Use the method of Separation of Variables to derive solutions to (1)–(3) having the form u(x, t) = ∞∑ n=1 An sin nπx L cos nπct L , (5)201MS Resit Coursework (2016–17) Page 10 of 22 0 b L µ 0 Figure 1: Initial displacement of string. where An, n = 1, 2, . . . are constants that need to be determined in order that u(x, t) also satisfies the constraint of (4). (Note: You may assume that the constant value referred to by k in lecture notes or by α in the handout ‘PDE’s I: The Wave Equation’ is negative.) Solution: We begin by looking for solutions to Equations (1)–(3), having the form u(x, t) = X(x)T(t). (6) Substituting this into the wave equation gives put result into box below (1 mark) Dividing through by u(x, t) = XT then gives 1 2c T′′ T = X′′ X . (7) The left-hand side is now a function of t only and the right-hand side is a function of x only. It follows that both sides must equal the same constant value k say. Let us first analyse the equation X′′ X = k ⇐⇒ X′′ − kX = 0. (8) We only consider the case where k < 0. Let us first note that Equation (2) implies that X(0)T(t) = X(L)T(t) = 0. Since we are not interested in the trivial solution where T(t) ≡ 0, we deduce that fill in boxes (1 mark) X(0) = X(L) = .201MS Resit Coursework (2016–17) Page 11 of 22 Since k < 0, we set k = −p2, : Equation (8) has the general solution put general solution into box below (2 marks) for arbitrary constants A and B. The boundary conditions X(0) = X(L) = 0 give rise to fill in working in box below for nontrivial solutions of X(x) (4 marks) for any constant (fill in A or B as appropriate) and any positive integer n. Returning to (7), we also have 1 2c T′′ T = k = −p2 = − nL π 2 , where L is as setup at the beginning of this part. Rearranging this produces T′′ + nπ L c2 T = 0. (9)201MS Resit Coursework (2016–17) Page 12 of 22 The general solution to this is given by T(t) = C cos nπct L + D sin nπ Lct where C and D are arbitrary constants. However, Equation (3) leads to D = 0. show working for this in the box below (3 marks) Therefore T(t) = C cos nπct L . for any constant C and any positive integer n. Combining our expression for X(x) and T(t), we arrive at the following solutions to equations (1)–(3) from the assignment: u(x, t) = BC sin nπx L cos nπct L . (10) Since the equation to be solved is linear we can add such solutions together to arrive at the desired expression (Equation (5)): u(x, t) = ∞∑ n=1 An sin nπx L cos nπct L . (11) where An, n = 1, 2, . . . are arbitrary constants until we add the constraint of Equation (4). It follows from (11) that (1 mark) u(x, 0) = u0(x) = (12) Assuming An can be found to satisfy (12), we must have An = 2 L Z0L u0(x) sin nπ Lx dx (13)201MS Resit Coursework (2016–17) Page 13 of 22 Step 2: Use this formula (in Matlab or Maple or by hand), together with the expression for u0(x), as given by (4), to obtain (4 marks) An = (14) If using Matlab, you will need to set up µ, b, x and n as symbolic variables. Just use mu to represent µ. Split the integral of (13) into two parts, given the change of formula for u0(x) at x = b, as given in (4). After adding the two parts together, simplify the resulting formula using the simplify command. In Maple, it is possible to specify that n is integer valued, but I don’t know if you can do this in Matlab. You will therefore need to simplify sin(nπ) as zero yourself. Putting this expression for An into (5) gives u(x, t) = 2µL2 π2b(L − b) ∞∑ n=1 1 2n sin nπb L sin nπ Lx cos nπ Lct (15) Step 3: We finish the assignment by constructing an animation showing the vibration of the string. Animation in Matlab and Maple was introduced in Question 3 in the ‘Introduction to Fourier Series’ handout. If you are using Matlab, then you can fill in the missing details in the following m-file. (4 marks) % m-file for simulating one-dimensional wave equation. % initialise constants L= ; mu = ; N = ; b = ; c = 100; % set up x-vector for plots x = Llinspace(0,1,500); % set up limits for x and y axes v = [0 L -1.1mu 1.1mu]; % set up value for time step T % T is the period of oscillation divided by 15201MS Resit Coursework (2016–17) Page 14 of 22 T = ; % set up for loop for incrementing t in 16 % equal steps of length T for m = 1:16 t = (m-1)T; % initialise u to be zero then sum first N terms u = 0; for n = 1:N u = u+2muL^2/(pi^2b(L-b)n^2) sin(npib/L)sin(npix/L)cos(npict/L); end % plot this sum and store for movie using getframe figure(1),plot(x,u),axis(v) M(m) = getframe; % replot in a 4X4 array for displaying frames figure(2),subplot(4,4,m),plot(x,u),axis(v) end % show movie with 10 periods figure(1),movie(M,10) % display the individual frames figure(2) Attach your printouts of the Figure 2 window for three different numbers of the terms of Fourier series given by N = 5, 10, 50. Comment on the effect of N on the accuracy of solution by comparing the initial deflection of string (4) with solution (15) at t = 0 for different values of N. (6 marks) Use the following steps (or your own variation on them). No further guidance on Maple is given here, but you can ask if you are interested. 1. Assign a meaningful value for µ. Do not copy a friend’s values, use your own!201MS Resit Coursework (2016–17) Page 15 of 22 2. Choose your value of b as b = L(lSID + 1)/11 using the last digit of your student ID , lSID. 3. Create a suitable vector of x-values for simulating the vibrating string. 4. Create a vector to be used in an axis command to specify the limits on the x and y axes for the simulation. 5. Set up a FOR loop, for example using the command >> for m = 1:16 The idea is that each iterate of the loop will plot and store the string profile for a specific time instant. These time instants should start at t = 0 and be equally spaced, say T seconds apart. Assuming you are going to produce an animation consisting of 16 time instants, the last plot will be for t = 15T. You should then pick T such that 15T is one complete period. By this means, looping the animation will produce a sensible animation for as many periods as you desire. You should calculate what the period is, which will depend on your chosen value for c. You may wish to assign this value for T in Matlab, before the FOR loop (since it doesn’t change, putting such a command inside a loop is poor programming). 6. Define the time instant. For example using t = (m-1)T; where T has been constructed as suggested earlier. 7. Set up another FOR loop (inside the current one). The idea here is that this loop will construct the sum of the first N terms of the Fourier sum given by (15) at the time instant t. This is similar in nature to that given in the Fourier series handout and so no further details are given. 8. Plot this truncated sum, using prescribed axes limits, into Figure 1 and store the result using getframe. Also store the result into an array of plots in Figure 2. This then ends the outer FOR loop. 9. Use the movie command to display the animation for 10 periods. Part 2 30 marks201MS Resit Coursework (2016–17) Page 16 of 22 Part 3: Multivariate Calculus A tent has volume V which is determined by the last digit of student ID, lSID, as V = =  4 m3 if lSID = < 5, 6 m3 if lSID = ≥ 5. (fill in boxes). For the rest of the assignment, use this value of V where required. The tent has the shape shown in the figure with ends but no floor. It is desired to make the tent with the least possible material. ✂✂✂✂✂✂✂ ❆ ✏ ❆❆ ✏ ❆❆ ✏ ❆❆ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏ ✏✏ ❆❆✏ ❆❆✏ ❆✏ ❆❆ 2w θ L θ 1. Derive a formula for the volume, V, in terms of w, L and θ. Then rewrite this formula with L as the subject. Solution: The area of each triangular end is given by Ae = wh, where the tent height is given by h = . (1 mark) It follows that the volume is given by V = . (1 mark) Rewriting with L as the subject gives (1 mark) L = . (16) 2. Derive a formula for the area of material, A, in terms of w, L and θ and then substitute for L as given by Eq. (16).201MS Resit Coursework (2016–17) Page 17 of 22 Solution: The area of each end is Ae and the area of each side is As = LH, where H is the length of each upper side of the triangular ends, i.e., H = . (1 mark) The total area is therefore given by A = 2Ae + 2As, which simplifies down to A = 2w2 tan θ + 2V w csc θ (where csc θ = (sin θ)−1). Show the working for this in the box below (3 marks) Given that V is a constant, this gives A as a function of the 2 variables, w and θ. 3. Using either hand-calculation or Matlab (with Symbolic Toolbox), find the two partial derivatives, ∂A ∂w and ∂ ∂θ A. (Aside: You can use Maple if you prefer, but you may have trouble using Maple for part (d)). You will find Matlab commands to answer question 5 from the Maxima and Minima handout (but without any explanatory comments) at the end of the assignment. You can modify these appropriately.)201MS Resit Coursework (2016–17) Page 18 of 22 Solution: Insert Matlab commands or hand calculations into the box below with the final results in the following boxes (4 marks) This gives ∂A ∂w = ∂A ∂θ = 4. Find the stationary point (subject to the constraints that w > 0 and 0 < θ < π 2 ). If using Matlab, use the solve command. By default, Matlab will also find stationary points that do not obey the above 2 constraints (even complex solutions). You will have to pick out the relevant one. Substitute the values w = w0, θ = θ0 say, as given by the stationary point for the volume of the tent specified above, into the formula for L in part (a) to find the corresponding value L = L0.201MS Resit Coursework (2016–17) Page 19 of 22 Solution: Attach Matlab commands and/or hand-calculations. (6 marks) θ0 = w0 = L0 = 5. Finally, use the second derivatives test (or otherwise) to show that the stationary point of part (d) is a local minimum. Again, you may use Matlab (or Maple) if you wish. Solution: First we calculate the second partial derivatives evaluated at the stationary point . Put the Matlab commands or the hand calculations for ∂2A ∂w2 , ∂2A ∂θ2 , and ∂∂w 2A ∂θ at (w0, θ0) into the box below. Solution: Put explanation in the box below (6 marks) Explain how the above can be used to show the stationary point is indeed a local minimum.201MS Resit Coursework (2016–17) Page 20 of 22 Solution: Put explanation in the box below (2 marks) Part 3 25 marks Total 85 marks201MS Resit Coursework (2016–17) Page 21 of 22 Possible Matlab commands for Part 3. Please note that the order of solutions returned by the Matlab ‘solve’ command matches the order of the unknowns specified in the command. The order of solutions is not affected by the names of variable on the left hand side. For example, the command ‘[a b]= solve(eq1,eq2,x,y)’ returns solution of two simultaneous equations ‘eq1’ and ‘eq2’ first for the unknown ‘x’ and then for the unknown ‘y’, which are stored as ‘a’ and ‘b’, respectively. >> syms w d V >> A = 3wd + 2V/w + 2V/d; >> Aw = diff(A,w); pretty(Aw) V 3 d - 2 ---- 2 w >> Ad = diff(A,d); pretty(Ad) V 3 w - 2 ---- 2 d >> [wsol dsol] = solve(Aw,Ad,w,d) wsol = [ 1/318^(1/3)V^(1/3)] [ -1/618^(1/3)V^(1/3)+1/6i3^(1/2)18^(1/3)V^(1/3)] [ -1/618^(1/3)V^(1/3)-1/6i3^(1/2)18^(1/3)V^(1/3)] dsol = [ 1/318^(1/3)V^(1/3)] [ -1/618^(1/3)V^(1/3)+1/6i3^(1/2)18^(1/3)V^(1/3)] [ -1/618^(1/3)V^(1/3)-1/6i3^(1/2)18^(1/3)V^(1/3)] >> w0 = wsol(1); pretty(w0)201MS Resit Coursework (2016–17) Page 22 of 22 1/3 1/3 1/3 18 V >> d0 = dsol(1); pretty(d0) 1/3 1/3 1/3 18 V >> h0 = V/(d0w0); pretty(h0) 1/3 1/3 1/2 18 V >> Aww = subs(diff(Aw,w),{w,d},{w0,d0}) Aww = 6 >> Add = subs(diff(Ad,d),{w,d},{w0,d0}) Add = 6 >> Adw = subs(diff(Ad,w),{w,d},{w0,d0}) Adw = 3 >> Delta = AwwAdd - Adw^2 Delta = 27