参考《数值分析与科学计算》一书。
matlab里有大量关于插值的命令。
1、介绍vander()和fliplr()两个与范德蒙有关的函数
>> x =[0 pi/2 pi 3*pi/2]；v =vander(x) v = 0 0 0 1.0000 3.8758 2.4674 1.5708 1.0000 31.0063 9.8696 3.1416 1.0000 104.6462 22.2066 4.7124 1.0000
>> v1 = fliplr(v) v1 = 1.0000 0 0 0 1.0000 1.5708 2.4674 3.8758 1.0000 3.1416 9.8696 31.0063 1.0000 4.7124 22.2066 104.6462
2、利用已知范德蒙函数和y，求系数
>> y =[0 1 0 -1]' y = 0 1 0 -1
>> p = v\y p = 0 1.6977 -0.8106 0.0860
>> p = flipud(p)' p = 0.0860 -0.8106 1.6977 0
3、用polyval命令计算给定系数的多项式的值
>> z= rand; p(1) *z^3 + p(2) *z^2 + p(3) *z +p(4) ans = 0.8916 >> polyval(p, z) ans = 0.8916
可见两个结果是一样的。
利用图形来对比
>> z = linspace(0, 2*pi, 100); >> zy = polyval(p, z); >> plot(z, zy, 'g', z, sin(z), 'r'), grid
其中红色是实际曲线，绿色是近似曲线。
改变：增加节点数
>> z = linspace(0, 4*pi, 200); >> zy = polyval(p, z); >> plot(z, zy, 'g', z, sin(z), 'r'), grid
可见大于3pi/2时，插值开始振荡，效果很差。
4、使用内部命令进行插值
>> x, y' x = 0 1.5708 3.1416 4.7124 ans = 0 1 0 -1 >> pnew = polyfit(x, y', 3) %x, y'是插值系数； 3为多项式系数。 pnew = 0.0860 -0.8106 1.6977 -0.0000 >> p p = 0.0860 -0.8106 1.6977 0
可见两者系数是相同的。
4、用10个等距节点进行拉格朗日插值、切比雪夫节点插值比较
>> x =linspace(0, 2*pi); n =9; >> plot(x ,abs(sin(x)), 'y') %图像在π附近不精确 >> nodes1 = linspace(0, 2*pi, n+1); %选择10个等距节点 >> ynodes1= abs(sin(nodes1)); %相应的y值 >> p= polyfit(nodes1, ynodes1, n); >> y1 = polyval(p, x); %在更细的划分上计算p(x)的值 >> hold on; plot(x, y1, 'g')
图中绿色图像为拉格朗日插值所做图像（10个节点），结果并不理想
下面进行改善（通过切比雪夫节点）
>> nodes1 >> plot(nodes1, zeros([1 n+1]),'b*') >> cheby= cos((2* (n-(0:n)) + 1) * pi/ (2*n +2)) >> nodes2 = cheby * pi + pi; >> hold on; plot(nodes2, zeros([1 n+1]), 'mx')
由图可知，切比雪夫节点并不是等距的，而是越接近端点节点就分布越密。
接下来对三图进行比较
>> ynodes2 = abs(sin(nodes2)); >> p2 =polyfit(nodes2, ynodes2, n); >> y2 = polyval(p2, x); >> plot(x ,abs(sin(x)), 'y') >> hold on; plot(x, y1, 'g') >> hold on; plot(x, y2, 'r')
可见，用切比雪夫节点，除了中间部分近似较差以外，其他位置近似都不错。
5、多项式求值精度
>> N=20; y= 10^N, y1=(10 * (1+2*eps)) ^N >> abs(y1 -y) %绝对误差 ans = 704512 >> ans/y %相对误差 ans = 7.0451e-15
绝对误差很小，而相对误差非常大。扰动和指数越大，结果越差。
这里简单介绍eps这个常量：传送门
简单来说：
eps是matlab中最小的正数。eps=2.22044604925031e-016 在matlab的数值计算中，当发现某个值小于eps时，就把这个数当做0来处理。 这也可以看做是matlab的精度值。
同时，抵消也是潜在的问题。
>> x = linspace(.995, 1.005, 200); >> plot(x, (x-1).^8) %f(x)
>> g = x.^8 - 8*x.^7 + 28*x.^6 - 56*x.^5 + 70*x.^4 -56*x.^3 + 28*x.^2 -8*x +1; >> hold on; plot(x, g), grid %g(x)
f(x)和g(x)是一样的，g为f的展开式，然而所得到的的图像却是不一样的。
但由减法产生的有效数字的抵消使得在x=1附近产生了很大的干扰函数。
当然抵消问题是可以得到解决的。在拉格朗日和牛顿形式中，如果要计算多项式p(x)在x=x0附近的值，一般写成(X-Xn)的形式，求(X-Xn)的幂而不是X的幂。
另一方面，可以用递归求值(霍纳方法)精确计算多项式的值。
>> a = 3*z^3 - 6*z^2 + z^2 + 4*z -5 a = 51.2372 >> b = ((3*z-6)*z+4)*z-5 b = 41.3676
其实a和b是一样的。
使用递归的形式可以使浮点数运算次数少一些。
？：polyfit的帮助中指出，当插值的结果很差时可以清除那些重复或接近重复的节点使数据集中（即移动数据，使横坐标或纵坐标的均值为零）以及重新调整数据的方法进行改进。
