Mathematical Modelling - Python knowledge required
Homework 2, Math 450, Spring 2018
Least Squares, Scaling laws, Differential equations
Instruction: Answer each of the following questions, showing and explaining your work as you go. Partial credit will be awarded based on how well I can follow your work and how far you get, so please use sentences, description, diagrams, and clear definitions to communicate your results as best you can. All diagrams and plots should be labelled, for example.
1. The calculations Chauvenet actually performs in his book are a little more complicated than described above. Based on other information, some measurements were more trustworthy than others. Specifically, the measurement of the angle between Burden to Jescelyne is believed to be less precise than the other three angles. It would be nice to be able to take this imprecision into account when calculating our least-squares solution. And there is indeed a way -- we can weight the errors produced by each equation such that equations measured to greater precision are given larger weights and equations measured to lower precision are given smaller weights.
(a) Find the weighted normal equations in matrix form for the vector xbest that minimizes the weighted square error
E(x) = ∑ i
Wii((Ax)i − bi)2,
where W is a non-negative diagonal matrix. Answer: Observing that
[ p1 p2
][m11 0 0 m22
][ p1 p2
] = m11p
2 1 + m22p
2 2,
it follows that in matrix-vector form, E(x) = (Ax− b)T W(Ax− b). So, using matrix algebra rules,
E(x) = (Ax− b)T W(Ax− b)
= (x T A
T − bT )W(Ax− b)
= x T A
T WAx−xT AT Wb− bT WAx + bT Wb
And since each term in this sum is a scalar (not a vector or matrix) xT AT Wb = (xT AT Wb)T = bT WAx, so
E(x) = x T A
T WAx−2xT AT Wb + bT Wb.
Now, if we calculate the gradient with respect to x, we find
dE(x)
dx = 2A
T WAx−2AT Wb.
Setting equal, and rearranging, we find the weighted normal equations
A T WAx = A
T Wb.
(b) Chauvenet weights the first 3 measurements as 3 times more precise than the 4th, and the observation that angles must sum to 360 degrees with infinitely more weight. Replicate the calculations in the published version of Chauvenet’s survey problem using the weighted normal equations obtained above. Explain any discrepancies with Chauvenet’s result. Answer: So, as described above, we want to resolve the linear least squares problem with weights 3, 3, 3, 1,∞. We cann’t actually use ∞, so let’s just pick a big number like 106 instead. Chauvenet makes a calculation mistake, and finds
b =
16.65916.509 16.457
1
instead of
b =
16.65916.509 16.543
2. In 1807, Nathaniel Bowditch observed a comet in the New England sky and made a series of astronomical observations. He then attempted to calculate the orbit of the comet using the methods described by Laplace. In the process, he developed a system of 56 equations for 5 unknowns. Use linear least squares to estimate the values of these 5 parameters and discuss the relationship between your estimates and Bowditch’s. Using your estimates, plot the residuals b − Axbest, plotted by equation number. The data has been extracted to here.
3. Table 2, Page 45 of J. Huxleys 1932 book Problems of Relative Growth provides a table of red deer antler statistics relative to body mass. Use python to recalculated a scaling relationship between the deer mass and the antler mass from this data. In your answer, present the data in a new table, your parameterization for the scaling law, your fitted parameter values, and a plot comparing the data to the fitted scaling law.
4. A string was wrapped around a metal bar. A 2 ounce weight was hung at one end of the string, and weights W at the other, so as to counterbalance the weight of 2 oz. The more revolutions of string used, the greater an imbalance the friction between string and bar can sustain. Find a law relating the number of revolutions to the maximum sustained weight. Possibilities include lines, power laws, and exponential laws - chose the one you think works the best. Include plots of the data and your fit using Cartesian, semilog and log-log plots. (See https://www.math.psu.edu/treluga/450/stringfriction.csv)
Answer: The data look most like a line when plotting the abscissa against the log of the ordinate. The best-fit exponential curve is given by
log y = 0.91237 + 1.14588x.
5. In a 2009 article, Head and other authors use a fossil vertebrae to deduce that about 60 million years ago, there existed a neotropical snake that grew to 13 meters in length and weighted more than 1,000 kilograms (42 feet and more than 2,000 pounds). They called this snake Titanoboa cerrejonensis. This discovery and others like it inspired the PBS Secrets of the Dead episode Graveyard of the giant beasts in 2016, also found on (youtube). Let’s see if we agree with the hype by looking through this scientific paper and check their numbers.
(a) What were the two regression lines the team obtained for the relationship between between vertebra width and total body length using the extreme values of 60% and 65% for vertebral position? Answer: The paper obtains 4 regression lines, but only two are for total body length: y = 100.7x + 436.2 and y = 106.0x + 390.
(b) Create a plot of width vs total body length using known 60% and 65% position data. Add to this plot the regression lines from part (a). Make sure your axes scales are large enough to show the predicted body length of Titanoboa, given that the discovered vertebra was 12 cm wide. Answer:
2
0
2000
4000
6000
8000
10000
12000
14000
0 20 40 60 80 100 120
T o
ta l b
o d
y l e
n g
th (
m m
)
Vertabral width (mm)
(c) Do you think Head et al’s conclusion about body length are supported by the evidence, are contradicted by the evidence, or that the evidence is inconclusive? Defend your opinion.
6. One day it started snowing at heavy and steady rate. A snowplow started out at noon, going 2 miles the first hour and 1 mile the second hour. What time did it start snowing? Answer: Let’s suppose it started snowing at time t0 at a rate a inches per hour. Then at time t, there were a(t− t0) inches on the ground. Now, from the description given, we now that the snowplow is slowed down the snow -- the more snow there is, the slower it will go. A natural conjecture, then, is that the speed of the snowplow (say, dx/dt) is inversely proportional to the depth of the snow:
dx
dt =
k
a(t− t0)
for proportionality constant k. Then, integrating and using x(12) = 0 since the snowplow starts out at noon,
x(t) = k
a log
( t− t0 12− t0
) Then, we can use our observations of x(13) and x(14) to determine t0.
7. (Letter by Thomas Godfrey to the Ben Franklin’s Pennsylvania Gazette, 1743) Mr. Franklin, As privateering is now so much in Fashion, the printing of the following question may be an amusement, if not to the privateers, yet to some of your correspondents or readers. Suppose a privateer, in the latitude of 10 degrees North, should at 6 in the morning spy a ship due south of her, distant 20 miles; upon which she steers directly for her, and runs at the rate of 8 miles an hour. The ship at the same time sees the privateer, but not being much afraid of her, keeps on her course due west, and sails at the rate of 6 miles an hour; how many hours will it be before the privateer overtakes the ship?
Answer: This problem is a ‘‘pursuit’’ problem, and the solution is called a ‘‘curve-of-pursuit’’. These problems were a bit of a rage back in 1740’s after Frenchman Pierre Bouguer had revealed their solution. It would have been extraordinary for Godfrey’s problem to have received any solutions from the newspaper’s readers. Suppose the merchant ship starts at position 0, and so, has travelled a distance 6t after t hours. Let (x(t), y(t)) be the position of the privateer at time t relative to the merchant ship’s initial position, so at t = 0, (x(0), y(0)) = (0, 20). Since the privateer always steers in the direction of the merchant ship, then by similar triangles,
−ẋ 6t + x
= −ẏ y
,
or, after cross-multiplying, −yẋ + ẏ(6t + x) = 0.
Now comes the fancy parts. Let’s reparameterize our problem in terms of y instead of t. We can divide through by ẏ to get
−y dx
dy + 6t(y) + x(y) = 0.
This looks okay, but we’ve got this weird t(y) term.
−y d2x
dy2 −
dx
dy + 6
dt
dy +
dx
dy = 0.
Since ẋ 2 + ẏ
2 = 8
2
3
1
8
√( dx
dy
)2 + 1 =
dt
dy
Substituting,
−y d2x
dy2 +
6
8
√( dx
dy
)2 + 1 = 0.
This is a second order equation, but can be transformed into a first order equation by taking f(y) = dx/dy. The resulting equation is separable and the solution can then be obtained by integration.
x(y) = y
2
[( y
A
)v − ( A
y
)v] + C
8. Plot the solution x(t) of our one-compartment acetaminophen model when x(0) = 0 and i(t) = 100(1 − cos(2πt/6)) over two days using the estimate that r = 0.4. (remember to always label your axes including units of measurement, include legends when needed, and make sure fonts are legible when printed)
9. Plot the solution x1(t) of our two-compartment acetaminophen model when x0(0) = x1(0) = i1(t) = 0 and i0(t) = 100(1 − cos(2πt/6)) over two days using the estimate that r0 = 0.15, r1 = 0.4, m01 = 2.8. What will be the estimated amount of acetaminophen in the gut and the body at t = 30 hours? (Give your answers to 5 decimal places.)
10. Draw a compartmental model for acetaminophen kinetics where a 3rd compartment (the liver) is included between the gut and blood. Be sure to label the compartments with their state variables and names, and label the arros with the parameters representing their rate constants. Then write out the corresponding system of ordinary differential equations.
11. Numerically solve our linearized HCV model
dI
dt = cδ
p V − δI
dV
dt = (1 − �)pI − cV
with parameter values c = 4.6, δ = 0.35, � = 0.95, and p = 50. Use initial conditions V0 = 10 6
and I0 = 9.2 × 104 over 14 hours. Plot I(t) and V (t) on separate subplots using ’semilogy()’. Also include a plot of phase-plane trajectory of I(t) vs V (t). Discuss how your numerical results compare with the patient-data presented in class.
4