Numerical Methods

Numerical Methods

Numerical methods provide powerful tools through computer algorithms to solve problems which cannot be solved by hand. The least squares method is a technique in scientific computing that solves the system of equations, Ax = b, when it is inconsistent. In this project, we use the Newton method to solve the inconsistent system of equations algorithmically. 

We are tasked to find the coordinates of a GPS receiver on earth given the locations of 3 satellites. At any given time, each satellite transmits a signal to the receiver. The receiver calculates the time difference between when the transmission is sent and when it is received through the use of clocks on the satellites and receiver. By calculating the time between when the transmission is sent and when it is received, the receiver can calculate the distance between itself and the satellites by multiplying the times by the speed of light, c. The 3 satellites can be viewed as spheres in the 3D plane. There are 2 intersections between these spheres, one being the location of the receiver, while the other being far from the Earth’s surface and is disregarded. 

However, one problem that arises is that the clock on the receiver is often cheaply made and not accurate to the degree of the satellite (for example, the clock on a cellphone). We solve this problem by introducing a 4th satellite, and thus a 4th equation to our system of equations. We explore the effects of which below: 

We are asked to solve the following system of equations, using the Multivariate Newton’s Method, with respect to the variables x, y, z, and d: 

^p(x − A₁)² + (y − B₁)² + (z − C₁)² = c(t₁ − d) 

^p(x − A₂)² + (y − B₁)² + (z − C₁)² = c(t₁ − d) 

^p(x − A₂)² + (y − B₁)² + (z − C₁)² = c(t₁ − d) 

^p(x − A₂)² + (y − B₁)² + (z − C₁)² = c(t₁ − d) 

In this system, c = 299,792.458 km/sec is the speed of light, and (A_i, B_i, C_i) are the coordinates of the four satellites, and t_i represents the transmission times of the signals from the satellites to the receiver. The variable d is the time drift of the receiver, and the variables (x, y, z) represent the coordinates of the position of the receiver. 

This system can also be written more conveniently as: 

(x − A₁)² + (y − B₁)² + (z − C₁)² = [c(t₁ − d)]² 

(x − A₂)² + (y − B₁)² + (z − C₁)² = [c(t₁ − d)]² 

(x − A₂)² + (y − B₁)² + (z − C₁)² = [c(t₁ − d)]² 

(x − A₂)² + (y − B₁)² + (z − C₁)² = [c(t₁ − d)]² 

We first solve this system for given inputs for the variables (A_i, B_i, C_i) and t_i, as described in the Methodology section (Section 2, below). 

Then we find arbitrary φ_is, for 0 < φ < π/2, and arbitrary θ_is, for 0 < θ < 2π. We set i = 1, …, 4. Now we use these φ_is θ_is to define (A_i, B_i, C_i) and t_i, using the following definitions: 

A_i = ρcos(φ_i)cos(θ_i) 

B_i = ρcos(φ_i)sin(θ_i) 

q 

R_i = 

C_i = ρsin(φ_i) 

A²_i + B²_i + (C_i − 6370)² 2

t_i = d + R_i/c 

We then solve the system of equations for this input. We then add an input error, randomized between ∆t_i = 10⁻⁸ and ∆t_i = −10⁻⁸, to the time intervals t_i. Then we run the Multivariate Newton’s Method to solve the system. When the Multivariate Newton’s Method returns, we are asked to find the forward error k∆x, ∆y, ∆zk_∞ and the error magnification factor k∆x,∆y,∆zk_∞ 

ck∆t₁,∆t₂,∆t₃,∆t₄k_∞. 

We are also asked to find the maximum forward error and the condition number. 

Methodology 

Multivariate Newton’s Method [2] 

We created a function called SatFunct – with inputs A_i, B_i, C_i, t_i, the speed of light, and our initial guess vector (x₀, y₀, z₀, d₀) – that implements the Multivariate Newton’s Method, and returns our solution. 

To implement the Multivariate Newton’s Method, we must prepare for its corrective step: x_i+1 = x_i − DF(X_i)⁻¹F(X_i) 

where F is our function vector, and DF is the Jacobian matrix of partial derivatives of our function, and X_iis the most up-to-date approximation of the solution-vector (from the previous iteration). 

Instead of calculating DF ⁻¹ we may define DF(X_i)⁻¹F(x_i) = s. When we multiply both sides by DF(X_i), we get F(X_i) = DF(X_i)s. Solving for s, using Matlab’s backslash operator, we can write our corrective step as: 

X_i+1 = X_i − s 

To solve for s, we define 

F(X_i) = 

and 

(x − A₁)² + (y − B₁)² + (z − C₁)² − [c(t₁ − d)]^{2 } 

(x − A₂)² + (y − B₂)² + (z − C₂)² − [c(t₂ − d)]² (x − A₃)² + (y − B₃)² + (z − C₃)² − [c(t₃ − d)]² (x − A₄)² + (y − B₄)² + (z − C₄)² − [c(t₄ − d)]² 

DF(X_i) = 

2(x − A₁) 2(y − B₁) 2(z − C₁) 2(c²(t₁ − d))

2(x − A₂) 2(y − B₂) 2(z − C₂) 2(c²(t₂ − d)) 2(x − A₃) 2(y − B₃) 2(z − C₃) 2(c²(t₃ − d)) 2(x − A₄) 2(y − B₄) 2(z − C₄) 2(c²(t₄ − d)) 

We set up a for-loop to run for many iterations, to refine the solution. However, we terminate the for-loop if the tolerance is small enough. In fact, I decided to code two tolerance-number criteria. One criterion tests how close to accurate our solution is. If kf(x_i)k₂ < tolerance then we say our solution is accurate enough, and we terminate the for-loop. I set this tolerance-number to 10⁻⁸. The other tolerance-number criterion tests how much our latest step refines the solution. If kx_i+1 − x_ik₂ < tolerance then we say that we have gotten to the point that refining our solution has too small of an impact for us, and we terminate the for-loop. I set this tolerance-number to 10⁻¹⁰. 

Computer Simulations and Results 

For part 1 of the assignment, I wrote a program called MNewMethAssign2. For this program, I set the values of (A_i, B_i, C_i) to equal (15600, 7540, 20140), (18760, 2750, 18610), (17610, 14630, 13480), and (19170, 610, 18390) km, respectively. The values of t_i equal 0.07074, 0.07220, 0.07690, and 0.07242 sec, respectively. Our initial vector, or initial guess, is (x₀, y₀, z₀, d₀) = (0, 0, 6370, 0). This program calls the Multivariate Newton’s Method function, which returns my solution. Then I print out the solution, and see what happens when I plug the solution into our F(x). When I run the code, I receive these results: 

A similar algorithm is shown below written in python. We use the Newton method to solve for the system of equations: 

Part 2 

We wrote a program that generates arbitrary (but consistent; i.e. the same arbitrary numbers are generated for each run) φ_is and θ_is, and plugged them into the definitions for A_i, B_i, C_i, t_i, defined above. When I ran the Multivariate Newton’s Method for these inputs, the program converged to a solution that seems to have an accuracy of at least 10⁻⁶. 

Then the code adds an input error, randomly adding or subtracting 10⁻⁸to the four ∆t_is. I set a for-loop with the following steps: apply the randomized ∆t_is (the input error), and run the Multivariate Newton’s Method. I define (∆x, ∆y, ∆z) as the difference between the values of (x, y, z) found by running the Multivariate Newton’s Method on a system with an input error and on a system without an input error. We then can calculate the forward error and the error magnification factor. 

While the for-loop repeats this process – applying randomized ∆t_is, running Multivariate Newton’s Method, and finding the errors – we select the maximum error magnification factor as our condition number, and we also select the maximum forward error. 

Similarly in python, we found each permutation of ∆t, by use of array q. Being that there are 4 satellites each with 2 possible errors, there are 2⁴ = 16 permutations of ∆t. We then found corresponding xyz coordinates for each permutation of ∆t. Then, we found the difference between the original coordinates, and each permutation of the error coordinates. Finally, we calculated the forward error of each permutation, the error magnification, and the max error magnification. 

Part 3 

We repeated part 2 with a slight modification. We set all the θ_is and φ_is to be within 5% of each other – meaning that we choose our first θ_i and φ_i arbitrarily, and then choose the subsequent θ_is and φ_is arbitrarily within 5% (on either side) of the first θ_i and φ_i chosen. the θ_is and φ_is are all within 5% of the first θ_i and φ_i on either side of the first θ_i and φ_i chosen. 

Without any input error, we received these vectors – which contain the solution and the result of the function when we input our solution, respectively. 

We also found the maximum forward error and the condition number (maximum error magnification factor): 

We observe that the condition number for loosely-bunched satellites (those of part 2) is much smaller than the condition number of the tightly-bunched satellites (those of part 3). The condition number for the loosely-bunched satellites has an order of magnitude of 10⁰, while the tightly-bunched satellites have an order of magnitude of 10³. 

We also observed several times that the Multivariate Newton’s Method did not converge in 1000 iterations for several runs for the tightly-bunched satellites with the input error. This occurred for 23/100 runs of the Multivariate Newton’s Method with an added input error. 

Similarly in python, we used the Newton equation to solve the system of 4 satellites with arbitrary phi and theta within 5% of each other. We followed the same steps as part 2, permutating the 16 possible dt’s, then calculating the corresponding positions and change in positions. Finally, we calculated the forward errors, found the maximum forward error, and the maximum error magnification (ConditionN). 

Numerical Methods Conclusion 

The Multivariate Newton’s Method is an effective way of solving an ill-conditioned system of equations. We solved this system for the input values specified in part 1, and for arbitrary input values. We added input errors, and identified the maximum forward error and the condition number of the system. We identified these values for both systems of tightly- and loosely- bunched satellites. We observed that when the satellites are tightly-bunched, the condition number grows significantly. 

Numerical Methods References 

[1] Marian Gidea. Assignment guidelines. Project-II-GPS 2021.pdf, 2021. [2] Marian Gidea. Systems of equations. Lecture 3 slides, pages 38–39, 2021. [3] Timothy Sauer. Numerical analysis, 3rd. edition. pages 248–251, 2017. 

Numerical Methods Sources

Yonah Moise did this project in Matlab, while Brian Livian did it in Python. Brian Livian wrote the Python codes. 

Brian Livian found the definitions of R_i and t_i, and corrected an error in generating the random φ_is. He also pointed out that we need to use arbitrary, as opposed to random numbers, in parts 2 and 3. 

Yonah Moise wrote the Matlab codes. They are named: SatFunct, MNewMethAssign2, part2, and part3. 

Yonah Moise wrote this assignment report, and Brian Livian contributed additional parts relevant to his Python codes, and the first few paragraphs of the introduction.

Numerical Methods

Physics Nobel Prize Winner MIT Prof Frank Wilczek on String Theory, Gravitation, Newton & Big Bang

Numerical Methods Numerical Methods Numerical Methods Numerical Methods Numerical Methods Numerical Methods Numerical Methods Numerical Methods