Newton polynomial interpolating points, matrix too big

Hi,

My code to make a newton polynomial and plot it does not work. I cannot figure out what the problem is. "Matrix dimensions must agree". It seems my matrix A grows too much.

f = @(x) x.^2.*sin(x);
x = [0 1 2 3 5 7 8];
y= f(x);

k = length(x)
ak=ones(k,1);
A = ak;

for column=2:k
    ak = ak.*(x-x(column-1));
    A =[A ak];
end
c=A\y
p =@(x) c(1) + c(2)*x + c(3)*x.^2+c(4)*x.^3;

tv = 0:0.1:6;
plot(tv, p(tv), 'r')
hold on;
plot(tv, f(tv), 'b')

Answer 1

[Alex Mcaulley]

Try changing this line

c=A\y

by

c=A\y'

[am]

Thank you Alex.

Do you know why my interpolation ended up completely flat?

[darova]

The problem is in dimensions

x = [0 1 2 3 5 7 8];        % row
ak=ones(k,1);               % column
ak = ak.*(x-x(column-1));   % multiplying column.*row    try: ak .* (x-x(column-1))'
c=A\y'
p =@(x) c(1) + c(2)*x + c(3)*x.^2+c(4)*x.^3;        % using only 4 constants of 7?

Little suggestion:

A = zeros(length(x));            % pre-allocating
A(:,column) = ak;
for column=2:k
    ak = ak.*(x-x(column-1))';
    A(:,column) = ak;
end

[am]

I tried to update the code according to your suggestions and I got a singular matrix unfortunately.

f = @(x) x.^2.*sin(x);
x = [0 1 2 3 5 7 8];
y= f(x);

k = length(x)
A = zeros(length(x));            
ak=ones(k,1);      
A(:,column) = ak;
for column=2:k
    ak = ak.*(x-x(column-1))';
    A(:,column) = ak;
end
A
c=A\y'
p =@(x) c(1) + c(2).*x + c(3).*x.^2+ c(4).*x.^3 + c(5).*x.^4+ c(6).*x.^5 + c(7).*x.^5;

tv = 0:0.1:6;
plot(tv, p(tv), 'r')
hold on;
plot(tv, f(tv), 'b')

[darova]

A(:,column) = ak;            % what 'column' equals to?
for column=2:k
    ak = ak.*(x-x(column-1))';
    A(:,column) = ak;
end 

There is a mistake near c(7) (must be c(7)*x.^6)

p = @(x) c(1) + c(2).*x + c(3).*x.^2+ c(4).*x.^3 + c(5).*x.^4+ c(6).*x.^5 + c(7).*x.^5;

[am]

Hi, it still does not give me a curve.

[darova]

You didn't answer my question: "What 'column' variable equals to" in line before loop?

A(:,column) = ak;

[am]

Sorry, this line is not my code so I don't know what it refers to?

My code:

k = length(x)
ak=ones(k,1);
A = ak;

[darova]

You change your code according to my suggestion but forgot to assign:

column = 1;
A(:,column) = ak;

When you run your code the value 'column' is probably: column = k = length(x) (from previous run)

Also, this part can be automated:

p = @(x) c(1) + c(2).*x + c(3).*x.^2+ c(4).*x.^3 + c(5).*x.^4+ c(6).*x.^5 + c(7).*x.^5;
% can be
p1 = tv*0;                % pre-allocating
for i = 1:length(p1)
    for j = 1:length(c)
        p1(i) = p1(i) + c(j)*x(i)^(j-1);
    end
end
plot(tv,p1,'r')

[am]

Thank you very much for this thorough answer.

[darova]

♥