Have you ever wondered how GPS works? Reality check 4 on page 238 of Numerical Analysis, Second Edition by Timothy Sauer challenges you to answer this question. GPS consists of 24 satellites carrying atomic clocks which orbit the earth at an altitude of 20,200 km.These satellites transmit carefully synchronized signals from space to GPS recievers on earth. These receivers take the signals and translate them via mathematics to determine the (x,y,z) coordinates of the receiver.
At a given instant, the receiver collects the signal from a satellite and determines the transmission time, the amount of time between when the signal was sent and when it was received. Theoretically, the speed of the signal is the speed of light so multiplying the transmission time by the speed of light should give the distance from the receiver to the satellite, putting the receiver on the surface of the sphere centered at the satallite position with a radius equal to the distance from the receiver to the satellite. If three satellite are used, then there are three spheres with the receiver at their intersection. However, the clocks in the receivers tend to be low accuracy which means that using only three satellites can result in position errors of several kilometers. Adding a fourth satellite to this problem fixes the error while complicating the mathematics.
Shown below is a simplified version of the GPS problem for four satellites. The position of each satellite in the sky is denoted by ( \(A_i,B_i,C_1 \)) while the transmission time is denoted \( t_i \). Define d to be difference between the synchronized time on the four satellites and the receiver clock. Then the intersection of the spheres, denoted (x,y,z), satisfies \[ \begin{align} r_1 (x,y,z,d) = \sqrt[2]{(x-A_1)^2 + (y- B_1)^2 + (z- C_1)^2}-c(t_1 -d) = 0 \\ r_2 (x,y,z,d) = \sqrt[2]{(x-A_2)^2 + (y- B_2)^2 + (z- C_2)^2}-c(t_2 -d) = 0 \\ r_3 (x,y,z,d) = \sqrt[2]{(x-A_3)^2 + (y- B_3)^2 + (z- C_3)^2}-c(t_3 -d) = 0 \\ r_4 (x,y,z,d) = \sqrt[2]{(x-A_4)^2 + (y- B_4)^2 + (z- C_4)^2}-c(t_4 -d) = 0 \end{align} \]
All work presented on this page was done jointly by Anna-Rose Wolff and Donnelly Phillips. Anna-Rose Wolff mainly encoded Newton's Method and related Jacobian code, debugged trials, developed website, and wrote the step descriptions. Donnelly Phillips coded step 2, produced trials for steps 4, 5, and 6, recorded data, created plots, and modified Newton's Method to create the Gauss-Newton Method.
First, the above system was solved using Multivariate Newtons Method.
The satellite positions in kilometers were set to be (15600, 7540, 20140), (18760, 2750, 18610), (17610, 14630, 13480), (19170, 610, 18390)
with transmission times in seconds of 0.07074, 0.07220, 0.07690, and 0.07242 respectively.
Lastly, the initial vector was set to \( (x_0,y_0,z_0,d_0) \) = (0, 0, 6370, 0).
The solution was found to be approximately (x,y,z) = (-41.77271, -16.78919, 6370.0596) kilometers with d = -0.003201566 seconds.
Our matlab code for Multivariate Newtons Method can be found here and
our corresponding code for the function we were solving and its jacobian can be found here
and here
It is actually possible to explicitly solve for a solution to this system of four equations, but quite a bit of algebra is involved. First, all of the quadratic terms can be eliminated by subtracting the first equation by the second, third, and fourth. This produces three equations with four unknowns, but without any quadratic terms. Next, determinant equations allow for x, y, and z to be solved for in terms of d. After this, by plugging in these equations for x, y, and z in the first equation, a quadratic equation is produced in terms of d. This can be solved by using the quadratic formula, bearing in mind that this will produce two solutions for d, which in turn will produce two coordinates for your position. However, only one of these coordinates will be on the surface of the Earth. when tested by using the previous settup, two coordinates were found, namely (-41.772709570837364, -6.789194106526850, 6370.059559223344) and (-39.747837348142937, -134.2741443606658, -9413.624553735754). The norm of the first was found to be 6370.218648080797, while the norm of the second is 9414.666041613707. The first one more closely resembles the radius of the Earth, so the first coordinate is the confirmed solution. The commented matlab code for this step with comments can be found here
Next we attempted to test the conditioning of the GPS problem. First, we defined the satellite positions \( (A_i,B_i,C_i) \) in terms of the spherical coordinates
\[
\begin{align}
A_i = \rho cos \phi_i cos \theta_i\\
B_i = \rho cos \phi_i sin \theta_i\\
C_i = \rho sin \phi_i
\end{align}
\]
where \( \rho = 26570\) kilometers was fixed and \( 0 \leq \phi_i \leq \pi/2 \) and \(0 \leq \theta_i \leq 2 \pi \) were chosen arbitrarily.
For the initial vector \( (x_0,y_0,z_0,d_0) \) = (0, 0, 6370, 0.0001) we defined the satellite ranges \( R_1 = \sqrt[2]{A_{i}^2 + B_{i}^2 + (C_1 - 6370)^2} \) and the
transmission times as \(t_i = d + R_i/c \).
Then the error magnification factor was defined to be the change in position divided by the speed of light times the change in transmission time.
First, a solution was found for four satellites with \( \phi = [3\pi/8,\pi/4,3\pi/8,\pi/4] \) and
\( \theta = [\pi/8,\pi/2,\pi,3\pi/2] \) defining the satellite's positions.
Then, the transmission times were changed by arbitrarily adding or subtracting \(10^{-8} \).
The table catalogs the changes and the results. The first column has plusses and minuses in the order that \(10^{-8} \) was either added or subtracted from the times. The second column has the changes of position, and the third column has the error magnification factor. All time is in seconds and all positions are in kilometers.
| Time Error | COP | EMF |
|---|---|---|
| + + + - | 0.008243337204476 | 2.749681313597343 |
| + + - - | 0.005669147821078 | 1.891024162998724 |
| + - - - | 0.009698382797069 | 3.235032281982867 |
| + - + + | 0.010415180993732 | 3.474130424499227 |
It is easily seen that the maximum positisition error found was 10.415180993732 meters and the condition number of the matrix seems to be around 3.
We repeated the procedure with a more tightly grouped set of satellites. In particular, we chose \(\phi \) and
\(\theta \) such that all the \( \phi_i \) were within five percent of each other and the \( \theta_i \) were within five percent of each other.
We solved with and without the input error as before.
The solution was calculated with \( \phi = [\pi/4,1.05\pi/4,1.02\pi/4,1.01\pi/4] \) and
\( \theta = [\pi/,1.05\pi,1.03\pi,1.01\pi] \) defining the satellite's positions.
| Time Error | COP | EMF |
|---|---|---|
| + + + - | 3.497561570176913 | 1.166660960555890e+03 |
| + + - - | 1.096393054388909 | 1.096393054388909e+03 |
| + - - - | 2.637937747787769 | 8.799213180298781e+02 |
| + - + + | 0.649185055055729 | 2.165448255058268e+02 |
The maximum position error in this case was found to be 3497.561570176913 meters which corresponded to an error magnification factor of about 1167. By bunching the satalites, the change in position is far more noticable. The observed error magnification factors and condition number are extremely large.
Next, three more satellites were added to the calculation in an attempt to reduce the position error and the condition number.
In order to be able to solve the resulting equation, the Gauss-Newton method was imployed and the associated
code can be found here.
We used \( \phi = [3 \pi/8,\pi/4,3 \pi/8,\pi/4,\pi/6,\pi/3,\pi/6]\) and
\( \theta = [\pi/8,\pi/2,\pi,3\pi/2,3\pi/4,5\pi/4,7\pi/4] \) with an initial vector of (x,y,z,d) = (0,0,6370,.0001).
| Time Error | COP | EMF |
|---|---|---|
| + + + - + + - | 0.005422503013506 | 1.808752311409307 |
| + + - - + + - | 0.003382873997695 | 1.128405304210651 |
| + - - - + + - | 0.001724615436017 | 0.575269787479879 |
| + - + + + + - | 0.004428214906511 | 1.477093498633355 |
The maximum position error was found to be 5.422503013506 meters with a error magnification value of 1.808752311409307.
By adding more satalites, the average change in position appears to be much smaller than
only using four satalites, and significantly smaller than using four bunched satalites.
The maximum error magnification factor found with seven satalites is close to the minimum found
with four. So the condition number of the seven satalites is much smaller than using four satalites.
Obviously, using seven unbunched satalites would produce the smallest error, and it can be speculated
that adding additional satalites would further reduce the error.
So, to summarize, if you were working off of four bunched satellites, you could experience a maximum GPS error of several
thousand meters but if you work off of seven satellites instead, your maximum GPS error drops to less than a
hundred meters.