--- date: '2023-10-24' description: gaussian elimination with partial pivoting, lu decomposition, polynomial interpolation using lagrange and newton basis, and hilbert matrix conditioning. id: A2 modified: 2026-06-06 00:10:38 GMT-04:00 tags: - swfr4x03 title: Gaussian elimination, LU decompositions, and errors LS solving created: '2023-10-24' published: '2023-10-24' pageLayout: default slug: thoughts/university/twenty-three-twenty-four/compsci-4x03/A2 permalink: https://aarnphm.xyz/thoughts/university/twenty-three-twenty-four/compsci-4x03/A2.md generator: quartz: v4.6.0 hostedProvider: Cloudflare baseUrl: aarnphm.xyz full: https://aarnphm.xyz/llms-full.txt --- **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} $$ > \[!tip\] 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: ```matlab title="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._ > \[!tip\] accuracy > > 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}$ > \[!tip\] results > > 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: ```matlab 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$ ```matlab % (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'); ``` ![[thoughts/university/twenty-three-twenty-four/compsci-4x03/a2-fig1.webp]] ![[thoughts/university/twenty-three-twenty-four/compsci-4x03/a2-fig2.webp]] b. The following is a snippet of `interp_problem.m` for Cheybyshev points ```matlab % (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'); ``` ![[thoughts/university/twenty-three-twenty-four/compsci-4x03/a2-fig3.webp]] ![[thoughts/university/twenty-three-twenty-four/compsci-4x03/a2-fig4.webp]] c. The following is a snippet of `interp_problem.m` through spline interpolation at $n + 1$ equally spaced points. ```matlab % (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 [[thoughts/representations|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\] 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\] 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\] 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.