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 $\mathrm{\mid }E(x)\mathrm{\mid }=\mathrm{\mid }f(x)-sin(x)\mathrm{\mid }\le \frac{M}{4(n+<}$