Phys 512 Problem Set 3

Phys 512 Problem Set 3
Problem 1: Write an RK4 integrator with prototype to take one step:
def rk4_step(fun,x,y,h):
Use this to integrate
dy/dx
=
y
1 + x
2
from x = −20 to x = 20 with y(−20) = 1 using 200 steps. Now write another
stepper
def rk4_stepd(fun,x,y,h):
that takes a step of length h, compares that to two steps of length h/2, and uses
them to cancel out the leading-order error term from RK4. How many function
evaluations per step does this one use? Use this modified stepper to carry out
the same ODE solution using the same number of function evaluations as the
original. Which is more accurate?
NB - the analytic solution to the equation can be found by separation of
variables, and is y = c0 exp(arctan(x)).
Problem 2: a) Write a program to solve for the decay products of U238
(refer to slides for the decay chain). You can use the ODE solver from scipy,
but you’ll need to set the problem up properly. Please make sure to include all
the decay prodcuts in the chain. Assume you start from a sample of pure U238
(in nature, this sort of separation happens chemically when rocks are formed).
Which solver would you use for this problem?
b) Plot the ratio of Pb206 to U238 as a function of time over a region where
it’s interesting. Does this make sense analytically? (If you look at the decay
chain, all the half-lives are short compared to U238, so you can approximate the
U238 decaying instantly to lead. Now plot the ratio of Thorium 230 to U234
over a region where that is interesting. Radioactive decay is frequently used
to date rocks, and these results point at how you can determine the age of a
uranium-bearing rock that is anywhere from thousands to billions of years old.
(Of course, in this case the starting ratio of U234 to U238 would probably have
already reached its long-term average when the rock was formed, but you could
still use the U234/Th230 ratio under that assumption.)
Problem 3: We’ll do a linear least-squares fit to some real data in this problem. Look at the file dish zenith.txt. This contains photogrammetry data for a
prototype telescope dish. Photogrammetry attempts to reconstruct surfaces by
working out the 3-dimensional positions of targets from many pictures (as an
1
aside, the algorithms behind photogrammetry are another fun least-squarestype problem, but beyond the scope of this class). The end result is that
dish zenith.txt contains the (x,y,z) positions in mm of a few hundred targets
placed on the dish. The ideal telescope dish should be a rotationally symmetric
paraboloid. We will try to measure the shape of that paraboloid, and see how
well we did.
a) Helpfully, I have oriented the points in the file so that the dish is pointing
in the +z direction (in the general problem, you would have to fit for direction
the dish is pointing in as well, but we will skip that here). For a rotationally
symmetric paraboloid, we know that
z − z0 = a

(x − x0)
2 + (y − y0)
2


and we need to solve for x0, y0, z0, and a. While at first glance this problem may
appear non-linear, show that we can pick a new set of parameters that make
the problem linear. What are these new parameters, and how do they relate to
the old ones?
b) Carry out the fit. What are your best-fit parameters?
c) Estimate the noise in the data, and from that, estimate the uncertainty in
a. Our target focal length was 1.5 metres. What did we actually get, and what
is the error bar? In case all facets of conic sections are not at your immediate
recall, a parabola that goes through (0, 0) can be written as y = x
2/(4f) where
f is the focal length. When calculating the error bar for the focal length, feel
free to approximate using a first-order Taylor expansion.
BONUS: Of course, we have just assumed that the dish is circularly symmetric. In real life, we’d obviously need to check that. The leading order correction
would give us a dish that looked like z = ax02+by02
if the vertex (bottom) of the
dish was at (0, 0, 0) and our coordinate system was aligned with the principal
axes of the dish. We won’t usually have the benefit of being aligned like that;
instead we’ll usually be rotated by some (unknown) angle θ, so our observed
coordinates x, y will be related to the original coordinates x
0
, y0 by a rotation:
x = cos(θ)x
0 + sin(θ)y
0 and y = − sin(θ)x
0 + cos(θ)y
0
. Find the focal lengths
of the two principal axes (and don’t forget we can still have arbitrary offsets
x0, y0, z0). Is the dish round?