多项式插值（一）

参考《数值分析与科学计算》一书。

 

matlab里有大量关于插值的命令。

一、介绍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

二、利用已知范德蒙函数和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

 

三、用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时，插值开始振荡，效果不好。

 

四、使用内部命令进行插值

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

可见二者系数是相同的。

 

四、用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')

可见，用切比雪夫节点，除了中间部分近似较差之外，其余位置近似都不错。

 

五、多项式求值精度

>> 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的帮助中指出，当插值的结果不好时能够清除那些重复或接近重复的节点使数据集中（即移动数据，使横坐标或纵坐标的均值为零）以及从新调整数据的方法进行改进。