matlab lagrange插值，Newton插值，Hermite插值，spline插值

时间  2020-07-14 标签 matlab lagrange 插值 newton hermite spline

标签（空格分隔）： matlab 插值 算法 数值 分析 lagrange Newton Hermite spline

文章目录

matlab lagrange插值，Newton插值，Hermite插值，spline插值

1 lagrange插值

1.1 测试lagrange插值函数lagrangeInsertValue()与lagrangeInsertValueWithError()

clc,clear,clf
f=@(x)(sqrt(x));
% x=[4 9 6.25 40]
x=[100 121 81 144]
y=f(x)
% xs=[7 8 100]
xs=[105]
% x=[-1 0 .5 1 ]
% y=[3  -.5 0 1]
% syms xs
%绘制插值函数图
hold on
plot(x,y,'o')
min_=min(x);max_=max(x);
xs_pic=min_:(max_-min_)/50:max_;
ys_pic=lagrangeInsertValue(x,y,xs_pic);
plot(xs_pic,ys_pic,'-.')
ys_pic=f(xs_pic);
plot(xs_pic,ys_pic)
xs_pic=[];ys_pic=[];
hold off
%%%%%%%%%%%%计算插值、偏差
% ys=lagrangeInsertValue(x,y,xs)
%增长的用于求偏差的点
xadd=prod(xs)^(1/length(xs));xadd=xadd+.1*max(xs);
yadd=f(xadd);
[ys,rs]=lagrangeInsertValueWithError(x,y,xs,xadd,yadd)
xadd,yadd
rs_true=f(xs)-ys

获得函数曲线与插值曲线以下：
在这里插入图片描述web

1.2 lagrangeInsertValue()

lagrange插值函数，不包含偏差计算功能算法

function ys=lagrangeInsertValue(x,y,xs)
%【程序特色】：考虑将内存和运算量同时达到极小化
%ys为插值xs的函数值
%rs为对应的插值偏差
%时间：2013/12/08 姓名：邓能财   
%%% 标记为显示行（可去）%#符号
len=length(x);
vlen=length(xs);
% %计算xi-xj造成的矩阵
% dx=zeros(len);
% for i=1:len
%     for j=i+1:len
%         dx(i,j)=x(i)-x(j);
%     end
% end
% dx=dx-dx'; %xj-xi=-(xi-xj)
% for i=1:len
%     dx(i,i)=1;
% end
%dx
%%%%%%%%%%%%计算w(xi)
wxi=ones(1,len);
for i=1:len
    for j=1:len
        if j~=i
            wxi(i)=wxi(i)*(x(i)-x(j));
        end
    end
end
wxi
%%%%%%%%%%%%计算函数值
    %用【卷计算的方式节约到vlen个double内存】
cmMultidx=zeros(1,vlen);
for j=1:vlen
    cmMultidx(j)=prod(xs(j)-x);
end
cmMultidx%#
ys=zeros(1,vlen);
sum_l=ys;%#
for i=1:len
    l=cmMultidx./((xs-x(i))*wxi(i))
    sum_l=sum_l+l;%#
    ys=ys+l.*y(i);
end
sum_l%#
end

1.3 lagrangeInsertValueWithError()

lagrange插值函数，包含偏差计算功能svg

function [ys,rs]=lagrangeInsertValueWithError(x,y,xs,xadd,yadd)
%lagrange插值，包含【过后偏差估计】版
%时间：2013/12/08 姓名：邓能财   
ys=lagrangeInsertValue(x,y,xs);
rs=ys-lagrangeInsertValue([x(2:end),xadd],[y(2:end),yadd],xs);
end

2 Newton插值

2.1 测试Newton插值函数newtonInsertValue()

clc,clear,clf
% x=[4 9 6.25]
% y=[2 3 2.5]
% xs=[7 8 8.5]
f=@(x)(x.^5+x.^3+1);
x=[-1 -.8 0 .5 1]
y=f(x)
%绘制插值函数图
hold on
plot(x,y,'o')
min_=min(x);max_=max(x);
xs_pic=min_:(max_-min_)/50:max_;
ys_pic=newtonInsertValue(x,y,xs_pic);
plot(xs_pic,ys_pic,'-.')
ys_pic=f(xs_pic);
plot(xs_pic,ys_pic)
xs_pic=[];ys_pic=[];
hold off
xs=[-.4 .7 4]
[ys,rs]=newtonInsertValue(x,y,xs)
xs

获得函数曲线与插值曲线以下：
在这里插入图片描述函数

2.2 newtonInsertValue()

function [ys,rs]=newtonInsertValue(x,y,xs)
%【程序特色】：考虑将内存和运算量同时达到极小化
%ys为插值xs的函数值
%rs为对应的插值偏差
%时间：2013/12/08 姓名：邓能财   
%%% 标记为显示行（可去）%#符号
len=length(x);
vlen=length(xs);
%%%%%%%%计算各阶差商，在计算差商中只记录：
    %1 插值多项式系数的差商df1
    %2 用于偏差计算的差商dfn
    %3 计算过程当中间的i阶差商dfi
df1=zeros(1,len); dfn=df1; dfi=df1;
dfi=y;
dfi_show=dfi(1:len)%#
df1(1)=dfi(1)
dfn(1)=dfi(len)
for i=1:len-1 %差商阶数
    for j=1:len-i %列标
        dfi(j)=(dfi(j+1)-dfi(j))/(x(j+i)-x(j));
    end
    dfi_show=dfi(1:len-i)%#
    df1(i+1)=dfi(1)
    dfn(i+1)=dfi(len-i)
end
dfi=[]; %释放内存
%%%%%%%%%计算插值多项式的值
    %用【卷计算的方式节约到vlen个double内存】
    %设置卷变量：wxi=(xs-x1)*...*(xs-xi)
wxi=ones(1,vlen);
ys=df1(1)*wxi;
for i=2:len
    wxi=wxi.*(xs-x(i-1));
    ys=ys+df1(i)*wxi;
end
df1=[]; %释放内存
%%%%%%%%%计算偏差项
%计算包含xs的差商df_nx
dfn_x=ys; %df_nx
for i=1:len
    dfn_x=(dfn_x-dfn(i))./(xs-x(len+1-i));
end
if length(dfn_x)<10, dfn_x,end%#
wxi=wxi.*(xs-x(len));
rs=dfn_x.*wxi;
end

3 Hermite插值

3.1 测试Hermite插值函数hermiteInsertValue()

clc,clear,clf
%%%%%%%%e1
% f=@(x)(sqrt(x));
% df=@(x)(.5*x.^(-.5));
% x=[4 9 6.25 40]
% xs=[7 8 100]
%%%%%%%%%e2
f=@(x)(sqrt(x).*sin(x));
df=@(x)(.5*x.^(-.5).*sin(x)+sqrt(x).*cos(x));
x=[4 9 16.5]
y=f(x)
dy=df(x)
xs=[7 8 20]
%绘制插值函数图
hold on
plot(x,y,'o')
min_=min(x);max_=max(x);
h_pic=(max_-min_)/50;
xs_pic=min_-1:h_pic:max_+1;
ys_pic=hermiteInsertValue(x,y,dy,xs_pic);
plot(xs_pic,ys_pic,'-.')
ys_pic=f(xs_pic);
plot(xs_pic,ys_pic)
xs_pic=[];ys_pic=[];
hold off
%%%%%%
disp('%%%%%%%%%%%%计算插值、偏差')
% [ys,rs]=
ys=hermiteInsertValue(x,y,dy,xs)

获得函数曲线与插值曲线以下：
在这里插入图片描述测试

3.2 hermiteInsertValue()

function ys=hermiteInsertValue(x,y,dy,xs)
%【程序特色】：考虑将内存和运算量同时达到极小化
%ys为插值xs的函数值
%rs为对应的插值偏差
%时间：2013/12/08 姓名：邓能财   
%%% 标记为显示行（可去）%#符号
max_show_numb=10;
len=length(x);
vlen=length(xs);
%%%%%%%%%%%%计算w(xi)、sum_dx(xi)
wxi=ones(1,len);
sum_dxi=zeros(1,len);
for i=1:len
    for j=1:len
        if j~=i
            wxi(i)=wxi(i)*(x(i)-x(j));
            sum_dxi(i)=sum_dxi(i)+1/(x(i)-x(j));
        end
    end
end
if len<max_show_numb, wxi,sum_dxi, end%#
%%%%%%%%%%%%计算函数值
    %用【卷计算的方式节约到vlen个double内存】
cmMultidx=zeros(1,vlen);
for j=1:vlen
    cmMultidx(j)=prod(xs(j)-x);
end
if vlen<max_show_numb, cmMultidx, end%#
ys=zeros(1,vlen);
sum_l=ys;%#
for i=1:len
    l=cmMultidx./((xs-x(i))*wxi(i));
    sum_l=sum_l+l;%#
    gx_f=(1-2*(xs-x(i))*sum_dxi(i)).*l.*l;
    hx_f_=(xs-x(i)).*l.*l;
    ys=ys+gx_f.*y(i)+hx_f_.*dy(i);
end
if vlen<max_show_numb, sum_l, end%#
end

4 spline插值

4.1 测试spline插值函数splineInterpolation()

clc,clear;
clf
% f=@(x)(x.^5+x.^3+1);
% f=@(x)(sqrt(x).*sin(x));
f=@(x)(log(x).*sin(x));
% x=[-1 -.8 0 .5 1]
% x=[.25 .30 .39 .45 .53]
% y=[.50 .54 .62 .67 .73]
x=[1 3 5.2 7.7 9]
y=f(x)
% % 绘制插值函数图
hold on
plot(x,y,'o')
min_=min(x);max_=max(x);
h_pic=(max_-min_)/50; d_pic=6;
xs_pic=min_-d_pic*h_pic:h_pic:max_+d_pic*h_pic;
ys_pic=splineInterpolation(x,y,xs_pic);
plot(xs_pic,ys_pic,'-.')
ys_pic=f(xs_pic);
plot(xs_pic,ys_pic)
xs_pic=[];ys_pic=[];
hold off
%%%%%%
xs=[-.4 .7 4]
m0=1,mn=.68;
ys=splineInterpolation(x,y,xs,m0,mn)
xs

获得函数曲线与插值曲线以下：
在这里插入图片描述3d

4.2 splineInterpolation()

function [ys,rs]=splineInterpolation(x,y,xs,mstart,mend)
%【程序特色】：考虑将内存和运算量同时达到极小化
%ys为插值xs的函数值
%rs为对应的插值偏差
%时间：2013/12/08 姓名：邓能财   
%%% 标记为显示行（可去）%#符号
%%%%%%%%%%%%默认端点条件为：两端导数为0
if nargin<4,  mstart=0;mend=0; end
len=length(x);
vlen=length(xs);
%%%%%%%%%%%%将x从小到大排序，y也随之排序
[x,index]=sort(x); y=y(index); index=[];
%%%%%%%%%%%%计算一阶导数m
%计算线性方程的系数
lamda=zeros(1,len-2); mu=lamda; c=lamda;
i=2;
df_memory=(y(i)-y(i-1))/(x(i)-x(i-1)); df=df_memory;
hi_memory=x(i)-x(i-1); hi=hi_memory;
for im=1:len-2
   i=im+1;
   hi_memory=hi; hi=x(i+1)-x(i);
   df_memory=df; df=(y(i+1)-y(i))/hi;
   lamda(im)=hi/(hi_memory+hi); mu(im)=1-lamda(im);
   c(im)=3*(lamda(im)*df_memory+mu(im)*df);
end
lamda,mu,c%#
%求解线性方程，计算m
c(1)=c(1)-lamda(1)*mstart; c(end)=c(end)-mu(end)*mend;
m=zeros(1,len); m(1)=mstart; m(end)=mend;
m(2:end-1)=bandMatrixSolve(2*ones(1,len-2),mu(1:end-1),lamda(2:end),c);
m%#
%%%%%%%%%%%%计算函数值
    %用【卷计算的方式节约到vlen个double内存】
%x已排序
%将xs排序,再查找xs(i)所在的区间[xi-1,xi]
[xs,index]=sort(xs);
ys=zeros(1,vlen);
j=1;
for i=1:vlen
   %从当前位置（非起始）开始查找xs(i)所在的区间
    while j+1<len && xs(i)>x(j+1), j=j+1; end
    %计算ys
    h=x(j+1)-x(j); h3_=1/(h*h*h);
    dx1=xs(i)-x(j); dx2=dx1*dx1;
    dx1_=xs(i)-x(j+1); dx2_=dx1_*dx1_;
    ys(index(i))=h3_*h*(dx2_*y(j)+dx2*y(j+1)...
			             +dx1*dx2_*m(j)+dx1_*dx2*m(j+1))...
	               +h3_*(2*dx1*dx2_*y(j)...
         				 -2*dx1_*dx2*y(j+1));
end
rs=0;
end

4.3 bandMatrixSolve()

function x=bandMatrixSolve(a,b,c,e)
%【程序特色】：考虑将内存和运算量同时达到极小化
%时间：2013/12/08 姓名：邓能财   
%带状矩阵方程组求解
%矩阵A由中a(dim维)、b、c，构成
%等式右边的向量为：e
%算法：追赶法
dim=length(a);
assert( dim==length(b)+1 &&dim==length(c)+1 &&dim==length(e),...
    ['Argument input error:   ',...
    'Matrix dimensions must agree.'])
%%%%%%%%%%%%%%%%%%%LU分解
alfa=zeros(1,dim); beta=zeros(1,dim-1);
%初始值
alfa(1)=a(1);
beta(1)=b(1)/alfa(1);
for i=2:dim-1
   alfa(i)=a(i)-c(i-1)*beta(i-1);
   beta(i)=b(i)/alfa(i);
end
i=i+1;
alfa(i)=a(i)-c(i-1)*beta(i-1);
a=[]; b=[];%释放内存
%%%%%%%%%%%%%%%%%%%解Ly=e
y=zeros(1,dim); 
%初始值
y(1)=e(1)/alfa(1);
for i=2:dim
   y(i)=(e(i)-c(i-1)*y(i-1))/alfa(i);
end
alfa=[]; c=[];%释放内存
%%%%%%%%%%%%%%%%%%%解Ux=y
x=zeros(1,dim); 
%初始值
x(dim)=y(dim);
for i=dim-1:-1:1
   x(i)=y(i)-beta(i)*x(i+1);
end
beta=[]; y=[];%释放内存
end

联系做者 definedone@163.com

endcode

联系我们 沪ICP备13005482号-10