LE/CIVL 3140
Assignment 3
Due : Wednesday, November 17, 2021
y 11:59 p.m.
Instructor: KD Papoulia
1. Given data points (x1, y1), . . . , (xn, yn) with x1 < x2 < · · · < xn, a differentiable
piecewise quadratic interpolant is a function q(x) defined so that q(x) = qi(x) on
the interval [xi, xi+1] and
qi(x) = ai(x− xi)2 + bi(x− xi) + ci
that interpolates the data points. For all except the first interval, the coefficients
(ai, bi, ci, i = 2, . . . , n are found from
qi(xi) = yi
qi(xi+1) = yi+1
q′i(xi) = q
′
i−1(xi)
These are not enough equations to uniquely specify q1, . . . , qn−1, i.e., specify the
coefficients (a1, b1, c1), . . . , (an−1, bn−1, cn−1. We need one more piece of data, so we
assume the user also specifies the value of q′(x1) = q
′
1(x1).
A differentiable (called C1 for continuous first derivative) piecewise quadratic, is
not widely used in practice because it has some pathological properties that will be
exhibited in the following questions1. The m-file called quadinterp from the e-class
website ca
ies out the C1 quadratic interpolation and plots the interpolant as well
as the original data points.
(a) Derive the formulas used in the m-file on lines 26–28 and line 30 to compute the
interpolant.
(b) Test the function on the data points (xi, cos(xi)), where x1, . . . , x20 are 20 evenly
spaced points in the interval [0, π]. Use the co
ect value of the first derivative of
the cosine function at x1 as the third argument to quadinterp.
1A better behaved alternative is a C0 quadratic, which does not require continuous derivatives but
interpolates three data points per segment. If C1 interpolation is needed, cubic splines are a better choice.
Cubic Hermite interpolation, covered in lecture, is also C1 but requires given values of the derivative at
all data points.
1
(c) Now try again with the same choice of (x1, y1), . . . , (x20, y20) with an inaccurate
specification of the third argument, namely, 1. Do you observe oscillations and how
do they behave as x increases?
(d) Here is a particularly pathological case: the x-coordinates are 20 points loga-
ithmically spaced in [10, 100] (see help on logspace), the y-coordinates are all 0
(so that that the true function is identically 0), and the initial slope is inaccurately
set to a nonzero value. Provide an explanation of the oscillatory plot that results.
For (b)–(d), hand in the statements that you typed in the command window or put
into a script to generate the requested plots, and also hand in your plots. Hand in
written answers to the questions posed in (a), (c), and (d).
2. Recall that a Vandermonde matrix is an n × n matrix formed from a vector x⃗ =
(x0, x1, x2, . . . , xn−1) as follows:
V (w⃗) =
xn−10
xn−11
xn−12
...
xn−1n−1
xn−20
xn−21
xn−22
...
xn−2n−1
· · ·
x20
x21
x22
...
x2n−1
x0
x1
x2
...
xn−1
1
1
1
...
1
. (1)
You may generate Vandermonde matrices using the MATLAB command vander,
which takes the vector x⃗ as its single parameter.
(a) Consider the Vandermonde matrix generated with x = [10.0 : 0.1 : 11.0]. Ob-
tain its LU decomposition and solve the system V (x)a = y, where
y = [1, 9, 1, 9, 1, 9, 1, 9, 1, 9, 1]T
and a are the coefficients of the polynomial interpolating data (xi, yi), using
the built-in Matlab commands [L,U]=lu(V) and V\ y. Compare LU and V
in some suitable norm and comment on the result; also compare Va, where
a=V\ y, and y and comment on the result. For the latter, normalize w.r.t.
oth ||y|| and ||V ||||a||. What do you observe?
(b) Obtain values of the interpolating polynomial at closely spaced points using
the Matlab function polyval and plot the interpolating polynomial against
the data. Compute the condition number of the matrix V (x) (you can use a
Matlab function for that) and use it to explain the observed discrepancy.
2
(c) Consider y⃗ = x⃗. What interpolation problem does V (x⃗)⃗a = y solve? Use the
same (⃗x) as in (a) and find the interpolating polynomial. Comment on the
esult in relation to the interpolation in (a),(b).
Submit all commands you used and plots of the data and of the interpolating
polynomial for (b) and (c).
3. (a) Show that the Gauss-Seidel method for solving a linear system Ax = b is a
fixed point splitting method x = g(x), g(x = M−1(b+Nx) with split matrices
M = D + AL and N = −AU , where D is a diagonal matrix whose entries are
the diagonal entries of A, and AL, AU are lower and upper triangular matrices,
espectively, whose entries are the below/above the diagonal entries of A and
have zeros on the diagonal.
(b) Implement the Gauss-Seidel method in Matlab in either block or indicial form.
Use the stopping criterion ∥Axk − b∥ < tol. Your algorithm should print out
the number of iterations it took to achieve the tolerance and the final estimate
of x. It should check for e
or cases, such as when a division by zero occurs, and
handle them appropriately, e.g., you could use the Matlab command e
or.
(c) Use your Gauss-Seidel function to solve the sustem Ax = b, where
A =
XXXXXXXXXX
−1 −3 1 1 0
2 1 5 −1 −1
−1 −1 3 4 0
0 2 −1 1 4
b =
6
6
6
6
6
.
with a tolerance of 10−6. Assuming that the “exact” solution is given by the
Matlab command A\b, find the norm of the e
or in the solution.
(d) Count floating point operations (flops) needed to solve a linear system Ax = b,
where A is a fully populated matrix of dimension n by n, by the method of
Gauss-Seidel.
(e) Perform three LU decomposition steps of the above matrix by hand. Download
the functions lutx, forward, and backsubs from the NCM li
ary
https:
it.mathworks.com/mole
chapters.html,
where you can download Matlab functions and chapters of Cleve Moler’s book
“Numerical Computing with Matlab” (functions forward, and backsubs are
to be found in the text of bslashtx. Paste the two functions separately in
your editor window and create files, which you then can call). Obtain the L,U
components and the permutation vector p from lutx and then obtain the A\
esult by calls to forward and backsubs.
3
https:
it.mathworks.com/mole
chapters.html