Write a routine in MATLAB that does the following: given an interval [a, b], a f

Question

Write a routine in MATLAB that does the following: given an interval [a, b], a function F, and a number n, your code should produce the following three interpolants to the function F: (i) a polynomial interpolant to F at n equidistant points on [a, b]; (ii) a polynomial interpolant at n Chebyshev points of that interval; (iii) a cubic spline interpolant at n equidistant points. You are allowed to look at our class demo subdir, but you must document this fact (and you may not use the matlab library, such as polyfit). Furthermore, the polynomial interpolation must be done via a table of divided differences.

Explanation / Answer

1.Interpolation with equidistant nodes

We want to approximate the function f(x) = 1/(1+x^2) on the interval [-5,5]

f = @(x) 1./(1+x.^2);          % define function f

a = -5; b = 5;                 % endpoints of interval

n = 11;                        % number of nodes for interpolation

xt = linspace(-5,5,1000);      % use these x values for plotting

xe = a + (0:n-1)*(b-a)/(n-1); % equidistant nodes, same as xe = linspace(a,b,n)

ye = f(xe);                    % find values of f at nodes xe

d = divdiff(xe,ye);            % use divided difference algorithm to find coefficients d

pt = evnewt(d,xe,xt);          % evaluate interpolating polynomial at points xt

plot(xt,f(xt),xt,pt,xe,ye,'o'); grid on % plot function f and interpolating polynomial

                                         % mark given points with 'o'

2. Interpolation with Chebyshev nodes

xc = (a+b)/2 + (b-a)/2*cos(pi/n*((1:n)-.5));      % Chebyshev nodes

yc = f(xc);                    % find values of f at nodes xc

d = divdiff(xc,yc);            % use divided difference algorithm to find coefficients d

pt = evnewt(d,xc,xt);          % evaluate interpolating polynomial at points xt

plot(xt,f(xt),xt,pt,xc,yc,'o'); grid on % plot function f and interpolating polynomial

                                         % mark given points with 'o'

3. A cubic spline interpolant at n equidistant points

I'm setting up a function to output to evaluate the spline function in certain points based on a given set of abscis points with its function values, but it's giving me very deviated results.

The inputs of my function are a vector of abscis points (x_a) with its function values (f) and also a vector of equidistant interpolation points in the interval (t) from which I want the function values of the spline as the outputs. So I test my code with the following vectors: x_a is a 1x10 double vector of 10 equidistant points in the interval [-2,2], f are the sin(x)-function values of x_a and t is a 1x20 double vector of 20 equidistant points in the interval [-2,2]:

x_a = linspace(-2,2,10);

f = sin(x);

t = linspace(-2,2,20);

This is my function (based on the natural interpolating cubic spline method):

function [ y ] = naturalspline( x_a, f, t )

if length(x_a) ~= length(f)

error('vectors x_a and f must have same length');

end

n = length(x_a);

d = zeros(n,1); %vector of differences between abscis points

d(1)=x_a(1);

for i = 2:n

d(i) = x_a(i) - x_a(i-1);

end

% Coefficientmatrix C (tridiagonal):

C = zeros(n);

% Natural Spline boundary conditions:

C(1,1)= 2*(d(1)+d(2));

C(n,n) = 2*(d(n-1)+d(n));

for j = 2:n-1

C(j,j-1) = d(j-1);

C(j,j) = 2*(d(j-1)+d(j));

C(j-1,j) = d(j);

end

% Vector b:

b = zeros(n,1);

for i = 2:n-1

b(i) = (6/d(i))*(f(i+1)-f(i)) - (6/d(i-1))*(f(i)-f(i-1));    

end

% Vector s (second derivatives, solution of tridiagonal system):

s = C;

% Integration constant c1 :

c1 = zeros(n,1);

for i = 1:n

c1(i) = (f(i)/d(i)) - d(i)/(6*s(i));

end

% Integration constant c2:

c2 = zeros(n,1);

for i = 1:n-1

c2(i) = (f(i)/d(i+1)) - d(i+1)/(6*s(i));

end

%Setting up polynomials for each interval

for i = 1:n-1

p{i} = @(x) [s(i+1)*((x-x_a(i))^3)/(6*d(i+1)) - s(i)*((x-x_a(i+1))^3)/(6*d(i+1)) ...

      + c1(i+1)*(x-x_a(i)) + c2(i+1)*(x_a(i+1)-x)];

end

% evaluating spline function values of the interpolation points

j = length(t);

y = zeros(j,1);

for i = 1:j

[~,~,k] = histcounts(t(i),x_a); %k decides which interval t(i) is located

y(i) = p{k}(t(i));

end

Related Questions

Navigate

• Browse (All)
• Browse W
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.