lab5

lab5function.m

function [ L,A,B,C,D ] = mydivdiff( x, f )
% INPUTS: x - Nx1 vector of x values at which data for the function is specified
% f - Nx1 vector of corresponding function values
% OUTPUTS: A - (N-1)x(N-1) table of divided differences

n = size(x,1);
L = zeros(n-1);

L(:,1) = (f(2:n)-f(1:n-1))./(x(2:n)-x(1:n-1));

for i = 2:n-1
    L(1:n-i,i) = (L(2:n-i+1,i-1)-L(1:n-i,i-1))./(x(2:n-i+1)-x(1:n-i));
end

n = size(x,1);
A = zeros(n,n+1);
A(:,1) = f;
for i = n-1:-1:1
    for j = 2:n-i+1
        A(i,j) = (A(i+1,j-1) - A(i,j-1)) / (x(i+j-1) - x(i));
    end
end
A = A(:,2:n);

B = zeros(n,n+1);
B(:,1) = f;
for j = 2:n
    for i = n-j+1:-1:1
        B(i,j) = (B(i+1,j-1) - B(i,j-1)) / (x(i+j-1) - x(i));
    end
end
B = B(:,2:n);

X = x;
Y = f'; 
n = length(X);
D = zeros(n,n);
D(:,1) = Y';
for j=2:n,
  for k=j:n,
      D(k,j) = (D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1));
  end
end
C = D(n,n);
for k=(n-1):-1:1,
  C = conv(C,poly(X(k)));
  m = length(C);
  C(m) = C(m) + D(k,k);
end

end

lab5script.m

clear all
close all

x = [0;5;10;20;30;40];
y = [1.787;1.519;1.307;1.002;0.7975;0.6529];

figure()
hold on
plot(x,y,'o')
xi = 7.5;
yi = mypolyint(x(2:3),y(2:3),1,7.5);
fprintf('interpolated value Q1a: %7.7f\n\n',yi);
plot(x(2):0.01:x(3),mypolyint(x(2:3),y(2:3),1,x(2):0.01:x(3)),'r-');
plot(xi,yi,'ro')
xlim([-5 45]);
ylim([0 2]);

figure()
hold on
plot(x,y,'o')
xi = 7.5;
yi = mypolyint(x(1:3),y(1:3),2,7.5);
fprintf('interpolated value Q1b: %7.7f\n\n',yi);
plot(x(1):0.001:x(3),mypolyint(x(1:3),y(1:3),2,x(1):0.001:x(3)),'r-');
plot(xi,yi,'ro')
xlim([-5 45]);
ylim([0 2]);

x = [88.9;94;99.1;104.0;108.5;116.8;127;139.7];
y = [14.6;16.1;16.6;15.3;16.7;18.1;19.5;23.2];
figure()
hold on
plot(x,y,'o')
X = 120;
P = polyfit(x,y,1);
Y = polyval(P,X);
fprintf('interpolated value Q2: %7.7f\n\n',Y);
X = min(x):0.1:max(x);
Y = polyval(P,X);
plot(X,Y,'r-')
xlim([80 145]);
ylim([12 24]);

A = mydivdiff(x,y);


a = [2;3;5];
b = [-3;-2;6];