數值分析實驗(一)之插值方法

實驗目的

通過設計、編制、調試2~3個多項式插值、拟合曲線的程式，加深對其數值計算方法及有關的基礎理論知識的了解。

實作拉格朗日（Lagrange）插值多項式、牛頓（Newton）插值、用線性函數拟合給定資料的程式。

問題

已知插值節點序列，用拉格朗日（Lagrange）插值多項式、牛頓（Newton）插值、用線性函數計算的函數在點的近似值。

我們選取函數y=3x+5

樣本點X=0,1,2,3,4,5,6,7,8,9；Y=5,8,11,14,17,20,23,26,29,32

拉格朗日（Lagrange）插值多項式

lagrange.m 檔案

用于編寫lagrange插值函數

%lagrange.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為m個插值構成的向量，R是誤差限 
function [y, R] = lagrange(X, Y, x, M)
n = length(X);
m = length(x); 
for i = 1:m 
    z = x(i); 
    s = 0.0; 
    for k = 1:n 
        p = 1.0; 
        q1 = 1.0; 
        c1 = 1.0; 
        for j = 1:n 
            if j~=k 
                p = p * (z - X(j)) / (X(k) - X(j)); 
            end
            q1 = abs(q1 * (z - X(j)));
            c1 = c1 * j; 
        end
        s = p * Y(k) + s; 
    end
    y(i) = s; 
    R(i) = M * q1 / c1; 
end
           

sy1_1.m檔案

用于運作求解得出結果

X = [0 1 2 3 4 5 6 7 8 9];
Y = [5 8 11 14 17 20 23 26 29 32]; 
x0 = 4.5; 
x = linspace(0,10,50);
M = 1; 
[y0, R] = lagrange(X, Y, x0, M); 
y = 3*x+5; 
errorbar(x0,y0,R,'.g')
hold on 
plot(X, Y, 'or', x0, y0, '.k', x, y, '-b'); 
legend('誤差','樣本點','拉格朗日插值估算','y=3x+5');
           

得出結果如圖：

牛頓（Newton）插值

newton.m 檔案

用于編寫牛頓插值函數

function yi=newton(x,y,xi)
%Newton基本插值公式
%x為向量，全部的插值節點
%y為向量，內插補點節點處的函數值
%xi為标量，是自變量
%yi為xi出的函數估計值
n=length(x);
m=length(y);
%計算均差表Y
Y=zeros(n);
Y(:,1)=y';
for k=1:n-1
    for i=1:n-k
        Y(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i));
    end
end
%計算牛頓插值公式
yi=0;
for i=1:n
    z=1;
    for k=1:i-1
        z=z*(xi-x(k));
    end
    yi=yi+Y(1,i)*z;
end
           

sy1_2.m檔案

用于運作求解得出結果

X = [0 1 2 3 4 5 6 7 8 9];
Y = [5 8 11 14 17 20 23 26 29 32]; 
x0 = 4.5; 
x = linspace(0,10,50);
y0 = newton(X, Y, x0); 
y = 3*x+5; 
R = y0-3*x0-5;
errorbar(x0,y0,R,'.g')
hold on 
plot(X, Y, 'or', x0, y0, '.k', x, y, '-b'); 
legend('誤差','樣本點','牛頓插值估算','y=3x+5');

           

得出結果如圖：

線性函數拟合

sy1_3.m檔案

用于運作求解得出結果

X = [0 1 2 3 4 5 6 7 8 9];
Y = [5 8 11 14 17 20 23 26 29 32]; 
x0 = 4.5;
x = linspace(0,10,50);
y = 3*x+5;
y1 = polyfit(X,Y,1);
y0 = x0*y1(1)+y1(2);
R = y0-3*x0-5;
errorbar(x0,y0,R,'.g')
hold on 
plot(X, Y, 'or', x0, y0, '.k', x, y, '-b'); 
legend('誤差','樣本點','線性函數拟合','y=3x+5');

           

得出結果如圖：