Mathematics of Global Positioning System
The basic mathematical idea described in terms of simple English goes something like this. First of think in two dimensions. We have three satellites represented as circles with known center coordinates and unknown radius. The point where the circles intercept will be our position. The radius is unknown since GPS receiver clocks are inaccurate. We thus want to solve the coordinate where the circles intercept as well as a time correction factor. To appreciate the importance of correcting for time one only need to remember that light travels around Earth several times in one second. Newtons method for solving systems of variables is one way of finding a solution to our position. This method has the advantage of quadratic convergence. For the real thing we need to think in three dimensions. Each satellite is represented by a sphere, we need four spheres to reduce the intercept to a point.
The actual mathematics go as follows. For each satellite we have an equation of the form:
(x_{1}  sx_{n})^{2} + (x_{2}  sy_{n})^{2} + (x_{3}  sz_{n})^{2} + c^{2}(sd_{n}  x_{4})^{2} = 0
sx, sy, sz and sd differs for each satellite, they represent the thee coordinates and measured time delay respectively. c is the speed of light. We need to solve our position; (x_{1}, x_{2}, x_{3}) as well as the error in the receiver clock x_{4}. That can be done with the following algorithm:

Arrange the satellite equations into a system S_{n,j}(x) with (n,j) = 1 ... 4. If we have five equations we can contruct a Coefficient Matrix by subtracting one of the equations from the rest.

Make an initial guess for position (x_{1}, x_{2}, x_{3}) and time correction factor x_{4}. Choosing position at the origin and correction so as not to cause any negative delays is appropriate.

Calculate the Jacobian Matrix J_{n,j} for S_{n} with (n,j) = 1 ... 4, more equations result in an overdetermined system. There are 24 satellites in both the GPS and GLONASS systems.

Solve the n by n Linear System J y = S for vector y, where we have substituted into S(x) our initial guess for x.

Update the variables: x = x + y.

Check for convergence, x_{new}^{2}  x_{old}^{2} < tolerance, otherwise return to step two.
On completion you know not only your position, you may update your clock as well.