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.)

© Copyright 2019 ExploreDoc