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

$$\begin{aligned} A|b &= \begin{bmatrix} 0.1 & 0.3 & 0.9 & | & 1.3\\ 0.3 & 0.9 & 2.7 & | & 3.9\\ 0.6 & 0.7 & 0.1 & | & 1.4 \end{bmatrix} \\\ R_{2} - R_{1} \rightarrow A|b &= \begin{bmatrix} 0.1 & 0.3 & 0.9 & | & 1.3\\ 0.2 & 0.6 & 1.8 & | & 2.6\\ 0.6 & 0.7 & 0.1 & | & 1.4 \end{bmatrix} \\\ R_{3} - 3*R_{1} \rightarrow A|b &= \begin{bmatrix} 0.1 & 0.3 & 0.9 & | & 1.3\\ 0.2 & 0.6 & 1.8 & | & 2.6\\ 0.3 & -0.2 & -2.6 & | & -2.5 \end{bmatrix} \\\ R_3 - \frac{1}{2}*R_2 \rightarrow A|b &= \begin{bmatrix} 0.1 & 0.3 & 0.9 & | & 1.3\\ 0.2 & 0.6 & 1.8 & | & 2.6\\ 0.2 & -0.5 & -3.5 & | & -3.8 \end{bmatrix} \\\ \text{Thus } \rightarrow A|b \leftarrow &\begin{bmatrix} 0.1 & 0.3 & 0.9 & | & 1.3\\ 0.2 & 0.6 & 1.8 & | & 2.6\\ 0.2 & -0.5 & -3.5 & | & -3.8 \end{bmatrix} \\\ & \\\ det(A) = a(ei−fh)−b(di−fg)+c(dh−eg), A &=\begin{bmatrix} a & b & c \\ d & e & f\\ g & h & i \end{bmatrix} \\\ & \\\ \rightarrow det(A) = 0.1*(-0.6*3.5+1.8*0.5) - & \\\ 0.3*(-0.2*3.5-1.8*0.2) + & \\\ 0.9*(-0.5*0.2-0.6*0.2) &= 0 \end{aligned}$$

Lemma

$A$ is singular

b. With partial pivoting:

$$\begin{align} A|b &=\begin{bmatrix} 0.1 & 0.3 & 0.9 & | & 1.3\\ 0.3 & 0.9 & 2.7 & | & 3.9\\ 0.6 & 0.7 & 0.1 & | & 1.4 \end{bmatrix} \\\ R3 \leftrightarrow R1 \leftarrow A|b&=\begin{bmatrix} 0.6 & 0.7 & 0.1 & | & 1.4\\ 0.3 & 0.9 & 2.7 & | & 3.9\\ 0.1 & 0.3 & 0.9 & | & 1.3 \end{bmatrix} \\\ R2 - \frac{1}{2}R1 \leftarrow A|b&=\begin{bmatrix} 0.6 & 0.7 & 0.1 & | & 1.4\\ 0 & 0.55 & 2.65 & | & 3.2\\ 0.1 & 0.3 & 0.9 & | & 1.3 \end{bmatrix} \\\ R3 - \frac{1}{6}R1 \leftarrow A|b&=\begin{bmatrix} 0.6 & 0.7 & 0.1 & | & 1.4\\ 0 & 0.55 & 2.65 & | & 3.2\\ 0 & 0.18333333 & 0.88333333 & | & 1.06666667 \end{bmatrix} \\\ R3 - \frac{1}{3}R2 \leftarrow A|b&=\begin{bmatrix} 0.6 & 0.7 & 0.1 & | & 1.4\\ 0 & 0.55 & 2.65 & | & 3.2\\ 0 & 0 & 0 & | & -0.3 \end{bmatrix} \end{align}$$

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

$$\begin{aligned} R_3 \leftrightarrow R_1 \rightarrow U &= \begin{bmatrix} 0.6 & 0.7 & 0.1\\ 0.3 & 0.9 & 2.7\\ 0.1 & 0.3 & 0.9 \end{bmatrix}, \quad P_1 = \begin{bmatrix} 0 & 0 & 1\\ 0 & 1 & 0\\ 1 & 0 & 0 \end{bmatrix} \\\ R_2 - \frac{1}{2}R_1 \rightarrow U &= \begin{bmatrix} 0.6 & 0.7 & 0.1\\ 0 & 0.55 & 2.6500000000000004\\ 0.1 & 0.3 & 0.9 \end{bmatrix} \\\ R_3 - \frac{1}{6}R_1 \rightarrow U &= \begin{bmatrix} 0.6 & 0.7 & 0.1\\ 0 & 0.55 & 2.6500000000000004\\ 0 & 0.18333333333333335 & 0.8833333333333333 \end{bmatrix} \\\ R_3 - \frac{1}{3}R_2 \rightarrow U &= \begin{bmatrix} 0.6 & 0.7 & 0.1\\ 0 & 0.55 & 2.6500000000000004\\ 0 & 0 & 4.8109664400423476 \times 10^{-17} \end{bmatrix} \end{aligned}$$

_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}$

$R5+R4\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 & 0 & 16 \end{bmatrix}$

Problem 3 [5 points] (a) (3 points) Let $A$, $B$, and $C$ be $n × n$ matrices, where $B$ and $C$ are nonsingular. For an $n-$vector $b$, describe how you would implement the formula $x = C^{-1} (A + I)(A + B^{−1})b.$ without computing any inverses. Here, $I$ is the $n × n$ identity matrix.

(b) (2 points) What is the complexity of your approach in terms of big-O notation?

Answer:

a. Given $B$ and $C$ are non-singular

1. Step 1: Using $LU$ decomposition of B, such that $B=LU$
2. Step 2: Solve for $y$ in $By=b$ (As $y=B^{-1}b$)
   1. solve for $z$ in $Lz=b$ via forward substitution
   2. solve for $y$ in $Uy=z$ via backward substitution
3. Step 3: Compute $z=(A+B^{-1})b$
   1. This becomes $z=Ab+y$
4. Step 4: Compute $w = (A+I)z$
   1. Via matrix multiplication $\rightarrow w=Az + z$
5. Step 5: using $LU$ decomposition of C, such that $C=LU$
6. Step 6: Solve for $x$ in $Cx=w$ (As $x=C^{-1}w$)
   1. Solve for $z'$ in $Lz'=w$ via forward substitution
   2. Solve for $x$ in $Ux=z'$ via backward substitution

With expansion, solved $x = C^{-1} (A + I)(A + B^{−1})b.$

b. Complexity analysis

Let total_cost be the big-O notation

Step 1 using $LU$ decomposition of $B$ $\rightarrow \text{total\_cost}=O(n^3)$

Step 2 solving each $Lz=b$ and $Uy=z$ takes $O(n^2)$ each, thus solving $Lz=b$ using $LU$ decomposition takes $O(2n^2)$ $\rightarrow \text{total\_cost}=O(n^3) + O(2n^2)$

Step 3 Compute $z=(A+B^{-1})b$

• MatmulOp of $Ab$ is $O(n^2)$
• AddOp of $Ab+y$ is $O(n)$
• Total for this step $O(n^2) + O(n)$ $\rightarrow \text{total\_cost}=O(n^3) + O(3n^2) + O(n)$

Step 4 Compute $w = (A+I)z$

• MatmulOp of $Ab$ is $O(n^2)$
• AddOp of $Ab+y$ is $O(n)$
• Total for this step $O(n^2) + O(n)$ $\rightarrow \text{total\_cost}=O(n^3) + O(4n^2) + O(2n)$

Step 5 using $LU$ decomposition of $C$ $\rightarrow \text{total\_cost}=O(2n^3) + O(4n^2) + O(2n)$

Step 6 solving each $Lz'=w$ and $Ux=z'$ using LU composition takes $O(2n^2)$ $\rightarrow \text{total\_cost}=O(2n^3) + O(6n^2) + O(2n)$

---

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:

hilb_problem.m

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}$

n	$\frac{\Vert{r}\Vert}{\Vert{b}\Vert}$	$\frac{\Vert{\hat{x} - x}\Vert}{\Vert{x}\Vert}$
1	0.00000000000000000000000000000000	0.00000000000000000000000000000000
2	0.00000000000000000000000000000000	0.00000000000000013220372219891702
3	0.00000000000000000000000000000000	0.00000000000000363350625815651572
4	0.00000000000000000000000000000000	0.00000000000006709266750580992637
5	0.00000000000000007733975117624287	0.00000000000747821082933078000054
6	0.00000000000000013934207506736382	0.00000000023960543432895825359428
7	0.00000000000000010660570398371085	0.00000000837749558262967895463873
8	0.00000000000000007165565184570407	0.00000009992506975169996005028294
9	0.00000000000000007076549838447114	0.00000608952488692639798140973303
10	0.00000000000000012662840530707719	0.00002450986238666613242472361311
11	0.00000000000000011997633780813789	0.00379971054180424641297242338567
12	0.00000000000000006503338066505365	0.25404291536273732043937911839748

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

x	1	4	9
$\sqrt{x}$	1	2	3

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