AMATH 301 – Spring 2020 Homework #8
Due on Friday, May 29, 2020
Instructions for submitting:
• Problems 1-3 should be submitted to MATLAB Grader. You have 3 attempts for each problem.
• Problems 4-5 should be submitted to Gradescope. The solutions and the code used to get those solutions should be submitted as a single pdf. All code should be at the end of the file. You must select which page each problem is on when you submit to Gradescope.
Note: All three MATLAB Grader problems in the homework have corresponding data held in a .mat file. You should download them to the same directory as your script file and access them by using the load command. MATLAB Grader has its own copy of the data files.
(25 points) Problem 1: MATLAB Grader
In order to determine the half-life of Plutonium-239, scientists start with a sample of approximately 50 kg of Plutonium-239 and measure the remaining amount each year for 40 years. The data is contained in the file Plutonium.mat which is included with the homework. Plutonium.mat contains two vectors, t and P. The vector t contains the number of years since the beginning of the experiment, and the vector P contains the corresponding amounts of Plutonium-239 remaining measured in kg. Let the function P (t) denote the amount of remaining Plutonium as a function of time.
(a) Use a central difference to approximate the derivative dP at time t = 10. Store the dt
result in the variable ans1.
(b) Use the second order difference formula that was derived in Activity 8 to approximate
the derivative dP at time t = 0. Store the result in the variable ans2. dt
(c) Use a second order accurate difference scheme to approximate the derivative dP at dt
time t = 40. Store the result in the variable ans3.
(d) The decay rate of Plutonium-239 at a time t is given by − 1 dP . Use a second order P dt
differenceschemetoapproximatedP atall41timesint.Youshoulduseyouranswers dt
to (b) and (c) for times t = 0 and t = 40 and a central difference for every other time. Then use these approximations to estimate the decay rate at all 41 times in t. Create a 1 × 41 row vector named ans4 with the estimates at each time in chronological order.
(e) If λ is the average of the decay rates that you found in part (d), then the half-life of Plutonium-239, denoted by t1/2, is given by the formula
t1/2 = ln(2). λ
Calculate the half-life, and store it in the variable ans5.
(15 points) Problem 2: MATLAB Grader
Consider blood flow through an artery or vein. For laminar flow (i.e. flow in which the fluid moves in parallel layers), the velocity of the blood is given by the equation
v(r)= ∆p(R2−r2), 4μL
where r is the distance from the center of the blood vessel, R is the radius of the blood vessel (distance from the center to the wall), ∆p is the change in pressure from the beginning of the blood vessel to the end, μ is the viscosity of the blood, and L is the length of the blood vessel. By examining this equation, we can see that the blood moves fastest in the center of the blood vessel (r = 0) and slowest near the walls (r is close to R). The volumetric flow rate Q (the volume of fluid that passes through a cross section per unit time) is given by
R 0
If all of the parameters in the function v(r) are known, this is a very easy integral to evaluate by hand. However, the parameters are often not known. Instead, measurements of the velocity can be taken at different values of r.
(a) The file BloodFlow.mat contains two vectors, r and v. The vector r contains the values of r at which measurements were taken, and the vector v contains the corre- sponding velocities in m/s. Using this data, use the trapezoidal rule to approximate the volumetric flow rate Q and store the result in the variable ans6.
Q =
2πrv(r)dr.
(b) The cross-sectional area A of the blood vessel can be calculated with the integral R
2πrdr.
Approximate this integral using the vector r and the right-sided rectangle rule. Store
the result in the variable ans7.
(c) The mean velocity is the volumetric flow rate Q divided by the cross-sectional area A. Calculate the mean velocity using your answers to part (a) and (b) and store the result in the variable ans8.
(10 points) Problem 3: MATLAB Grader
The cardiac output C0 is the volume of blood pumped by the heart per unit time. One way to measure cardiac output is to inject dye into the right atrium and then measure the concentration of dye in the blood that is leaving the heart. If A is the amount of injected dye and c(t) is the concentration of dye as a function of time, then the cardiac output is
A
C0 = T c(t)dt
0
where T is the final measurement time. The file Dye.mat contains two vectors. The vector t contains the times (in seconds) at which the dye concentration was measured, and the vector c contains the concentration of dye in the blood at those times. Use Simpson’s rule to evaluate the integral in the denominator of the formula above. Store the answer in the variable ans9. Then use the value A = 3 ml to calculate the cardiac output in ml/s. Store the result in the variable ans10.
(45 points) Problem 4: Gradescope
Adult Alaskan Malamutes (i.e. Huskies) have weights that are normally distributed with a mean of 85 pounds and a standard deviation of 5 pounds. To compute the probability that a randomly selected Malamute has a weight between 76 and 86 pounds, you would compute the integral
86 1 −(x−85)2/50 P=√e dx
50π
This integral cannot be evaluated exactly by using any of the methods you learned in
Calculus class so we will evaluate it using numerical integration.
A =
0
76
(a) Use the integral function to calculate the “exact” value of P .
(b) Use the left-sided rectangle rule to approximate P with step sizes of
h = 1, 2−1, 2−2, . . . , 2−16. Calculate the error for each h by taking the absolute value of the exact solution minus the approximation. Store the errors in a vector.
(c) Do the same for the right-sided rectangle rule, the midpoint rule, the trapezoidal rule, and Simpson’s rule.
(d) Plot the errors versus h on a log-log plot. Use a different color and marker type for each method. Plot a trend line that represents O(h) by plotting c · h versus h on the log-log plot. Choose the constant c so that the trend line falls near your error points. Also include trend lines for O(h2) and any other orders that are represented by the numerical integration methods in your plot. Use different line styles for each trend line. Below is a sample of what the plot might look like for just two different integration schemes and two trend lines.
100
10-5
10-10
10-15
10-5 10-4
10-3 10-2 Grid Spacing, h
10-1 100
Convergence of Numerical Integration Schemes
Method Name 1 Method Name 2 O(h)
O(h2)
Also add a horizontal line at 10−16 which is (approximately) “machine precision”. This is the lowest you could reasonably expect the error of one of the methods to be because of rounding error. Add appropriate labels to the x and y axes, a legend, and a title. You will be graded on how easy it is to see and interpret your plot and how well it illustrates the orders of each method.
Error
(5 points) Problem 5: Gradescope
Which of the following finite difference approximations have error that is O(∆x2)? Choose all that apply.
(a) f′(x) = f(x+∆x)−f(x−∆x) 2∆x
(b) f′(x) = f(x)−f(x−∆x) ∆x
(c) f′(x) = f(x+∆x)−f(x) ∆x
(d) f′(x) = −3f(x)+4f(x+∆x)−f(x+2∆x) 2∆x
(e) There is not enough information because the problem doesn’t specify whether it is local error or global error.
Gradescope Deliverables Your Gradescope writeup should contain the following:
• Problem 4: The plot
• Problem 5: A letter or multiple letters corresponding to your answer choice • Code: Code for problem 4