max3function val = max3(a, b, c)
if (a >= b) & (a >= c)
val = a;
elseif (b >= a) & (b >= c)
val = b;
else
val = c;
end
end
costaylorfunction val = costaylor(x, tol, kmax)
% costaylor Computes an approximation to cos(x) using Taylor series
%
% Input arguments
%
% x: The value of x at which the approximation is to be computed.
% tol: The tolerance for the computation. The function will
% return the approximation when the magnitude of the current
% term being summed is <= this value.
% kmax: The maximum number of terms that will be summed.
%
% Output argument
%
% val: The approximate value of cos(x).
%
% Also illustrates the use of the return statement to exit a function
% immediately (i.e. before the end of the function definition).
% Initialize the output argument.
val = 0;
for k = 0 : kmax - 1
% Compute k-th term in series.
term = (-1)^(k) * x^(2*k) / factorial(2*k);
% Update the approximation.
val = val + term;
% If the magnitude of the current term is smaller than
% the tolerance, exit the procedure using the return statement.
if abs(term) <= tol;
return;
end
end
% If execution reaches this point then convergence to the requested
% tolerance has not been achieved using the specified maximum
% number of terms. Set the output argument to NaN and return.
val = NaN;
end
We will be concerned with a version of the barycentric form of polynomial interpolation.
Recall from class notes that given the \( n \) data points
\[ (x_j, f_j), \quad j = 1, 2, \ldots, n \]
we first construct the barycentric weights, \( w_j \) via
\[ w_j = \frac{1}{\prod_{k=1, k\ne j}^n \left( x_j - x_k \right)} \]
Then defining
\[ l(x) = (x - x_1)(x - x_2) \ldots (x - x_n) \]
the interpolating polynomial \( p(x) \) is
\[ p(x) = l(x) \sum_{j=1}^n \frac{w_j}{x - x_j} f_j \quad\quad\quad\quad (1) \]
We now get an alternate form of the barycentric interpolating polynomial by letting the data values \( f_j \) be the constant function 1. Then from equation (1) above, we have
\[ 1 = l(x) \sum_{j=1}^n\frac{w_j}{x - x_j} \]
Dividing equation (1) by this last result and cancelling \( l(x) \) we obtain our final result for the barycentric interpolating polynomial, \( p(x) \):
\[ p(x) = \frac{\sum_{j=1}^n\frac{w_j}{x-x_j}f_j} {\sum_{j=1}^n\frac{w_j}{x-x_j}} \quad\quad\quad\quad (2) \]
Create a MATLAB function with the following header
function fto = barycentric(xfr, ffr, xto)
that implements barycentric interpolation according to equation (2). The input and output arguments to and from the function are defined as follows:
% Input arguments
%
% xfr: x values to fit (length n vector)
% ffr: y values to fit (length n vector)
% xto: x values at which to evaluate fit (length nto vector)
%
% Output argument
%
% yto: Interpolated values (length nto vector)
Although it is possible to write this function without the use of explicit loops, I recommend that you do use loops unless you are a MATLAB expert.
As usual, the source code for your function should be created in the folder
~/Documents/MATLAB
with a filename barycentric.m
For initial testing I suggest that you choose some low order polynomial, say a quadratic, and then generate vectors of test input data, \( (x_j, f_j) \), of appropriate length (3 for the case of a quadratic), and a short vector of output \( x \) values to feed to your routine.
Once you are reasonably satisfied that your implementation is working, test it with the following code:
function p = tpoly(x)
p = 1 - x + x.^2 - x.^3 + x.^4 - x.^5 + x.^6 - x.^7;
end
(define in file tpoly.m).
clf;
x = linspace(-0.9,0.9,101);
fxct = tpoly(x);
plot(x,fxct);
hold on;
xfr = linspace(-1,1,8);
ffr = tpoly(xfr);
xto = linspace(-0.9,0.9,101);
fto = barycentric(xfr,ffr,xto);
df = fxct - fto
plot(x,1.0e15 .* df);
title('Barycentric polynomial interpolation');
xlabel('x');
legend('Exact values', '10^{15} * Error in interpolated values');
(define in file tbarycentric.m).
Execute tbarycentric and analyze what you see in the plot
that is created. If your implementation is sound, the two curves appearing
in the plot should have approximately the same vertical range.