We Helped With This Mathematics Assignment: Have A Similar One?

Direct Methods and Least Squares

M426 - Peter Monk

Version of September 26, 2017

Introduction

In order to get maximum credit for this and future assignments, you must do as follows:

1.    This assignment should be completed using MATLAB.

2.    Hand in listings of each MATLAB script/function you use to obtain your results.

3.    Use format compact. Hand in output from MATLAB (you may cut and paste code, output and graphs to Word to create a report)

4.    Clear answers to all the questions.

5.    All code and written answers must be your own work. You may discuss the assignment but must submit your own code, results and writeup.

1             Work estimates for

Let us investigate if the work estimates for Gaussian elimination (matlab: ) correctly predict the time for computing the solution of linear systems in MATLAB. The idea is that Gaussian elimination for an n × n matrix uses (2/3)n³ + O(n²) operations. If each operation takes the same short time then the overall time for Gaussian elimination should be O(n³). This ignores the time needed to access memory and other things taking place on the computer that might slow it down.

For each n = 10,···,1000 create a random n × n matrix A and a random n × 1 vector b (use the command rand). Then using the matlab commands tic and toc (see the help page for tic) to measure the time taken to compute x = Ab and store the time in an array. This code may take a while to run!

You now have a vector of matrix sizes n = 10 : 1000 and corresponding vector of times t. According to theory we expect t = O(n³). Find the best fit to the data using the model t = a₁ +a₂n+a₃n² +a₄n³ (use the MATLAB polyfit command). Repeat with a quadratic model t = b₁ + b₂n + b₃n². Plot the two best fits and the actual time on one graph axis (use suitable axis labels and add a legend). Does your experiment support the theoretical estimate?

Hand in your graph with the comment on it and add printouts of the code.

2             Computing the determinant

If a matrix is poorly conditioned, the computed value of the determinant may not be very accurate. For example, set

U=round(100*rand(10))

U=triu(U,1)+0.1*eye(10)

In theory det(U) = det(U^T ) = 10⁻¹⁰ and det(UU^T ) = det(U)det(U^T ) = 10⁻²⁰. Compute det(U), det(U^T ) and det(UU^T ) using MATLAB. Also compute the condition number of U. Do the computed values of the determinant match the theoretical prediction?

3             Ill-conditioned matrix

This question is motivated by one in Numerical Computing with MATLAB by C. Moler. The question shows that a badly conditioned matrix does not necessarily lead to small pivots in Gaussian elimination (and can be a very simple matrix). The matrix is the n × n upper

triangular matrix A with elements





a_i,j =



Show how to generate this matrix in Matlab for any given n with eye, ones, and triu. Using the cond(A,1) command compute κ₁(A) for various n and plot a graph of the logarithm of the condition number against n (semilogy). For what n does κ₁(A) exceed 1/eps (this should be the final value of n for your graph)? Because this matrix is already upper triangular, Gaussian elimination with partial pivoting has no work to do. What are the pivots?

What is det(A)? This matrix is not singular, so Ax cannot be zero unless x is zero. Consider Ax = b and choose b = (0,0,···,1)^T . For n = 53 compute Ab. What is kxk and kbk (see the MATLAB norm command). This vector x is such that Ax is much smaller than x which suggests poor conditioning.

4             Least Squares

You are going to analyze some temperature data (in fact the average daily temperature in degrees Fahrenheit at McGuire AFB in New Jersey (the data is from NOAA). The data starts on Jan 1, 1955 and runs to Aug 13, 2010 (more is available online). There are two columns of data. The first is the date in the format YYYYMMDD. The second is the termperature.

Read in the data and graph it

1.    Download the file McGuireAFB.dat from Sakai.

2.    Use MATLAB’s “Import Data” (on the HOME menu) to read the data from the fileinto MATLAB. This will first show two column vectors of data. Make sure to change the column 1 header so this column is imported as TEXT (but column 2 is imported as NUMERIC).

3.    Use the datenum command to convert the dates in the first column to day numbers since the first day of data (which should be day zero). This should be stored in a variable denoted t. Store the temperature date in an array called temp. The following code will do the date conversion assuming your dates are stored in a cell array of text data called VarName1

day=datenum(VarName1,’yyyymmdd’); day=day-day(1)

4.    Plot a graph all the data (temperature against days from Jan 1, 1955) in a graph withsuitable annotation (i.e. axis labels, title),

5.    Plot two years worth of the data (I mean 730 days) again with suitable title and labels.The graph should show the rise and fall of temperature throughout the year

6.    Are there any missing days of data (how many)?

Preliminary analysis As you will have seen the data is quite erratic. To make sense of the data we can try to model the oscillatory fluctuations over a year by using a sine wave, plus an overall average temperature for the entire period. We are also interested in seeing if there is a trend (warming or cooling) in the data. Let t = 0,1,2,···,N (ie “days”) denote the days in the dataset. Since there are 365.25 days in a year the sinusoidal fluctuation due to seasonal change might be modeled by

f(t) = Asin(2πt/365.35 + φ)

where A is the amplitude of the fluctuation φ is the phase. This is inconvenient because the expression is a non-linear function of φ. Using trigonometry this is equivalent to

c₃ cos(2πt/365.25) + c₄ sin(2πt/365.25)

where c₃ = Asin(φ) and c₄ = Acos(φ). So our first model for the temperature is

f(t) = c₁ + c₂t + c₃ cos(2πt/365.25) + c₄ sin(2πt/365.25)

The coefficient c₁ (in ^◦F) will give the average temperature during the period and c₂ will be the trend (in units of ^◦F/day).

Find the best least squares fit of the temperature data using this model. Note that c₂ has units of degrees Fahrenheit per day. Convert this to degrees Fahrenheit per century and record also the converted quantity.

We are not sure if the fluctuation is exactly sinusoidal. Maybe higher order Fourier modes are present. Also fit the model f(t) = c₁+c₂t+c₃ cos(2πt/365.25)+c₄ sin(2πt/365.25)+c₅ cos(4πt/365.25)+c₆ sin(4πt/365.25)

Report the new coefficients. Do you think the higher Fourier modes help model the fluctuations?