2

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?

Schytheron
  • 715
  • 8
  • 28
  • When I run the code my matrix `A` has nonzero entries. Putting ones is the right thing to do for the constant, because whatever the corresponding entry of c is will tell you the coefficient of 1. Also, since you are trying to solve the least squares problem `Ac = y` you could just use `A\y` where `A` doesn't have the `y` appended to it. This might give better results in generalA\y – overfull hbox Jan 23 '19 at 15:44
  • The column of your matrix are non linearly independant. So you get an "almost" singular matrix, so of course matlab cannot solve that efficiently. – obchardon Jan 23 '19 at 15:57
  • if your replace `A = [x.^6 x.^5 x.^4 x.^3 x.^2 x ones(numel(x),1) y]` by `A = [x.^6 x.^5 x.^4 x.^3 x.^2 x ones(numel(x),1)];` and compute c with `c = A\y` matlab should give you a warning. – obchardon Jan 23 '19 at 15:59
  • @tch I have also written code for a "least squares" method (but not posted it). I am now meant to write code strictly for polynomial interpolation. Read my edit in my post. It explains everything a bit better :) – Schytheron Jan 23 '19 at 16:17

1 Answers1

4

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));

Try it online!

flawr
  • 10,814
  • 3
  • 41
  • 71
  • Thanks for the explanation and thanks for the alternative solution. I have one small question regarding your code though: What does `[q,r] = qr(A);` mean in Matlab? You aren't assigning the result to one variable but instead to two variables? Also, you should read the edit I just posted... forgot to add some vital info. After reading the edit... do you think there exists an even better solution (how do I represent infinity?) ? – Schytheron Jan 23 '19 at 16:13
  • You're trying to solve linear system of equations (here written in matrix form) `A*c == y`. There are various ways to do this. The one you seem to be familiar with is called gaussian elimination, this is also what is used in `rref` and it is equivalent to solving this equation using the so called [LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition). This algorithm is numerically not very good for certain problems, which is why I used the so called [QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition) which is sligthyl slower but has some better numerical properties. – flawr Jan 23 '19 at 18:31
  • (If we were using exact arithmetic it would not matter which algorithm we would use.) – flawr Jan 23 '19 at 18:32
  • Now to the problem of infinity: In matlab there is the constant `Inf` in case you'll ever need it. But here we have a different problem: *Every* polynomial function `f(x) = an * x^n + .... + a1 * x + a0` that is not constant will diverge as x goes to infinity. There is no polynomial function that converges to a finite value as x approaches infinity. – flawr Jan 23 '19 at 18:35
  • Regarding the output of `qr`: Yes in MATLAB it is possible that a function returns *multiple* values. (Search for `output arguments` to find information on that.) – flawr Jan 23 '19 at 18:40
  • Thanks for the answers! :) – Schytheron Jan 23 '19 at 21:53