I am trying to write some code that does polynomial interpolation but I can't get it to work. I am a student and I am following the logic from this video https://www.youtube.com/watch?v=VpI-wC94RKw in order to recreate it it in code-form in Matlab, but whenever I try to create my own version of the matrix shown in the video I instead get a matrix almost exclusively filled with zeros (except one element). I can't understand why this is happening.
My code:
x=[150, 200, 300, 500, 1000, 2000, 99999]';
y=[2, 3, 4, 5, 6, 7, 8]';
function f = interPoly(x,y)
% Skapar en matris A var varje rad är [x_1^6, x_1^5,..., 1] till [x_n^6, x_n^5,..., 1]
A = [x.^6 x.^5 x.^4 x.^3 x.^2 x ones(numel(x),1) y];
% Gaussar matris A
R = rref(A);
% Plockar sista kolumnen ur R
c = R(:,end);
f = c(1)*x.^6+c(2)*x.^5+c(3)*x.^4+c(4)*x.^3+c(5)*x.^2+c(6)*x+c(7);
end
(The matrix 'A' is the one that's being problematic here. The function I get in the end is also just filled with zeros as values. Also sorry for the comments being in swedish)
I have 7 values in x and y and therefore a polynomial of the 6th order but I don't really know what the constant should be in the second last column so I just put a bunch of ones there (I am new to this so I am a little bit unsure about the logic).
Anyway I have tried the using the same function with some other input data and it has worked fine.
Alt input data:
x=[0, 0.5, 1, 1.5, 2, 2.99, 3]';
y=[0, 0.52, 1.09, 1.75, 2.45, 3.5, 4]';
Does it give me zeros because the elements overflow (for example 99999^6 is a very high number)? I don't really understand what is going on here and why it's working just fine with a different set of input data. Help?
Thanks!
EDIT: The entire point of this task (given by my school) is to compare the "least squares" method (which I have also written code for but not posted) with a polynomial interpolation method (the one in the code above). The last value in 'x' above is supposed to be infinity (f(inf)=8) so I just replaced it with a really high number, hence why it isn't "evenly" distributed. Is there a better way to do this?
In your example you get following matrix for R:
1.00000 0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 1.6925e-05
0.00000 1.00000 0.00051 0.00000 0.00000 0.00000 0.00000 0.00000 7.1286e-05
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 5.4078e-04
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 6.9406e-03
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 2.2098e-01
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 7.0000e+00
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 8.0000e+00
If the system of equation you define has one unique solution, then you should have an identity matrix (with this additional column on the right). In your case you have (appart from the first two lines) the equations
0 == 5.4078e-04
0 == 6.9406e-03
0 == 2.2098e-01
0 == 7.0000e+00
0 == 8.0000e+00
which indicate that your system does not admit a solution. This is due to your choice of inputs that in theory might have a solution, but in practice the numerical properties are very bad. And indeed, if we compute the condition number cond(A) we get a value of 4.3691e+06 which is very bad for such a small matrix. This is most likely due to the uneven spread of your x values. (You're trying to solve a Vandermonde system.) If you them in a more "evenly" distributed way things look a lot better.
Apart from that your code looks fine, but you probably want to evaluate the interpolation function in values other than the ones used to define them. Furthermore this "exact" approach is known to be not very ideal numerically, so it might be worth using a more stable algorithm for solving the system than rref, for example a qr decomposition. I suggest following changes:
function f = interPoly(x,y,xnew)
A = [x.^6 x.^5 x.^4 x.^3 x.^2 x ones(numel(x),1)];
[q,r] = qr(A);
c = r\(q' * y);
disp(cond(A));
f = c(1)*xnew.^6+c(2)*xnew.^5+c(3)*xnew.^4+c(4)*xnew.^3+c(5)*xnew.^2+c(6)*xnew+c(7);
end
x=[150, 200, 300, 500, 1000, 2000, 99999]';
y=[2, 3, 4, 5, 6, 7, 8]';
disp(interPoly(x,y,x+eps));