數值計算方法實驗之Newton 多項式插值(MATLAB代碼)

一、實驗目的

  在己知f(x),x∈[a,b]的表達式，但函數值不便計算或不知f(x),x∈[a,b]而又需要給出其在[a,b]上的值時,按插值原則f(x_i)=y_i (i=0,1,……, n)求出簡單函數P(x)(常是多項式)，使其在插值基點x_i處成立(x_i)= y_i(i=0,1,……,n),而在[a,b]上的其它點處成立f(x)≈P(x).

 二、實驗原理

   

三、實驗內容

    求f(x)=x⁴在[0,2]上按5個等距節點確定的Lagrange插值多項式

四、實驗程序

   (1).m文件

 

%輸入的量:X是n+1個節點(x_i,y_i)(i = 1,2, ... , n+1)橫坐標，Y是縱坐標，
%x是以向量形式輸入的m個插值點，M在[a,b]上滿足｜f~(n+1)(x)｜≤M
%注：f~(n+1)(x)表示f(x)的n+1階導數
%輸出的量：向量y是向量x處的插值，誤差限R，n次牛頓插值多項式L及其系數向量C，
%差商的矩陣A
function[y,R,A,C,L] = newton(X,Y,x,M)
n = length(X);
m = length(x);
for t = 1 : m
    z = x(t);
    A = zeros(n,n);
    A(:,1) = Y';
    s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;
        for j = 2 : n
            for i = j : n
                A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
            end
            q1 = abs(q1*(z-X(j-1)));
            c1 = c1 * j;
        end
        C = A(n, n); q1 = abs(q1*(z-X(n)));
        for k = (n-1):-1:1
            C = conv(C, poly(X(k)));
            d = length(C);
            C(d) = C(d) + A(k,k);%在最后一維，也就是常數項加上新的差商
        end
        y(t) = polyval(C,z);
        R(t) = M * q1 / c1;
end
L = poly2sym(C);

 

　（2）命令窗口輸入

X = [0 0.5 1.0 1.5 2.0];  
Y = [0 0.0625 1 5.0625 16];  
x = linspace(0,pi,50);  
M = 1;  
[y,R,A,C,L] = newton(X, Y, x, M);  
y1 = x.*x.*x.*x;  %可根據所給函數更改
errorbar(x,y,R,'.g')  
hold on  
plot(X, Y, 'or', x, y, '.k', x, y1, '-b');  
legend('誤差','樣本點','牛頓插值估算','x^4');

五、運算結果

    (1) 圖像

      

    (2) 運算結果

        第一列為所得多項式系數：