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?

## Reviews

There are no reviews yet.