MatLab: Lagrange interpolation function graph deviates when nodes are greater than 28

Question

I don't know if this is the right place to ask this question. But I have the following code that interpolates a function $y$ base on the node $x$, which based on the termNumber. The function is $exp(3x)*cos(3x)$ and the nodes are base on $x=(2j/n)-1$ where $n=$termNumber and $j=0,...,n$

termNumber = 28;
x = (2.*(0:termNumber)./termNumber)-1;
y = exp(3.*x).*cos(3.*x);

sum=0;
for i=1:length(x)
    p=1;
    for j=1:length(x)
        if j~=i
            c = poly(x(j))/(x(i)-x(j));
            p = conv(p,c);
        end
    end
    term = p*y(i);
    sum= sum + term;
end
disp(sum);
K=flip(sum);

for  g=length(x):-1:2 %display the coefficient on polynomial form
    fprintf('+(%.2fx^%d)',K(g),g-1);
end

fprintf('+(%.2f)',K(1));

v=linspace(-1,1,100);
b = polyval(sum,v);
b1 = exp(3.*v).*cos(3.*v);
plot(v,b,'o') %display lagrange polynomial
hold on
plot(v,b1) %display actual function

The graph is displayed correctly using the above method when the termNumber is less than 28: see below:graph when termNumber = 28 where the red curve is the actual function(blue is the Lagrange interpolating polynomial). But when termNumber >= 29, the two curves started to separate: termNumber = 29 . And when termNumber gets large, both curves rea not recognizable anymore termNumber = 40 . Since I coded the method according to my textbook so I can't seem to figure out what I did wrong. Any inputs are appreciated.

Answer 1

There is no theoretical reason for the error, as that should fall like $(6\sqrt2/n)^{n+1}/(4(n+1))$. The observed error is thus due to cancellation errors and other floating point problems in the construction of the interpolation polynomial.

With a different, more compact computation, the deviation only starts to become visible at $32-35$, the error at $40+1$ nodes is $1/10$ of the error reported in the question.

  q = poly(x);
  dq = polyder(q)
  sum=0;
  for i=1:length(x)
      sum = sum + y(i)*deconv(q,poly(x(i)))/polyval(dq,x(i))
  end

One can get some nodes better by computing the interpolation polynomial of the error and add that to the polynomial, and then iterate that. This helps up to $41$. One might perhaps get to higher node densities using the Newton interpolating scheme, starting outside on both interval ends and moving alternatingly inwards. Or use some dedicated divide-and-conquer approach.