Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

My submission for PHY450H1S, Relativistic Electrodynamics, Problem Set 1 (partial)

Posted by peeterjoot on February 1, 2011

[Click here for a PDF of this post with nicer formatting]


This problem set is as yet ungraded. It is also incomplete since I ended up hand-writing some parts.

Problem 2.


From the Lorentz transformations of space and time coordinates.

\item Derive the transformation of velocities. With a particle moving with \mathbf{v} in the unprimed (stationary) frame, find its velocity \mathbf{v}' in the primed frame. The primed frame is moving with some \mathbf{V} with respect to the unprimed one. Make sure to finally derive the general “addition of velocities” equation in terms of vectors and dot products, as given in [1].
\item Then, use the addition of velocities rule to show that: a) if v < c in one frame, then v'  c in one frame, than v' > c in any other frame.


Part 1.

We need a vector form of the Lorentz transform to start with. Let’s write \boldsymbol{\sigma} for a unit vector colinear with the primed frame velocity \mathbf{V}, so that \mathbf{V} = (\mathbf{V} \cdot \boldsymbol{\sigma}) \boldsymbol{\sigma}. When our boost was in the x direction, our Lorentz transformation was in terms of x = \mathbf{x} \cdot \hat{\mathbf{x}}. The component in the direction of the boost is now \mathbf{x} \cdot \boldsymbol{\sigma}, and we have


\begin{aligned}c t' &= \gamma \left( ct - (\mathbf{x} \cdot \boldsymbol{\sigma}) \frac{\mathbf{V} \cdot \boldsymbol{\sigma}}{c} \right) \\ \mathbf{x}' \cdot \boldsymbol{\sigma} &= \gamma \left( \mathbf{x} \cdot \boldsymbol{\sigma} - \frac{\mathbf{V} \cdot \boldsymbol{\sigma}}{c} c t \right) \\ \mathbf{x}' \wedge \boldsymbol{\sigma} &= \mathbf{x} \wedge \boldsymbol{\sigma} .\end{aligned} \hspace{\stretch{1}}(2.1a)


We can add the vector components using \mathbf{x} = (\mathbf{x} \cdot \boldsymbol{\sigma}) \boldsymbol{\sigma} + (\mathbf{x} \wedge \boldsymbol{\sigma}) \boldsymbol{\sigma}, leaving


\begin{aligned}c t' &= \gamma \left( ct - (\mathbf{x} \cdot \boldsymbol{\sigma}) \frac{\mathbf{V} \cdot \boldsymbol{\sigma}}{c} \right) \\ \mathbf{x}' &= (\mathbf{x} \wedge \boldsymbol{\sigma}) \boldsymbol{\sigma} + \gamma \left( (\mathbf{x} \cdot \boldsymbol{\sigma}) \boldsymbol{\sigma} - \frac{\mathbf{V}}{c} c t \right) .\end{aligned} \hspace{\stretch{1}}(2.2a)


Writing (\mathbf{x} \wedge \boldsymbol{\sigma}) \boldsymbol{\sigma} = \mathbf{x} - (\mathbf{x} \cdot \boldsymbol{\sigma})\boldsymbol{\sigma} we have for the spatial component transformation

\begin{aligned}\mathbf{x}' = \mathbf{x} + (\mathbf{x} \cdot \boldsymbol{\sigma}) \boldsymbol{\sigma} (\gamma - 1) - \gamma \frac{\mathbf{V}}{c} c t.\end{aligned} \hspace{\stretch{1}}(2.3)

Now we are set to take derivatives to calculate the velocities. This gives us

\begin{aligned}\frac{dt'}{dt} &= \gamma \left( 1 - \left( \frac{d\mathbf{x}}{dt} \cdot \boldsymbol{\sigma} \right) \frac{\mathbf{V} \cdot \boldsymbol{\sigma}}{c^2} \right) \\ \frac{d\mathbf{x}'}{dt'} \frac{d t'}{dt} &= \frac{d\mathbf{x}}{dt} + \left(\frac{d\mathbf{x}}{dt} \cdot \boldsymbol{\sigma}\right) \boldsymbol{\sigma} (\gamma - 1) - \gamma \frac{\mathbf{V}}{c} c .\end{aligned} \hspace{\stretch{1}}(2.4a)


Dividing this pair of equations, and using \mathbf{v} = d\mathbf{x}/dt, and \mathbf{v}' = d\mathbf{x}'/dt', this is

\begin{aligned}\mathbf{v}' = \frac{\gamma^{-1} \mathbf{v} + (\mathbf{v} \cdot \boldsymbol{\sigma}) \boldsymbol{\sigma} (1 - \gamma^{-1}) - \mathbf{V}}{ 1 - \left( \mathbf{v} \cdot \boldsymbol{\sigma} \right) (\mathbf{V} \cdot \boldsymbol{\sigma})/c^2 }.\end{aligned} \hspace{\stretch{1}}(2.5)

Since \mathbf{V} and our direction vector \boldsymbol{\sigma} are colinear, we have (\mathbf{v} \cdot \boldsymbol{\sigma}) (\mathbf{V} \cdot \boldsymbol{\sigma}) = \mathbf{v} \cdot \boldsymbol{\sigma}, and can simplify this last expression slightly

\begin{aligned}\mathbf{v}' = \frac{\gamma^{-1} \mathbf{v} + (\mathbf{v} \cdot \boldsymbol{\sigma}) \boldsymbol{\sigma} (1 - \gamma^{-1}) - \mathbf{V}}{ 1 - \mathbf{v} \cdot \mathbf{V}/c^2 }.\end{aligned} \hspace{\stretch{1}}(2.6)

Finally, if we are to compare to the text, we note that the inverse expression requires replacement of \mathbf{V} with -\mathbf{V} and switching \mathbf{v} with \mathbf{v}'. That gives us

\begin{aligned}\mathbf{v} = \frac{\gamma^{-1} \mathbf{v}' + (\mathbf{v}' \cdot \boldsymbol{\sigma}) \boldsymbol{\sigma} (1 - \gamma^{-1}) + \mathbf{V}}{ 1 + \mathbf{v}' \cdot \mathbf{V}/c^2 }.\end{aligned} \hspace{\stretch{1}}(2.7)

The expression in the text is also a small velocity approximation. For {\left\lvert{\mathbf{V}}\right\rvert} \ll c, we have \gamma^{-1} \approx 1, and (1 + \mathbf{v}' \cdot \mathbf{V}/c^2)^{-1} \approx 1 - \mathbf{v}' \cdot \mathbf{V}/c^2. This gives us

\begin{aligned}\mathbf{v} \approx (\mathbf{v}' + \mathbf{V})( 1 - \mathbf{v}' \cdot \mathbf{V}/c^2 ) \approx \mathbf{V} + \mathbf{v}' - \mathbf{v}' (\mathbf{v}' \cdot \mathbf{V})/c^2\end{aligned} \hspace{\stretch{1}}(2.8)

One additional approximation was made dropping the \mathbf{V} (\mathbf{v}' \cdot \mathbf{V})/c^2 term which is quadratic in \mathbf{V}/c, which leave us with equation 5.3 in the text as desired.

Part 2.

In 2.7, let’s write \mathbf{v}' = u \mathbf{u}, where \mathbf{u} is a unit vector, V = \mathbf{V} \cdot \boldsymbol{\sigma}, and \alpha = \mathbf{u} \cdot \boldsymbol{\sigma} for the direction cosine between the primed frame’s direction of motion and the particle’s velocity direction (also in the unprimed frame). The stationary frame’s particle velocity is then

\begin{aligned}\mathbf{v} = \frac{\gamma^{-1} u \mathbf{u} + u \alpha \boldsymbol{\sigma} (1 - \gamma^{-1}) + V \boldsymbol{\sigma}}{ 1 + \alpha u V/c^2 }.\end{aligned} \hspace{\stretch{1}}(2.9)

As a check, note that for 1 = \alpha = \mathbf{u} \cdot \boldsymbol{\sigma} = \cos(0), we recover the familiar addition of velocities formula

\begin{aligned}\mathbf{v} = \mathbf{u} \frac{u + V}{ 1 + u V/c^2 }.\end{aligned} \hspace{\stretch{1}}(2.10)

We want to put 2.9 into a form that renders it more tractable for general angles too. Factoring out the \gamma^{-1} term appears to do the job, yielding

\begin{aligned}\mathbf{v} = \frac{u \gamma^{-1} (\mathbf{u} -\alpha \boldsymbol{\sigma}) + (u \alpha + V) \boldsymbol{\sigma}}{ 1 + \alpha u V/c^2 }.\end{aligned} \hspace{\stretch{1}}(2.11)

After a bit of reduction and rearranging we can dot this with itself to calculate

\begin{aligned}\mathbf{v}^2 = \frac{V^2(1 - \alpha^2)(1 - u^2/c^2) + (u + \alpha V)^2}{ (1 + \alpha u V/c^2)^2 }\end{aligned} \hspace{\stretch{1}}(2.12)

Note that for u = c, we have \mathbf{v}^2 = c^2, regardless of the direction of \mathbf{V} with respect to the motion of the particle in the unprimed frame. This shouldn’t be surprising since this light like invariance is exactly what the Lorentz transformation is designed to maintain. It is however slightly comforting to know that the algebra appears to be still be kosher after all this. This also answers part (b) of this question, since we have tackled the v = c case in the primed frame, and seen that the speed remains v = c in the unprimed frame (and thus any frame moving at constant speed relative to another).

Observe that since 1 - \alpha^2 = \sin^2\theta, and u \le c, this is positive definite as expected. If one allowed u > c in some frame, our speed could go imaginary!

For the u  c cases, let x = u/c and y = V/c. This allows 2.12 to be casted in a simpler form

\begin{aligned}\mathbf{v}^2 = c^2 \frac{y^2 (1 - \alpha^2)(1 - x^2) + (x + \alpha y)^2}{ (1 + \alpha x y)^2 }\end{aligned} \hspace{\stretch{1}}(2.13)

We wish to verify that (a) given any x \in (-1,1), we have \mathbf{v}^2  1, we have \mathbf{v}^2 > c^2 for all y \in (-1,1), \alpha \in (-1,1).

Considering (a) first, this requires a demonstration that

\begin{aligned}y^2 (1 - \alpha^2)(1 - x^2) + (x + \alpha y)^2 < (1 + \alpha x y)^2 .\end{aligned} \hspace{\stretch{1}}(2.14)

Expanding out the products and canceling terms, we want to show that for (a) that if {\left\lvert{x}\right\rvert},{\left\lvert{y}\right\rvert} < 1 we have

\begin{aligned}x^2 (1 - y^2) + y^2 < 1,\end{aligned} \hspace{\stretch{1}}(2.15)

and for (c) that if {\left\lvert{x}\right\rvert} > 1, we have for any {\left\lvert{y}\right\rvert} < 1

\begin{aligned}x^2 (1 - y^2) + y^2 > 1.\end{aligned} \hspace{\stretch{1}}(2.16)

Observe that the frame velocity orientation direction cosines have completely dropped out, leaving just the (relative to c) velocity terms.

To get an initial feel for this function f(x,y) = x^2 (1 - y^2) + y^2, notice that when graphed we have a bowl with a minimum (zero) at the origin, and what appears to be a uniform value of one on the boundary (case (b)). Then provided {\left\lvert{y}\right\rvert} < 1 it appears that the function f increases monotonically to a value greater than one (case (c)). While looking at a plot is not any sort of rigorous proof, let's move on to some of the other problems for now, and return to this last loose thread later if time permits.

Problem 3.


A toy model of a GPS system has satellites moving in a straight line with constant velocity V_x and at a constant height h (measured, e.g., along the y-axis) above “ground” (the x-axis). The satellites broadcast the time in their rest frame as well as their location at a time of broadcast. Imagine a person on the ground receives simultaneously broadcasts from two satellites, A and B, reporting their locations x_A' and x_B' as well as times of broadcast (which happen to be equal), t_A' = t_B'.

\item Find a condition determining your position in x. Evaluate it to find your deviation from the midpoint between the satellites to first order in V_x/c.
\item For some real numbers, note that in reality there are 24 satellites, moving with V ~4 \text{km}/s, a distance R \approx 2.7 \times 10^4 \text{km}. Use these numbers and the result from the previous problem (assuming a flat Earth, to be sure…) to get an idea whether (special) relativistic effects are important for the typical modern GPS accuracy of order 10 m (or less)?

Discussion of the non-toy model.

It is fairly easy to find interesting info about the mechanisms that real GPS works using. NASA has a nice How does GPS work page [2], and How stuff works has a nice How GPS Receivers Work article [3]. Reading these one finds that the GPS clocks are actually kept synchronized. The typical GPS receiver obviously has a clock, since we have countdown timers for time until arrival, is that clock accurate enough compared to the satellite atomic clocks to be used for the GPS location algorithm? What is done in fact is to use the local receiver time to seed the iterative algorithms, allowing the local time to be calculated eventually with an accuracy that actually approaches that of the satellite’s atomic clocks. Some of the sources of error, like reflection of the signals, delaying them, and interference by atmosphere are also discussed in these articles. Also interesting is that there is a table lookup of the satellites position implemented in the GPS receivers. This table lookup is used to seed the iterative algorithms, and can be used to reduce calculation error.

Our basic GPS problem is to calculate the intersection of a number of “spherical” hypersurfaces. This is made more interesting by the fact that this is both a non-linear and an over-specified problem. Let’s consider the geometric problem to get an idea of how to set up this problem. Suppose that we have a set of k satellites, located at position \mathbf{p}_i, i \in [1,k], and we know that these are located with distance d_i from our position \mathbf{x}.

Our problem is then to find the simultaneous solution to the following set of equations

\begin{aligned}\begin{aligned}(\mathbf{x} - \mathbf{p}_1)^2 &= d_1^2 \\ (\mathbf{x} - \mathbf{p}_2)^2 &= d_2^2 \\ \vdots \\ (\mathbf{x} - \mathbf{k}_2)^2 &= d_k^2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.17)

Observe that even if we reduce this to a one dimensional problem in a single variable x, we still have a non-linear system

\begin{aligned}\begin{aligned}0 &= (x - p_1)^2 - d_1^2 \\ 0 &= (x - p_2)^2 - d_2^2 \\ \vdots \\ 0 &= (x - k_2)^2 - d_k^2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.18)

We also need to be aware of the fact that each of the positions p_i, and the respective distances d_i will in reality both have associated errors, so there is not likely any specific single value of x that “solves” this problem, unless it is setup in a contrived and perfect fashion. This intrinsic error, and the k equations, one unknown nature of the problem (or three unknowns for spatial, or four for spatial and time position) suggests a least squares approach, but it will have to be one that also incorporates iteration.

We can setup our problem in matrix form, where we are looking for a solution to

\begin{aligned}F(x) = \begin{bmatrix}F_i(x)\end{bmatrix}=\begin{bmatrix}{\left\lvert{\mathbf{x} - \mathbf{p}_1}\right\rvert} - c {\left\lvert{ t - t_1 }\right\rvert} \\ {\left\lvert{\mathbf{x} - \mathbf{p}_2}\right\rvert} - c {\left\lvert{ t - t_2 }\right\rvert} \\ \vdots \\ {\left\lvert{\mathbf{x} - \mathbf{p}_k}\right\rvert} - c {\left\lvert{ t - t_k }\right\rvert} \\ \end{bmatrix}= 0.\end{aligned} \hspace{\stretch{1}}(3.19)

We seek the spacetime event vector x = (c t, \mathbf{x}) for the spatial location and the exact local time at the location of the GPS receiver. Given any approximation of the solution, we can refine the solution using Newton’s root finding method by taking partials, forming the Jacobian matrix for our function F(x). That is

\begin{aligned}F(x_0 + \Delta x) \approx F(x_0) + {\left.{{\frac{\partial {F_i(x)}}{\partial {x^j}}}}\right\vert}_{{x_0}} \Delta x^j = F(x_0) + J(x_0) \Delta x = 0\end{aligned} \hspace{\stretch{1}}(3.20)

This leaves us with the our least squares problem, requiring the generalized inverse to the matrix equation

\begin{aligned}x_1 = x_0 - J^\dagger (x_0) F(x_0),\end{aligned} \hspace{\stretch{1}}(3.21)


\begin{aligned}J^\dagger = (J^\text{T} J)^{-1} J^\text{T}.\end{aligned} \hspace{\stretch{1}}(3.22)

This is a solution in the least squares sense that given b = J x, the norm {\left\lvert{ J \bar{x} - b}\right\rvert} is minimized by \bar{x} = J^\dagger b.

This iterative method of solution, in the context of finding fitting circles and ellipses can be found discussed in detail in [4].

Solution. Part 1.

For our toy model we have two satellites A and B both moving in the positive x-axis direction at velocity V_x at height h. As seen above, we do not require the velocity of the satellites to setup the problem, and could express the problem to solve as the numerical solution of the set of equations

\begin{aligned}F(x) = \begin{bmatrix}\sqrt{ h^2 + (x - x_A')^2 } - c {\left\lvert{t - t_A'}\right\rvert} \\ \sqrt{ h^2 + (x - x_B')^2 } - c {\left\lvert{t - t_B'}\right\rvert} \\ \end{bmatrix} = 0.\end{aligned} \hspace{\stretch{1}}(3.23)

If we assume that our GPS receiver’s clock is synchronized sufficiently with satellites A and B, this single variable problem admits a closed form for one iteration of the least squares process. However, since we are asked for a result that includes a V_x/c term, we can augment our matrix equation by two additional rows, with a secondary set of data points introduced at an offset time interval. That is

\begin{aligned}F(x) = \begin{bmatrix}\sqrt{ h^2 + (x - x_A')^2 } - {\left\lvert{c t - c t_A'}\right\rvert} \\ \sqrt{ h^2 + (x - x_B')^2 } - {\left\lvert{c t - c t_B'}\right\rvert} \\ \sqrt{ h^2 + (x - x_A' - (V_x/c) c \delta t)^2 } - {\left\lvert{c t - c t_A' - c \delta t}\right\rvert} \\ \sqrt{ h^2 + (x - x_B' - (V_x/c) c \delta t)^2 } - {\left\lvert{c t - c t_B' - c \delta t}\right\rvert} \\ \end{bmatrix} = 0.\end{aligned} \hspace{\stretch{1}}(3.24)

With t, \delta t, x_A', and x_B' given, and an initial seed value for the iterative procedure assumed to be the midpoint x_0 = (x_A' + x_B')/2, we can calculate a first approximation to the receiver position x_1 = x_0 + \Delta x using the Newton’s procedure outlined above.

For this system our Jacobian elements are all differentials of the following form

\begin{aligned}\frac{\partial {}}{\partial {x}} \sqrt{ h^2 + (x - p)^2 } - \Delta = \frac{x-p}{\sqrt{h^2 + (x-p)^2}},\end{aligned} \hspace{\stretch{1}}(3.25)

so, our Jacobian is

\begin{aligned}J = \begin{bmatrix}\frac{x - x_A'}{\sqrt{ h^2 + (x - x_A')^2 } } \\ \frac{x - x_B'}{\sqrt{ h^2 + (x - x_B')^2 } } \\ \frac{x - x_A' - (V_x/c) c \delta t}{\sqrt{ h^2 + (x - x_A' - (V_x/c) c \delta t)^2 } } \\ \frac{x - x_B' - (V_x/c) c \delta t}{\sqrt{ h^2 + (x - x_B' - (V_x/c) c \delta t)^2 } } \\ \end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.26)

Our deviation from the midpoint to first order in V_x/c is

\begin{aligned}\Delta x = - {\left.{{ \frac{J^\text{T} F}{J^T J} }}\right\vert}_{{x = (x_A' + x_B')/2}}\end{aligned} \hspace{\stretch{1}}(3.27)

To tidy this up, let

\begin{aligned}s = \frac{1}{{2}} (x_B' - x_A')\end{aligned} \hspace{\stretch{1}}(3.28)

\begin{aligned}D = (V_x/c) c \delta t\end{aligned} \hspace{\stretch{1}}(3.29)

\begin{aligned}{\left.{{J^\text{T} J}}\right\vert}_{{(x_A' + x_B')/2}} &= \frac{s^2}{h^2 + s^2}+\frac{(-s)^2}{h^2 + s^2}+\frac{(s - D)^2}{h^2 + (s -D)^2}+\frac{(-s - D)^2}{h^2 + (-s -D)^2} \\ &= \frac{2 s^2}{h^2 + s^2}+\frac{(s - D)^2}{h^2 + (s -D)^2}+\frac{(s + D)^2}{h^2 + (s +D)^2}\end{aligned}


\begin{aligned}{\left.{{J^\text{T} F}}\right\vert}_{{(x_A' + x_B')/2}} &=\begin{bmatrix}\frac{s}{\sqrt{ h^2 + (s)^2 } } &\frac{-s}{\sqrt{ h^2 + (-s)^2 } } &\frac{s - D}{\sqrt{ h^2 + (s - D)^2 } } &\frac{-s - D}{\sqrt{ h^2 + (-s - D)^2 } } \end{bmatrix}\begin{bmatrix}\sqrt{ h^2 + (s)^2 } - {\left\lvert{c t - c t_A'}\right\rvert} \\ \sqrt{ h^2 + (-s)^2 } - {\left\lvert{c t - c t_B'}\right\rvert} \\ \sqrt{ h^2 + (s - D)^2 } - {\left\lvert{c t - c t_A' - c \delta t}\right\rvert} \\ \sqrt{ h^2 + (-s - D)^2 } - {\left\lvert{c t - c t_B' - c \delta t}\right\rvert} \\ \end{bmatrix} \\ &=\frac{s}{\sqrt{s^2 + D^2}} \left( {\left\lvert{ ct - c t_B' }\right\rvert} - {\left\lvert{ c t - c t_A'}\right\rvert} \right) -2 D- \frac{(s - D){\left\lvert{c t - c t_A' - c \delta t}\right\rvert} }{ \sqrt{ h^2 + (s - D)^2 } } + \frac{(s + D){\left\lvert{c t - c t_B' - c \delta t}\right\rvert} }{ \sqrt{ h^2 + (s + D)^2 } } \end{aligned}

The final beastly ugly result, utilizing the helper variables of 3.29, and 3.28, we have for the deviation from the midpoint (after one iteration of this least squares Newton’s method):

\begin{aligned}\Delta x =\frac{   -\frac{s}{\sqrt{s^2 + D^2}} \left( {\left\lvert{ ct - c t_B' }\right\rvert} - {\left\lvert{ c t - c t_A'}\right\rvert} \right) +2 D   + \frac{(s - D){\left\lvert{c t - c t_A' - c \delta t}\right\rvert} }{ \sqrt{ h^2 + (s - D)^2 } }    - \frac{(s + D){\left\lvert{c t - c t_B' - c \delta t}\right\rvert} }{ \sqrt{ h^2 + (s + D)^2 } } }{   \frac{2 s^2}{h^2 + s^2}   +\frac{(s - D)^2}{h^2 + (s -D)^2}   +\frac{(s + D)^2}{h^2 + (s +D)^2}}\end{aligned} \hspace{\stretch{1}}(3.30)

In practice it doesn’t make much sense to compute this. You’ll want to use a computer, and assuming the availability of a pre-canned SVD routine to compute the generalized inverse, the toy model wouldn’t be any easier to solve than the real thing.


[1] L.D. Landau and E.M. Lifshits. The classical theory of fields. 1971.

[2] unknown. How does gps work [online]. 2011. [Online; accessed 24-Jan-2011].

[3] unknown. How gps receivers work [online]. 2011. [Online; accessed 24-Jan-2011].

[4] W. Gander, G.H. Golub, and R. Strebel. Least-squares fitting of circles and ellipses. BIT Numerical Mathematics, 34(4):558–578, 1994.


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: