PHYS512 Assignment 3 -Solved

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 modifified 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 fifit to some real data in this prob
lem. Look at the fifile 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 
1aside, 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 fifile so that the dish is pointing 
in the +z direction (in the general problem, you would have to fifit 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