Assignment 3, due 18 Nov 2014

Scientific Computing for DPhil Students I — Assignment 3
Due at 9:00 on Tuesday, 18 November 2014. This is the third of four assignments this term.
1.
Fitting ellipses via least-squares. Suppose we have n data points (x1 , y1 ), . . . , (xn , yn ) in the plane
and we want to find an ellipse that fits them well. Finding the geometrically closest fit is a nonlinear
problem, but we can come close by a linear formulation. An equation for an ellipse centred at (0, 0) is
bx2 + cxy + dy 2 = 1.
Let us view b, c and d as unknowns and find them by solving a linear least-squares problem.
(a) Write down in matrix form an n × 3 least-squares problem whose unknown vector is (b, c, d)T .
(b) Write a Matlab function [b,c,d] = ellipse(x,y) which uses “ \ ” to solve this problem. Write
driver code to call ellipse for the data
(3, 3), (1, −2), (0, 3), (−1, 2), (−2, −2), (0, −4), (−2, 0), (2, 0).
and then print b, c, d and also plot the data points and the fitting ellipse. (Hint: to plot the ellipse you
may find it helpful to take a range of angles θ and work with corresponding ratios y/x = tan θ.)
(c) Here’s a little code to let you put in points interactively with the mouse. Try it for some data points
of your own choosing and turn in the resulting plot of a fitted ellipse.
hold off, axis([-3 3 -3 3]), axis manual, hold on, grid on
x = []; y = []; button = 1;
disp(’input points with mouse, button >= 2 for final point’)
while button == 1
[xx,yy,button] = ginput(1)
x = [x; xx]; y = [y; yy]; plot(xx,yy,’x’)
end
2.
Least-squares fitting in Chebfun. Read Chapters 1 and 6 of the Chebfun Guide (www.chebfun.org).
Download Chebfun and put it in your Matlab path. Using Chebfun, Compute the degree 5 polynomial
p that is closest to f (x) = tanh(10x − 5) in the L2 norm on [−1, 1]. Plot the error and report the 2and ∞-norms. Now use the Chebfun command remez to compute the degree 5 polynomial p∗ that is
closest to f in the L∞ norm. Again plot the error and report the two norms.
3.
Eigenvalues of the Laplacian on a trapezoid. Consider the trapezoidal region in the x-y plane bounded
by the points (0, 0), (1, 0), (1, 1), and (−1, 1), and consider the Laplace eigenvalue problem
(1)
∂2u ∂2u
+ 2 = −λ2 u
∂x2
∂y
with u = 0 along the boundary. Our mission is to find the lowest three eigenvalues of this problem
using least-squares fitting. The idea is this. For any λ > 0, n, and coefficients cj , the sum
(2)
w(r, θ) = c1 J4/3 (λr) sin(4θ/3) + c2 J8/3 (λr) sin(8θ/3) + · · · + cn J4n/3 (λr) sin(4nθ/3)
satisfies (1) except for the boundary condition along the two boundary segments running from (1, 0)
to (1, 1) and (−1, 1) to (1, 1). (Here θ is the angle up from the x-axis and Jν is a Bessel function; in
Matlab, besselj.) If we can find λ, n and cj such that w is (nearly) zero on this boundary, then we
(nearly) have an eigenvalue and eigenvector.
Given n, discretize the two boundary segments by plenty of points, say, m = 3n. Assuming λ is given,
describe the m × n system of equations Ax = 0 we would like a vector of coefficients c to satisfy.
Explain why we cannot expect a solution in general but why it makes sense more generally to take c to
be the singular vector σmin (column of V ) associated with the minimum singular value of A. Plot σmin
against λ for various values of n and a suitable range of λ. What’s your best estimate of the first three
eigenvalues? (You may use Matlab’s fminbnd to help with this, or other methods.)

Download Report