数值计算方法——插值公式的实现

lagrange插值公式

function Ln=lagrange(x,y,w)
% x=input('请输入x向量：');
% y=input('请输入y向量：');
% w=input('请输入待求值：');
Ln=0;
for i=1:n
    z=1;%基函数，每次循环前置为1,这个很重要
    for j=1:n
        if i~=j
            z=z.*(w-x(j))/(x(i)-x(j));
        end   
    end
    Ln=Ln+z.*y(i);
end

%调用函数脚本文件，运行可输出。
clear;
clc;
fprintf('输入一组节点（x,y)以及待求x0。其中y为函数值')
while 1
    x=input('x');
    y=input('y');
    w=input('x0');
    m=length(x);
    n=length(y);
    if m==n
        break
    end
end
%参考数据
% x=[100 121 144];
% y=[10 11 12];
%w=125
Ln=lagrange(x,y,w);
disp(Ln)

两点Hermite插值公式

%函数名：phi0.m
function phi0=phi0(t)
phi0=(1-t)*(1-t)*(1+2*t);

%函数名：phi1.m
function phi1=phi1(t)
phi1=t*(1-t)^2;

%函数名：Hermite.m
function H3=Hermite(x,y0,y1,w)
m=length(x);
n=length(y0);
r=length(y1);
if m==n&&n==r
%     phi0(t)=(1-t)^2*(1+2t);
%     phi1(t)=t*(1-t)^2;
    h=x(2)-x(1);
    t=(w-x(1))/h;
    h00(w)=phi0(t);
    t=(x(2)-w)/h;
    h10(w)=phi0(t);
    t=(w-x(1))/h;
    h01(w)=h*phi1(t);
    t=(x(2)-w)/h;
    h11(w)=-h*phi1(t);
    H3=y0(1)*h00(w)+y0(2)*h10(w)+y1(1)*h01(w)+y1(2)*h11(w);
end

%脚本文件运行可输出，调用了前三个函数
clear;
clc;
fprintf('输入一组节点（x,y0,y1)以及待求x0。其中y0为函数值，y1为一阶导数值')
while 1 %提高健壮性，避免输入错误的数
    x=input('x');
    y0=input('y0');
    y1=input('y1');
    w=input('x0');
    m=length(x);
    n=length(y0);
    r=length(y1);
    if m==2&&m==n&&n==r
        break
    end
end 
%参考数据
% x=[122 144];
% y0=[11 12];
% y1=[1/22 1/24];
% w=125;
H3=Hermite(x,y0,y1,w);
disp(H3)

分段低次插值公式

function s=dici(x,y,x0)
t1=x0-x(2);
t2=x0-x(1);
b1=x(1)-x(2);
b2=x(2)-x(1);
s=(t1/b1)*y(1)+(t2/b2)*y(2);

function f=FenD(x,y,x0)
%输入一组节点（x,y)
%输入待求值x0
n=length(x);
for a=0:n-1
    if x0<x(a+1)&&x0>x(a)
        x1=[x(a),x(a+1)];
        y1=[y(a),y(a+1)];
        f=dici(x1,y1,x0);
    end
end

%测试代码
clear;
clc;
fprintf('输入一组节点（x,y)以及待求x0。其中y为函数值')
while 1
    x=input('x');
    y=input('y');
    x0=input('x0');
    m=length(x);
    n=length(y);
    if m==n
        break
    end
end
% x=[11 12 34 45];
% y=[23 34 66 45];
% x0=33;
y0=FenD(x,y,x0);
disp(y0)

分段高次插值公式

function s=gaoci(x,y0,y1,x0)
%函数名gaoci.m
h=zeros(2,2);
d=x(2)-x(1);
t1=(x0-x(1))/d;
h(1,1)=phi0(t1);
t2=(x(2)-x0)/d;
h(2,1)=phi0(t2);
t3=(x0-x(1))/d;
h(1,2)=d*phi1(t3);
t4=(x(2)-x0)/d;
h(2,2)=-d*phi1(t4);
s=y0(1)*h(1,1)+y0(2)*h(2,1)+y1(1)*h(1,2)+y1(2)*h(2,2);

phi0和phi1函数已经在两点Hermite插值公式代码处给出！

function f=FenG(x,y0,y1,x0)
%输入一组节点（x,y0,y1)其中y0为函数值，y1为一阶导数值
%输入待求值x0
n=length(x);
for a=1:n-1
    if x0<x(a+1)&&x0>x(a)
        b=[x(a),x(a+1)];
        b0=[y0(a),y0(a+1)];
        b1=[y1(a),y1(a+1)];
        f=gaoci(b,b0,b1,x0);
    end
end

clear;
clc;
%分段高次插值公式测试代码
fprintf('输入一组节点（x,y0,y1)以及待求x0。其中y0为函数值，y1为一阶导数值')
while 1 %提高健壮性，避免输入错误的数
    x=input('x');
    y0=input('y0');
    y1=input('y1');
    x0=input('x0');
    m=length(x);
    n=length(y0);
    r=length(y1);
    if m==n&&n==r
        break
    end
end 
% x=[12 14 17 87 122 144];
% y0=[ 1 2 4 3 11 12];
% y1=[2 3  1 1 1/22 1/24];
% x0=125;
y0=FenG(x,y0,y1,x0);
disp(y0)