Description
We investigate interpolation of the function
f : [−1, 1] 7→ R, f(x) = 1
1 + 30x
2
,
very similar to a famous example, the “Runge” function, used to exhibit – yes you guessed it – the
“Runge phenomenon”.
f = @(x) 1./(1 + 30*x.^2 );
M1 Polynomial interpolation with equidistant points, xk = −1+2k/N, k = 0, 1, 2, . . . N. Weights
for the barycentric formula are wk = (−1)k
N
k
!
.
M2 Polynomial interpolation with Chebyshev points of the second kind,
xk = cos(kπ/N), k = 0, 1, 2, . . . N.
Weights for the barycentric formula are wk = (−1)k
, except w0 =
1
2
, wN = (−1)N 1
2
.
M3 Mix and match: Use barycentric formula with equidistant points, but with Chebyshev
weights. In this case the interpolating function is no longer a polynomial, but a rational
function.
Your coding tasks.
a. You will need to write code to compute the weights for the barycentric formula. Matlab has a
built-in command nchoosek; You may use that, or compute binomial coefficients recursively,
N
0
= 1,
N
k+1
=
N
k
N−k
k+1 .
b. You will have to write code implementing the barycentric formula,
p(x) =
PN
m=0
wmfm
P
x−xm
N
m=0
wm
x−xm
.
This can be done to varying degree of efficiency and flexibility. Remember, that Matlab likes
vectors; for loops not so much. Ideally, your function can take an array of evaluation points
as inputs, and return the values at all of those points simultaneously.
One detail you will have to pay attention to is the case x − xk = 0, i.e., an evaluation point
coinciding with an interpolation point. In that case your function should simply assign the
given value fk to p(xk). If your barycentric code allows vectors as input arguments, then the
Matlab function find is useful for this task. Otherwise, a simple if statement will do.
mr
Macm 316 Computing Assignment #5
For your report compare the interpolated values to the exact values of the function f at 321 evenly
space points on the interval [−1, 1]. So finding the maximum error on the interval, translates into
finding the largest error at the points t(j) below:
t = linspace(-1,1,321); fex = f(t);
For all methods plot the maximum errors against N.
M1 p1 Behaviour on the whole interval [−1, 1]. Interpolate with polynomials of degree N =
8, 12, 16, 20; compute the maximum error on the interval [-1,1]. Does the error decrease
or increase with N? What is the rate of decay/growth?
p2 Behaviour on the smaller interval [−
1
2
,
1
2
]: Here we use exactly the same computation as
in part [p1], i.e., we are still interpolating on the interval [−1, 1]. We are, however, only
looking at the errors on the smaller interval [−
1
2
,
1
2
], i.e., only at t(j), with 81 ≤ j ≤ 241.
Confirm numerically that there is convergence on this interval, by interpolating with
some higher-degree polynomials, say N = 40, 80.
M2 Interpolate with polynomials of degree N = 8, 16, 20, 80, 200, 250; compute the maximum
error on the interval [-1,1]. Does the error decrease or increase with N? What is the observed
rate of decay/growth? What is the lowest error achieved? Do you gain accuracy by increasing
N from 200 to 250?
M3 Compute again with N = 8, 16, 20, 80, 200, 250; compute the maximum error on the interval [-
1,1]. Does the error decrease or increase with N? What is the observed rate of decay/growth.
What is the lowest error achieved?
Beyond this laundry list of things to do, I am not giving you any additional instructions as how to
present your report, just suggestions. You probably want to consider:
1. Including at least one plot that shows f and some of its interpolating function, but certainly
not more than four. One might be fine.
2. What kind of plot should you use to graph the maximum errors vs N: plot; semilogy;
or loglog? Generally, when graphing quantities of very different magnitudes, a logarithmic
plot is better; when graphing the function f and its approximations using plot might be
preferable. Some or perhaps all your error plots could be combined into one plot.
3. The key part of your coding in this assignment is the barycentric formula. Make sure to
feature it prominently, and include comments in your code.
mrt 2
