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
ε(θ) = m1sin(θ) + m2cos(θ) + b
and solve for m and φ in terms of m1 and m2.
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 m1, m2, 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||| = maxi ∑j|aij|
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?