Gaussian elimination, LU decompositions, and errors LS solving

Problem 1 [8 points] Consider the system $Ax = b$, where $A=\begin{bmatrix} 0.1 & 0.3 & 0.9\\ 0.3 & 0.9 & 2.7\\ 0.6 & 0.7 & 0.1 \end{bmatrix}$ and $b = \begin{bmatrix} 1.3 & 3.9 & 1.4\end{bmatrix}^T$

a. [2 points] Show that $A$ is singular.

b. [2 points] If we were to use Gaussian elimination with partial pivoting to solve this system using exact arithmetic, show where the process fails.

c. [2 points] Solve this system in double precision using partial pivoting. Do not use Matlab’s functions. What is the solution that you obtain?

d. [2 points] Matlab’s A\b produces NaN -Inf Inf as a solution. Explain why NaN, -Inf and Inf.

Answer:

a. For $A$ to be singular, prove $det(A) = 0$

Using Gaussian elimination without partial pivoting

Lemma

$A$ is singular

b. With partial pivoting:

We notice that $R3-\frac{1}{3}R2 \rightarrow 0=-0.3$, thus invalid.

c. With partial pivoting in double precision

The $LU$ decomposition of $A=\begin{bmatrix} 0.1 & 0.3 & 0.9\\ 0.3 & 0.9 & 2.7\\ 0.6 & 0.7 & 0.1 \end{bmatrix}$

The following portray steps to calculate $U$ (lower triangular):

_note: the $a_{33}$ is close to zero, hence consistent with previous finding_

$L=\begin{bmatrix} 1 & 0 & 0\\ 0.5 & 1 & 0\\ 0.16666666666666669 & 0.33333333333333326 & 1 \end{bmatrix}$

To solve for $x$ with $LU$ decomposition, We solve $L(Ux)=Pb$

$\rightarrow x=\begin{bmatrix} 14.006993006993 & -10.48951048951048 & 3.3846153846153832\end{bmatrix}$

d. Since A is singular, it doesn’t have an inverse. Matlab uses LU decomposition, and as we explored above, a pivot element is found to be zero or close to zero (matrix is probably ill-conditioned), which leads to $0x_1 + 0x_2 + 0x_3=\text{non negative value}$, which results in NaN. For the second value -Inf, the division is small. Inf is due to division by zero

Problem 2 [2 points] Apply Gaussian elimination with partial pivoting on the following matrix $A=\begin{bmatrix} 1 & 0 & 0 & 0 & 1\\ −1 & 1 & 0 & 0 & 1\\ −1 & −1 & 1 & 0 & 1\\ −1 & −1 & −1 & 1 & 1\\ −1 & −1 & −1 & −1 & 1 \end{bmatrix}$

Show all the steps.

Answer:

$A=\begin{bmatrix} 1 & 0 & 0 & 0 & 1\\ −1 & 1 & 0 & 0 & 1\\ −1 & −1 & 1 & 0 & 1\\ −1 & −1 & −1 & 1 & 1\\ −1 & −1 & −1 & −1 & 1 \end{bmatrix}$

$R2+R1 \text{ and } R3+R1\text{ and } R4+R1\text{ and } R5+R1\rightarrow A=\begin{bmatrix} 1 & 0 & 0 & 0 & 1\\ 0 & 1 & 0 & 0 & 2\\ 0 & −1 & 1 & 0 & 2\\ 0 & −1 & −1 & 1 & 2\\ 0 & −1 & −1 & −1 & 2 \end{bmatrix}$

$R3+R2 \text{ and } R4+R2\text{ and } R5+R2\rightarrow A=\begin{bmatrix} 1 & 0 & 0 & 0 & 1\\ 0 & 1 & 0 & 0 & 2\\ 0 & 0 & 1 & 0 & 4\\ 0 & 0 & −1 & 1 & 4\\ 0 & 0 & −1 & −1 & 4 \end{bmatrix}$

$R4+R3 \text{ and } R5+R3\rightarrow A=\begin{bmatrix} 1 & 0 & 0 & 0 & 1\\ 0 & 1 & 0 & 0 & 2\\ 0 & 0 & 1 & 0 & 4\\ 0 & 0 & 0 & 1 & 8\\ 0 & 0 & 0 & −1 & 8 \end{bmatrix}$

Problem 4 [6 points] An $n × n$ Hilbert matrix, denote it by $H$, has entries $h_{ij} = \frac{1}{(i+j-1)}, i, j = 1, . . . , n.$

For $n = 2, 3, . . .$ , generate the Hilbert matrix of order $n$, and also generate the $n-$vector $b = Hx$, where $x$ is a random vector. Solve the resulting system $Hx = b$ to obtain an approximate solution $\hat{x}$. (See the functions hilb and rand.)

(a) [2 points] How large can you take $n$ before the error $\frac{\Vert{\hat{x} - x}\Vert}{\Vert{x}\Vert}$ is 100 percent?

(b) [2 points] For $n$ up to the value you find in (a), report $\frac{\Vert{r}\Vert}{\Vert{b}\Vert}$ , where $r = b − H\hat{x}$, and $\frac{\Vert{\hat{x} - x}\Vert}{\Vert{x}\Vert}$.

(c) [2 points] As $n$ increases, how does the number of correct digits in the computed solution relate to the condition number of the matrix? See the cond function.

Submit your Matlab program producing the above results. Name the Matlab file hilb_problem.m.

Answer:

The following hilb_problem.m is used:

function hilb_problem()
    n = 1;
    while true
        % Generate Hilbert matrix of order n
        H = hilb(n);
 
        % Generate random vector x
        x = rand(n, 1);
 
        % Compute b = Hx
        b = H * x;
 
        % Solve the system Hx = b
        x_hat = H \ b;
 
        % Compute the relative error
        error = norm(x_hat - x) / norm(x);
        fprintf("error=%d, n=%d\n", error, n)
        % If the error is 100 percent, break
        if error >= 1
            break;
        end
 
        n = n + 1;
    end
 
    fprintf('\n=============\n\nThe largest n before the error is 100 percent is: %d\n\n=============\n', n-1);
 
    for i = 1:n-1
        H = hilb(i);
        x = rand(i, 1);
        b = H * x;
        x_hat = H \ b;
 
        r = b - H * x_hat;
        rel_resid = norm(r) / norm(b);
        rel_error = norm(x_hat - x) / norm(x);
 
        %fprintf('%d %.16f\n',i, rel_resid)
        fprintf('| %d | %.32f | %.32f |\n', i, rel_resid, rel_error);
    end
 
    cond_num = cond(H);
    fprintf('The condition number of the matrix for n = %d is: %f\n', n-1, cond_num);
end

a. largest $n=12$ before the error $\frac{\Vert{\hat{x} - x}\Vert}{\Vert{x}\Vert}$ is 100 percent.

b. The following entails the value of $\frac{\Vert{r}\Vert}{\Vert{b}\Vert}$ and $\frac{\Vert{\hat{x} - x}\Vert}{\Vert{x}\Vert}$

c. As $n$ increases, the condition number increases, which means the matrix becomes more ill-conditioned. This means fewer digits in the computed solution are correct.

IMPORTANT

The number of correct digits in the computed solution decreases due to the increase in the condition number as $n$ increases

Problem 5 [4 points] You have to interpolate $sin(x)$ by a polynomial of degree five using equally spaced points in [0, 1].

(a) [2 points] What (absolute) error would you expect if you use this polynomial?

(b) [2 points] Using equally spaced points, what degree polynomial would you use to achieve a maximum error of $10^{-8}$?

Answer:

a. Interpolate $sin(x)$ by a polynomial of degree five using equally spaced on in $[0,1]$, Error as follow

$f(x) - p_n(x) = E(x) = \frac{f^{n+1}(\xi)}{(n+1)!}\prod_{i=0}^{n}{(x-x_i)}$

where

• $n$ is the degree of the polynomial ($n=5$)
• $f^{n+1}(\xi)$ is $(n+1)\text{-th}$ derivate of $f$

Derivate of $sin(x)$ every 4 terms is $sin(x), cos(x), -sin(x), -cos(x)$. Therefore the 6th derivative is $-cos(x)$

Here $h=\frac{b-a}{n}=\frac{1}{5}$ and $M = max_{0\leq t\leq 1}|-cos(t)| = 1 - cos(1) = 2sin^2(\frac{1}{2})$

Therefore $|E(x)| = |f(x) - sin(x)| \leq \frac{M}{4(n+1)}h^{n+1}=\frac{2sin^2(\frac{1}{2})}{4(6)}(1/5)^6 \approx 1.225860517684960×10^{−6}$

b. To achieve maximum error of $10^{-8}$, We have

$|f(x) - sin(x)| \leq\frac{max_{0\leq t\leq 1}|sin^{(n+1)}(t)|}{4(n+1)*n^{n+1}} = 10^{-8}$

derivative of $sin(x)$ cycles every 4 term, thus the max value of $|sin^{(n+1)}(t)|$ over $[0,1]$ is 1

Thus we need to solve for $n$ in $\frac{1}{4(n+1)n^{n+1}}=10^{-8} \rightarrow n\approx 7 \text{ (through trial and error)}$

Hence considering to use polynomial degree seven to achieve the desired error bound.

Problem 6 [3 points] You are given the values of $\sqrt{x}$ at three points

(a) [2 points] Construct the interpolating polynomial interpolating these data.

(b) [1 points] Using this polynomial, approximate $\sqrt{1.5}$.

Answer:

a. To construct the interpolating polynomial for these data, we will use Lagrange basis

$P(x)=\sum_{i=0}^{n-1}{y_i}{L_i(x)}$

where $L_i(x)$ is the $i\text{-th}$ Lagrange basis polynomial, defined as

$L_i(x) = \prod_{j=0,j\neq i}^{n-1}\frac{x-x_j}{x_i-x_j}$

With $y(x) = \sqrt{x}$, and data point $x_0=1,y_0=1;x_1=4,y_1=2;x_2=9,y_2=3$

$P(x)=\sum_{i=0}^{2}{y_i}{L_i(x)} \text{ where } L_i(x) = \prod_{j=0,j\neq i}^{2}\frac{x-x_j}{x_i-x_j}$

$L_0(x) = \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} = \frac{(x-4)(x-9)}{(1-4)(1-9)} = \frac{(x-4)(x-9)}{24}$ $L_1(x) = \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}=\frac{(x-1)(x-9)}{(4-1)(4-9)}=\frac{(x-1)(9-x)}{15}$ $L_2(x) = \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}=\frac{(x-1)(x-4)}{(9-1)(4-1)} = \frac{(x-4)(x-1)}{40}$

$P(x) = y_0L_0(x) + y_1L_1(x) + y_2L_2(x) = 1 * \frac{(x-4)(x-9)}{24} + 2*\frac{(x-1)(9-x)}{15} + 3*\frac{(x-4)(x-1)}{40}$

IMPORTANT

The interpolating polynomial $P(x)=\frac{(x-4)(x-9)}{24} + \frac{2(x-1)(9-x)}{15} + \frac{3(x-4)(x-1)}{40}$

b. The approximation of $P(\sqrt{1.5})=\frac{(1.5-4)(1.5-9)}{24} + \frac{2(1.5-1)(9-1.5)}{15} + \frac{3(1.5-4)(1.5-1)}{40}=1.1875$

Problem 7 [7 points] Let $f(x) = \frac{sin(x)}{(1+20x)^2}$. Interpolate this function over $x \in [−1, 1]$ using

(a) [2 points] polynomial interpolation of degree $n = 15$ at equally spaced points. Then evaluate this polynomial at $N = 100$ equally spaced points. Denote the interpolating polynomial by $p(x)$. Plot

• $f(x)$ and $p(x)$ versus $x$ at the interpolation points and at the $N$ points (on the same plot);
• $|f(x) − p(x)|$ versus $x$ at the $N$ points. You can use the polyfit function. See the linspace function.

(b) [2 points] Repeat (a) but now using Chebyshev points.

(c) [2 points] Repeat (a) but now using spline interpolation at $n + 1$ equally spaced points. See the spline function.

(d) [1 points] Discuss the accuracies of your results.

Submit your plots (6 in total) and the Matlab code producing them. Name your Matlab file interp_problem.m.

Answer

$f(x)$ implementation in matlab are as follow:

f = @(x) sin(x)./((1 + 20*x).^2);

a. The following is a snippet of interp_problem.m for polynomial interpolation of degree $n=15$

% (a) Polynomial interpolation of degree n = 15 at equally spaced points
 
% Define the number of interpolation points and the degree of the polynomial
n = 15;
N = 100;
 
% Generate n+1 equally spaced points in the interval [-1, 1]
x = linspace(-1, 1, n+1);
y = f(x);
 
% Interpolate using polyfit
p_coeff = polyfit(x, y, n);
 
% Evaluate the interpolating polynomial at N equally spaced points
x_N = linspace(-1, 1, N);
p_N = polyval(p_coeff, x_N);
 
% Plot f(x) and p(x) on the same graph
figure;
plot(x_N, f(x_N), 'b-', x_N, p_N, 'r--', x, y, 'go');
legend('f(x)', 'p(x)', 'Interpolation Points');
title('f(x) and p(x) vs. x');
xlabel('x');
ylabel('y');
 
% Plot the absolute error |f(x) - p(x)| at the N points
figure;
plot(x_N, abs(f(x_N) - p_N), 'm-');
title('Absolute Error |f(x) - p(x)| vs. x');
xlabel('x');
ylabel('Error');

b. The following is a snippet of interp_problem.m for Cheybyshev points

% (b) Polynomial interpolation using Chebyshev points
 
% Generate Chebyshev points in the interval [-1, 1]
x_cheb = cos((2*(1:n+1)-1)*pi/(2*n));
y_cheb = f(x_cheb);
 
% Interpolate using polyfit
p_cheb_coeff = polyfit(x_cheb, y_cheb, n);
 
% Evaluate the interpolating polynomial at N equally spaced points
p_cheb_N = polyval(p_cheb_coeff, x_N);
 
% Plot f(x) and p(x) using Chebyshev points on the same graph
figure;
plot(x_N, f(x_N), 'b-', x_N, p_cheb_N, 'r--', x_cheb, y_cheb, 'go');
legend('f(x)', 'p(x) with Chebyshev', 'Interpolation Points');
title('f(x) and p(x) with Chebyshev vs. x');
xlabel('x');
ylabel('y');
 
% Plot the absolute error |f(x) - p(x)| using Chebyshev points at the N points
figure;
plot(x_N, abs(f(x_N) - p_cheb_N), 'm-');
title('Absolute Error |f(x) - p(x) with Chebyshev| vs. x');
xlabel('x');
ylabel('Error');

c. The following is a snippet of interp_problem.m through spline interpolation at $n + 1$ equally spaced points.

% (c) Spline interpolation at n+1 equally spaced points
 
% Evaluate the function at n+1 equally spaced points
y_spline = f(x);
 
% Use the spline function to get the piecewise polynomial representation
pp = spline(x, y_spline);
 
% Evaluate the spline at N equally spaced points
spline_N = ppval(pp, x_N);
 
% Plot f(x) and the spline on the same graph
figure;
plot(x_N, f(x_N), 'b-', x_N, spline_N, 'r--', x, y_spline, 'go');
legend('f(x)', 'spline(x)', 'Interpolation Points');
title('f(x) and spline(x) vs. x');
xlabel('x');
ylabel('y');
 
% Plot the absolute error |f(x) - spline(x)| at the N points
figure;
plot(x_N, abs(f(x_N) - spline_N), 'm-');
title('Absolute Error |f(x) - spline(x)| vs. x');
xlabel('x');
ylabel('Error');

d. Discussion

1. The polynomial interpolation using equally spaced points might show oscillations near endpoints due to Runge phenomenon (oscillations near the endpoints of the interpolated interval become pronounced). We saw oscillation in the error graph here.
2. Polynomial interpolation using Chebyshev points should mitigate the oscillations
3. The spline interpolation will provide a piecewise polynomial that should fit the function smoothly and might offer better accuracy than polynomial interpolation

Problem 8 [4 points] Given the three data points $(−1, 1), (0, 0), (1, 1)$, determine the interpolating polynomial of degree two using:

a. [1 point] monomial basis

b. [1 point] Lagrange basis

c. [1 point] Newton basis

[1 point] Show that the three representations give the same polynomial.

Answer:

a. Monomial basis

The monomial basis for a polynomial of degree two is given by: $p(x)=a_0+a_1*x+a_2*x^2$

The linear system as follow

$a_0-a_1+a_2=1$ $a_0=0$ $a_0+a_1+a_2=1$

Solving this system to obtain the $a_0=0,a_1=0, a_2=1$

NOTE

Thus monomial basis of this polynomial of degree two is $p(x) = x^2$

b. Lagrange basis

The Lagrange basis for a polynomial of degree two is given by: $p(x)=\sum_{j=0}^{2}{y_j}{L_j(x)} = f(x_0)L_0{(x)} + f(x_1)L_1{(x)} + f(x_2)L_2{(x)}$ where $L_0(x) = \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} = \frac{x(x-1)}{2}$ $L_1(x) = \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}=-x(x-1)$ $L_2(x) = \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}=\frac{x(x+1)}{2}$

Thus $p(x) = 1*\frac{x(x-1)}{2} + 0*(-x(x-1)) + \frac{x(x+1)}{2} = x^2$

NOTE

Thus Lagrange basis of this polynomial of degree two is $p(x) = x^2$

c. Newton basis

The Newton basis for a polynomial of degree two is given by: $p(x)=f(x_0)+(x-x_0)f[x_0, x_1] + (x-x_0)(x-x_1)f[x_0, x_1, x_2]$ where $f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0} = \frac{0-1}{0+1} = -1$ $f[x_0,x_1,x_2]=\frac{f[x_1, x_2]-f[x_0, x_1]}{x_2-x_0} = \frac{1+1}{1+1} = 1$

We have $f[x_1,x_2]=\frac{f(x_2)-f(x_1)}{x_2-x_1} = \frac{1-0}{1-0} = 1$

Thus $p(x)=1+(x+1)(−1)+(x+1)(x)*2 =1 - x-1 + (x^2+x)=x^2$

NOTE

Thus Newton basis of this polynomial of degree two is $p(x) = x^2$

Therefore, we prove that all three basis yield the same polynomial for degree two.