MATLAB例項：Hermite插值多項式

作者：凱魯嘎吉 - 部落格園http://www.cnblogs.com/kailugaji/

問題描述：清華大學出版社李慶揚《數值分析》第五版教材P49習題14：

求次數小於等於3的多項式P(x)，使其滿足條件P(0)=0, P’(0)=1, P(1)=1, P’(1)=2

1. MATLAB程式

Hermite_kailugaji.m

function f = Hermite_kailugaji(t,y,y_1,t0)
%t y為座標向量 y_1為一階導的值 t0為插值點的t座標|| f0為t0對應的值
syms x;
f = 0.0;

if(length(t) == length(y))
    if(length(y) == length(y_1))
        n = length(t);
    else
        disp('y和y的導數的維數不相等！');
        return;
    end
else
    disp('x和y的維數不相等！');
    return;
end

for i=1:n
    h = 1.0;
    a = 0.0;
    for j=1:n
        if( j ~= i)
            h = h*(x-t(j))^2/((t(i)-t(j))^2);
            a = a + 1/(t(i)-t(j));
        end
    end
    
    f = f + h*((t(i)-x)*(2*a*y(i)-y_1(i))+y(i));

    
    if(i==n)
        if(nargin == 4)
            f = subs(f,'x',t0);
        else
            f = vpa(f,6);
        end
    end
end
 

demo.m

clear
clc
% 求次數小於等於3的多項式P(x)，使其滿足條件
% P(0)=0, P’(0)=1, P(1)=1, P’(1)=2
x=[0 1];
y=[0 1];
y_1=[1 2];
%x, y為座標向量, y_1為一階導的值, x0為插值點的x座標|| f0為x0對應的值
f=Hermite_kailugaji(x,y,y_1);
f=collect(f); % Matlab 合併同類項，將多項式化為一般式
fprintf('Hermite插值多項式為: y=%s\n', f);
f=Hermite_kailugaji(x,y,y_1,2);
fprintf('在x=2處的Hermite插值為: %f\n', f);
 

2. 結果

Hermite插值多項式為: y=x^3 - 1.0*x^2 + 1.0*x
在x=2處的Hermite插值為: 6.000000

提示：MATLAB多項式整理小技巧：collect()是化為一般式，factor()是因式分解，合併同類項。