Numerical Analysis Newton’s Divided Difference and cubic spline interpolation 1) Write a function divdiff in matlab. The function should read in two vectors x and y, and returns a table (a matrix) of the divided difference values. The prototype of the function is: function a = divdiff(x, y) finput: x, y: the data set to be interpolated foutput: a: a table for Newton’s divided differences 2) Write a function called polyvalue, that reads in the table of divided difference values generated by the function divdiff, the x-vector, and a vector t. The output of the function should be the values of Newton’t polynomial computed at points gives in the vector t. The prototype of the polyvalue function is function v = polyvalue (a, x, t) $ input: a = Newton’s divided differences x = the ponts for the dataset, same as in divdiff t = the points where the polynomial should be evaluated 음 음 3) Test your Matlab functions in a and b using 11 equally spaced nodes on the interval [-5, 5] >>x=-5:5), find the interpolating polynomial p of degree 10 for the function f(x) = (x2 +1)-1. Plot the function f(x) using the matlab ezplot function. Then on the same plot (use hold on) plot p(x) at the 11 point. Does it work well or poor? explain why? 4) Now, use the cubic spline interpolation provided by the matlab functions spline and ppval to generate the interpolation assuming natural cubic spline. Plot the result on the same plot as in three. 5) Repeat (4), but now use clamped cubic spline instead. You can do that by differentiating the function and finding f(–5) and f(.5). Plot the result on the same plot as in 3 and 4. Your final result should look something like the following figure. 1/(x2 + 1) 1 0.9 actual nodes interpolated natural cubic …. clamped 0.8 0.7 mmmmmmm 0.6 f(x) 0.5 0.4 0.3 0.2 0.1 -3 -2 – 1 0 1 3 N x What to turn in: Submit a zip file that contains the required matlab files that would generate a figure similar to the above figure. Explain which method is the best for this problem, natural cubic spline or clamped cubic spline or Newton’s divided difference.

42 1

Get full Expert solution in seconds

$1.97 ONLY

EXPERT ANSWER

%%Matlab code for Newton divide difference formulation
clear all
close all
%function for which interpolation have to do
f=@(x) (x.^2+1).^-1;
%all nodes
xx=linspace(-5,5,11);
yy=f(xx);

t=linspace(-5,5,100);

a = divdiff(xx,yy);
v = polyvalue(a,xx,t);

%function after interpolation
p=@(x1) 0.02036*x1 + 0.01041*(x1 + 4.0)*(x1 + 5.0) + 0.006335*(x1 + 4.0)*(x1 + 5.0)*(x1 + 3.0) + 0.004299*(x1 + 2.0)*(x1 + 4.0)*(x1 + 5.0)*(x1 + 3.0) – 0.002036*(x1 + 1.0)*(x1 + 2.0)*(x1 + 4.0)*(x1 + 5.0)*(x1 + 3.0) – 0.001131*x1*(x1 + 1.0)*(x1 + 2.0)*(x1 + 4.0)*(x1 + 5.0)*(x1 + 3.0) + 0.001086*x1*(x1 – 1.0)*(x1 + 1.0)*(x1 + 2.0)*(x1 + 4.0)*(x1 + 5.0)*(x1 + 3.0) – 0.0004299*x1*(x1 – 1.0)*(x1 + 1.0)*(x1 – 2.0)*(x1 + 2.0)*(x1 + 4.0)*(x1 + 5.0)*(x1 + 3.0) + 0.0001131*x1*(x1 – 1.0)*(x1 + 1.0)*(x1 – 2.0)*(x1 + 2.0)*(x1 + 4.0)*(x1 + 5.0)*(x1 – 3.0)*(x1 + 3.0) – 2.262e-5*x1*(x1 – 1.0)*(x1 + 1.0)*(x1 – 2.0)*(x1 + 2.0)*(x1 – 4.0)*(x1 + 4.0)*(x1 + 5.0)*(x1 – 3.0)*(x1 + 3.0) + 0.1403;
ezplot(f)
hold on
scatter(xx,yy,’g’,’filled’,’linewidth’,3)
ezplot(p)
pp = spline(xx,yy);
yy1 = ppval(pp, t);
plot(t,yy1)
yy2 = pchip(xx,yy,t);
plot(t,yy2)

xlim([-3 3])
ylim([0 1])
legend(‘Actual’,’nodes’,’interpolated’,’natural cubic’,’clamped cubic’)
xlabel(‘x’)
ylabel(‘f(x)’)
title(‘x vs. f(x) plot’)
box on

%Matlab function for Newton’s Divided Difference
function a = divdiff(x,y)
    n=length(x);
    y2=y;
    zz=zeros(n,n+1);

    zz(:,1)=x’; zz(:,2)=y’;

    for i=1:n-1

        n1=length(y2);
        for j=1:n1-1
            y3(j)=(y2(j+1)-y2(j))/(x(i+j)-x(j));
            zz(j,i+2)=y3(j);
        end
        a(i)=y3(1);
        y2=y3;
        clear y3;
    end
    nn=length(zz);
    fprintf(‘Newton divide difference table.\n’)
    for i=1:n
        fprintf(‘\t%d\t%4.2f\t%4.2f\t%4.2f\t%4.2f\t%2.2f\t%2.2f\n’,…
            zz(i,1),zz(i,2),zz(i,3),zz(i,4),zz(i,5),zz(i,6),zz(i,7))
    end
    a=[a y(1)];
end


function v = polyvalue(a,x,t)

    syms x1
    %loop for creating the function
    n=length(a);
        for i=1:n-1
            s1=1;
            for j=1:i
                    s1=(x1-x(j))*s1;
             end
             zz1(i)=(s1)*a(i);
        end
        f_newton(x1)=sum(zz1)+a(n);

    fprintf(‘\nThe Interpolation function using Newton Divided Difference Method of order %d\n’,n)
    disp(vpa(f_newton,4))
    v=double(f_newton(t));
  
end