多項式插值（一）

參考《數值分析與科學計算》一書。

 

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的幫助中指出，當插值的結果不好時能夠清除那些重複或接近重複的節點使數據集中（即移動數據，使橫座標或縱座標的均值爲零）以及從新調整數據的方法進行改進。