插值MATLAB实现（牛顿差商、插值误差、龙格现象、切比雪夫插值）

牛顿差商

function [c,y] = newtondd(a,b,x)
n=length(a);
for i=1:n
    v(i,1)=b(i);
end
for j=2:n
    for i=1:n-j+1
        v(i,j)=(v(i+1,j-1)-v(i,j-1))/(a(i+j-1)-a(i));
    end
end
for j=1:n
c(j)=v(1,j);
end
k=length(x);
y=zeros(1,k);
for w=1:k
for i=n:-1:2
    f=1;
    for j=1:i-1
        f=f*(x(w)-a(j));
    end
    f=f*c(i);
    y(w)=y(w)+f;
end
y(w)=y(w)+c(1);
end
end

说明：a为多项式插值基点横坐标，b为多项式插值基点纵坐标，x为输入待求值。y为利用牛顿差商公式得到的多项式的x对应值，c为牛顿差商公式系数。

插值误差

function z = chazhierror(f,h,x0)
n=length(h);
d=1;
for i=1:n
    d=d*i;
end
syms x
g=f(x);
fn=diff(g,n);
for i=1:n
    a(i)=subs(fn,x,h(i));
end
b=max(abs(a));
k=length(x0);
for j=1:k
    l=1;
    for i=1:n
        l=l*(x0(j)-h(i));
    end
    z(j)=b/d*l;
end
z=double(z);
end

f为模拟插值函数，h为模拟插值横坐标，z为在x0点处的误差上限。

切比雪夫插值

function [d,d1,cha] = qiebixuefu(f,n,a,b)
deta=(b-a)/n;
h=[a:deta:b];
h1=(a+b)/2+(b-a)/2*cos((1:2:2*n-1)*pi/(2*n));
syms x
g=f(x);
y=subs(g,x,h);
y1=subs(g,x,h1);
y=double(y);
y1=double(y1);
l=a:0.01:b;
[c,d] = newtondd(h,y,l);
[c1,d1] = newtondd(h1,y1,l);
d2=subs(g,x,l);
cha=d-d1;
plot(l,d,'-',l,d1,'o',l,d2,'*');
end

f为模拟插值函数，n为插值点数，a，b为基本域。d为等间距插值得到的多项式对应横坐标l的函数值，d1为切比雪夫插值得到的多项式对应横坐标l的函数值。

龙格现象图示

f =

  包含以下值的 function_handle:

    @(x)exp(abs(x))

>> a=-1

a =

    -1

>> b=1

b =

     1

>> n=10

n =

    10

图一为上述参数的等间距插值与切比雪夫插值对比，图二为局部放大。可以看到切比雪夫插值更接近原始函数，避免了龙格现象。