Description
Phys 512 Problem Set 3
Due on github Friday September 30 at 11:59 PM. You may discuss problems, but everyone must write their own code.
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+x2
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 prob- lem. 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-squares- type 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 = x2/(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 symmet- ric. In real life, we’d obviously need to check that. The leading order correction would give us a dish that looked like z = ax′2 + by′2 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′,y′ by a rotation: x = cos(θ)x′ + sin(θ)y′ and y = − sin(θ)x′ + cos(θ)y′. 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?
2
