Lagrange、Newton和Hermite插值例题

污浊的双黑

已于 2023-04-12 17:37:13 修改

阅读量567

点赞数

文章标签： matlab Powered by 金山文档

于 2023-03-25 09:58:16 首次发布

本文链接：https://blog.csdn.net/qq_53011270/article/details/129763230

版权

问题1：

分别使用Lagrange和Newton插值公式计算函数的值。对于Newton插值要求先输出差商表，在计算x=7.5时4次插值多项式的值，并与四次Lagrange插值结果进行比较。

x	4	5	6	7	8
y	2.0000	2.23607	2.44949	2.64575	2.82843

Lagrange插值

代码1

function y = Lagrange(a,b,x,n)
%Lagrange插值法 a,b分别为x以及对应的f(x),x为插入数据，n为插值多项式次数
% a=4:8
% b=[2 2.23607 2.44949 2.64575 2.82843]
if length(a)~=length(b)
    error('输入数据不同维')
end
if n>=length(a)
    error('n值太大')
end
d=abs(x-a);
[~,idx]=sort(d);
S=0;
for i=1:n+1
    t1=1;
    t2=1;
    for j=1:n+1
        if j~=i
            t1=t1*(x-a(idx(j)));
            t2=t2*(a(idx(i))-a(idx(j)));
        end
    end
    S=S+t1/t2*b(idx(i));
end
y=S;
end

代码2

function y=lagrange_(xx,yy,x)
n=1:length(xx);
y=zeros(size(x));
for i=n
    ij=find(n~=i);
    y1=1;
    for j=1:length(ij),y1=y1.*(x-xx(ij(j)));
    end
    y=y+y1*yy(i)/prod(xx(i)-xx(ij));
end
end

% Command Window 输入
% x=sym('x');
% xx=[4:8];
% yy=[2 2.23607 2.44949 2.64575 2.82843];
% y = lagrange_(xx,yy,x)
% y = lagrange_(xx,yy,7.5)

代码3

function y=lagrangey(xx,yy,x)
n = length(xx);
m = length(x);
y = zeros(1,m);
c1 = ones(n-1,1);
c2 = ones(1,m);
for i=1:n
    xp = xx([1:i-1 i+1:n]);
    y = y + yy(i)*prod((c1*x-xp'*c2)./(xx(i)-xp'*c2));
end
end

% Command Window 输入
% x=sym('x');
% xx=[4:8];
% yy=[2 2.23607 2.44949 2.64575 2.82843];
% y = lagrangey(xx,yy,x)
% y = lagrangey(xx,yy,7.5)

代码4

function y = lagrange_fun(x,xi,yi)
y=0;
n=length(xi);
for i=1:n
    t=1;
    for j=1:n
        if j~=i
            t=t*(x-xi(j))/(xi(i)-xi(j));
        end
    end
    y=y+yi(i)*t;
end
end

% Command Window 输入
% x=sym('x');
% xi=[4:8];
% yi=[2 2.23607 2.44949 2.64575 2.82843];
% y = lagrange_fun(x,xi,yi)
% y = lagrange_fun(7.5,xi,yi)

Newton插值

差商表代码

function A = Chashang_List(a,b)
%Newton插商计算函数
if length(a)~=length(b)
    error('输入数据不同维')
end
N=length(a);
A=zeros(N,N+1);
A(:,1)=a;
A(:,2)=b;
for j=3:N+1
    for i=j-1:N
        A(i,j)=(A(i,j-1)-A(i-1,j-1))/(A(i,1)-A(i-(j-2),1));
    end
end
ChaShang_List=A
end

Newton差商表

x	f(x)	一阶差商	二阶差商	三阶差商	四阶插商
4	2	0	0	0	0
5	2.23607	0.23607	0	0	0
6	2.44949	0.21342	-0.011325	0	0
7	2.64575	0.19626	-0.00858	0.000915	0
8	2.82843	0.18268	-0.00679	0.000597	-0.00008

插值脚本

代码5

function y = Newton(a,b,x,n)
% Newton插值法 a,b分别为x以及对应的f(x),x为插入数据，n为插值多项式次数
% a=4:8
% b=[2 2.23607 2.44949 2.64575 2.82843]
if n>=length(a)
    error('n值太大')
end
A=Chashang_List(a,b);
S=A(1,2);
X=cumprod(x-a);
for i=1:n
    S=S+A(i+1,i+2)*X(i);
end
y=S;
end

代码6

function p= Newton_fun(x,xi,yi)
n=length(xi);
f=zeros(n,n);

% 对差商表第一列赋值
for k=1:n
    f(k)=yi(k);
end
% 求差商表
for i=2:n       % 差商表从0阶开始；但是矩阵是从1维开始存储！！！！！！
    for k=i:n
        f(k,i)=(f(k,i-1)-f(k-1,i-1))/(xi(k)-xi(k+1-i));
    end
end
disp('差商表如下：');
disp(f);

%求插值多项式
p=0;
for k=2:n
    t=1;
    for j=1:k-1
        t=t*(x-xi(j));
        disp(t)
    end
    p=f(k,k)*t+p;
    disp(p)
end
p=f(1,1)+p;

end

% Command Window 输入
% x=sym('x');
% xi=[4:8];
% yi=[2 2.23607 2.44949 2.64575 2.82843];
% p= Newton_fun(x,xi,yi)
% p= Newton_fun(7.5,xi,yi)

Lagrange与Newton插值公式分别得到的多项式值为

Lagrange插值	2.73863
Newton插值	2.73863

问题2

设 $\text{[math]}$ ，试用 $\text{[math]}$

确定二点三次Hermite插值多项式H3(x)并计算H3(pi/12)的值。

构建hermite函数

function y=hermite(x0,y0,dy,x)
n=length(x0);
m=length(x);
for k=1:m 
    yy=0.0;
    for i=1:n
        h=1.0;
        a=0.0;
        for j=1:n
            if j~=i
                h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;
                a=1/(x0(i)-x0(j))+a;
            end
        end
        yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-dy(i))+y0(i));
    end
    y(k)=yy;
end
end

插值多项式H3(x)并计算H3(pi/12)的值

clc;clear;
x=sym('x');
x0=[0,pi/6];
y0=[0,1/2]; %函数值
y1=[1,3^0.5/2]; %导数值
format long;
y=hermite(x0,y0,y1,pi/12)
y=hermite(x0,y0,y1,x)

画图与标准值对比

x1=[0:pi/60:pi/6];y=hermite(x0,y0,y1,x1);
plot(x1,y,'b') %蓝色为插值计算结果
y2=sin(x1);  %标准值
hold on
plot(x1,y2,'r--') %红色为标准值