##### Reference no: EM131097701

**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_{1}sin(θ) + m_{2}cos(θ) + b

and solve for m and φ in terms of m_{1} and m_{2}.

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_{1}, m_{2}, 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 δ_{1} = (ε, ε) and δ_{2} = (ε, -ε). 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