E19: Numerical Methods for Engineering Applications Spring 2016 - HOMEWORK 4

1. Compass calibration

In the presence of magnets or circuitry, a magnetic compass may be subject to heading dependent heading error of the form

Ε(θ) = m sin(θ + φ) + b                                                                                   (1)

where m is the amplitude of the heading error, φ is the phase, and b is the DC offset. Given the true heading θ_T the measured heading θ_M is given by

θ_M = θ_T + ε(θ_T)                                                                                                    (2)

Given a set of true and measured compass data, we wish to estimate the unknown parameters m, φ, and b of the error function, and to be able to compute the true heading from the measured heading.

a. Make the error linear in the unknowns. Show that the error function (1) can also be written as

ε(θ) = m₁sin(θ) + m₂cos(θ) + b

and solve for m and φ in terms of m₁ and m₂.

b. Fit the data from the zip file. The file compass_data.txt is a simple text file with each line providing a (θ_T , θ_M) pair. You can load it into Python as an array using the command

data = numpy.genfromtxt('compass_data.txt')

Set up and solve a least-squares problem to obtain m₁, m₂, and b, and use your answer above to determine m and φ as well.

c. Create a function to estimate true heading. Equation (2) above tells how to compute θ_M from θ_T but we typically want to do the opposite. Use the fixed point iteration method to write a function which takes as input θ_M and then computes a number x such that x + ε(x) = θ_M. Verify your function by checking its output for 1.5, 2.0, and 2.5 radians.

2. Iterative solution methods

Modify the heat_1d.py example from class to solve the system using either Jacobi iteration or Gauss-Seidel iteration, but with a twist: you should not explicitly represent the A matrix or b vector at all. Since the structure of the A matrix in this problem is so regular, you should be able to replace any reference to a diagonal element of A with α, reference to to any element one step from the diagonal with -1, and totally remove any references to entries of A more than one step from the diagonal. You can reason similarly about the elements of b.

Your solver should terminate when the approximate relative error

ε_a = ||x^(k+1) - x^(k)||/||x^(k+1)||

falls below 1e-5 (in my experience, this takes on the order of 150 iterations for n = 50).

3. Condition number of a matrix

Note: we will finish covering this material on 2/15.

For a small number ε > 0, the matrix

is nearly un-invertible. How much so? We can find out using the condition number of A.

Using the matrix norm

|||A||| = max_i ∑_j|a_ij|

and the formula for the inverse of a 2 x 2 matrix, compute

cond(A) = |||A|||·|||A^-1|||

analytically, in terms of ε.

Next, let's examine this in practice. For ε = 10^-6, numerically compute the solution to Ax = b for b = (1, 1). Now observe how the solution changes as you perturb b with each of two small perturbations of δ₁ = (ε, ε) and δ₂ = (ε, -ε). Explain what you observe in terms of the matrix condition number. Recall that the condition number was originally defined as an inequality - how tight a bound is it in these two cases?

Attachment:- Assignment.rar