matlab求解插值与拟合问题(Lagrange插值,chebyshev插值,hermite插值,三次样条插值,曲线拟合)

1.Lagrange插值

function y0 = lagrange(x,y,x0)
%x:插值节点的横坐标
%y:插值节点的纵坐标
%x0:估值点的横坐标,可以是一个数组
%y0:估计出的x0处的函数值,是和x0同型的数组

%x,y转为向量
x=x(:)';
y=y(:)';

%容错处理开始,判断x与y是否同型,x是否包含重复的值
N=length(x);
if N~=length(y)
    error('x和y不同型!');
end
if length(unique(x))~=N&&min(diff(sort(x)))<eps %先判断有没有完全相同的,如果没有,再判断有没有差非常非常小的
    error('x包含相同的元素!');
end
%容错处理结束

sx=size(x0);
x0=x0(:);
ii=1:N;
lenx0=length(x0);
tmp=zeros(lenx0,N);
t=repmat(x0,1,N-1);
for i=1:N
    indx=find(i~=ii);
    %插值公式,prod为连乘函数
    tmp(:,i)=y(i)*prod(t-repmat(x(indx),lenx0,1),2)./...
        prod(repmat(x(i),lenx0,N-1)-repmat(x(indx),lenx0,1),2);
end

y0=sum(tmp,2);
y0=reshape(y0,sx)
end


clear;
clc;
x=[1,1.2,1.8,2.5,4];
y=[0.8415,0.9320,0.9738,0.5985,-0.7568];
y0=lagrange(x,y,1.6);

2.chebyshev插值

clear;
clc;
obj=@exp;
x=-2:.1:2;
y0=obj(x);
sx1=chebyshev2(obj,-2,2,1)
sx2=chebyshev2(obj,-2,2,2)
sx3=chebyshev2(obj,-2,2,3)
y1=polyval(fliplr(sx1),x);
y2=polyval(fliplr(sx2),x);
y3=polyval(fliplr(sx3),x);
plot(x,y0);
hold on
plot(x,y1,'k-');
plot(x,y2,'b--');
plot(x,y3,'m--');
legend('exp(x)','1-order','2-order','3-order')
function che_co = che_gen(n)
che_co=cell(1,n+1);
che_co{1}=1;
che_co{2}=[0,1];
for i=3:n+1
    che_co{i}=zeros(1,i);
    che_co{i}(2:i)=che_co{i-1}*2;
    che_co{i}(1:i-2)=che_co{i}(1:i-2)-che_co{i-2};
end
end
function y=che_i(x,i,che_co)
cof=che_co{i+1};
y=zeros(size(x));
for j=1:i+1
    y=y+cof(j).*x.^(j-1);
end
end
function [s,sx] = chebyshev1(obj,n)
%s:chebyshev多项式系数
%sx:1,x,x^2…的系数
che_co=che_gen(n);
a=zeros(n+1,1);
for i=1:n+1
    a(i)=quad(@(x)fun_phi(x,i-1,obj,che_co),-1,1)*2/pi;
end
s=a;
s(1)=s(1)/2;
sx=zeros(n+1,n+1);
for i=1:n+1
    sx(i,1:i)=che_co{i}*s(i);
end
sx=sum(sx)';

end

function sx = chebyshev2(obj,a,b,n)
%obj:要逼近的函数
%[a,b]:obj的定义域,a默认为-1,b缺省值为1,n缺省值为3
if ~exist('a','var')
    a=-1;
end
if ~exist('b','var')
    b=1;
end
if ~exist('a','var')
    n=3;
end

%若定义区间不是[-1,1],则进行变量代换
g=@(t)obj((t+1)./2*(b-a)+a);

%进行逼近
[~,sx]=chebyshev1(g,n);

%系数变换
co=[2./(b-a),-2*a./(b-a)-1];

%t=cx+d
c=co(1);
d=co(2);

%((sx(n)(cx+d)+sx(n-1))(cx+d)+sx(n-2))(cx+d)++sx(1)
tmp=sx(end);
for i=length(sx):-1:2
    tmp=conv(tmp,[d,c]);
    tmp(1)=tmp(1)+sx(i-1);
end
sx=tmp;
end
function y= fun_phi(x,i,obj,che_co)
y=weightx(x).*che_i(x,i,che_co).*obj(x);
end
function y=weightx(x)
y=1./sqrt(1-x.^2);
end

3.hermite插值

clear;
clc;
%给定样本点的值
x=0:.125:.5;
y=[0,6.24,7.75,4.85,0];
%给定待插值点
xi=0:.01:.5;
%设置随机数种子为0
rand('state',0);
dy=rand(1,5)-0.5;
y0=hermite(x,y,dy,xi);
plot(xi,y0,'r-');
set(gcf,'color','w');
hold on;
plot(x,y,'ob')
function y0 = hermite(x,y,dy,x0)
%x:插值节点的横坐标
%y:插值节点的纵坐标
%dy:x处的一阶导数值
%x0:估值点的横坐标,可以是一个数组
%y0:估计出的x0处的函数值,是和x0同型的数组

%x,y转为向量
x=x(:)';
y=y(:)';

%如果dy=[],利用x和y来计算出dy
if isempty(dy)
    dy=gradient(y,x);
end

%容错处理开始,判断x与y是否同型,x是否包含重复的值
N=length(x);
if N~=length(y)||N~=length(dy)
    error('x、y、dy不同型!');
end
if length(unique(x))~=N&&min(diff(sort(x)))<eps %先判断有没有完全相同的,如果没有,再判断有没有差非常非常小的
    error('x包含相同的元素!');
end
%容错处理结束

sx=size(x0);
x0=x0(:);
lenx0=length(x0);
%迭代开始
ii=1:N;
l=zeros(lenx0,N);
a=l;
y0=l;
t=repmat(x0,1,N-1);

for i=1:N
    indx=find(i~=ii);
    xt=repmat(x(i),lenx0,N-1);
    %基函数
    l(:,i)=prod(t-repmat(x(indx),lenx0,1),2)./...
        prod(xt-repmat(x(indx),lenx0,1),2);
    a(:,i)=sum(1./(xt-repmat(x(indx),lenx0,1)),2);
    y0(:,i)=l(:,i).^2.*((1-2*(x0-repmat(x(i),lenx0,1)).*a(:,i)).*y(i)+(x0-repmat(x(i),lenx0,1))*dy(i));
end
y0(:,i)=l(:,i).^2.*((1-2*(x0-repmat(x(i),lenx0,1)).*a(:,i)).*y(i)+(x0-repmat(x(i),lenx0,1))*dy(i));
y0=sum(y0,2);
y0=reshape(y0,sx)

end

4.三次样条插值

clear;
clc;
x=0:8;
y=[5,9,8,17,26,30,25,39,38];
x0=0:.2:8;
y0=sanci(x,y,x0)
plot(x,y,'o');
hold on;
plot(x0,y0,'r-')
set(gcf,'color','w');
function [x,y,sizey,endslopes] = chckxy(x,y)
if length(find(size(x)>1))>1 
  error(message('MATLAB:chckxy:XNotVector')) 
end

if any(~isreal(x)) 
  error(message('MATLAB:chckxy:XComplex')) 
end

nanx = find(isnan(x));
if ~isempty(nanx)
   x(nanx) = [];
   warning(message('MATLAB:chckxy:nan'))
end

n=length(x);
if n<2 
  error(message('MATLAB:chckxy:NotEnoughPts')) 
end

x=x(:).'; 
dx = diff(x);

if any(dx<0), [x,ind] = sort(x); dx = diff(x); else ind=1:n; end

if ~all(dx), error(message('MATLAB:chckxy:RepeatedSites')), end

sizey = size(y);


while length(sizey)>2&&sizey(end)==1, sizey(end) = []; end


yn = sizey(end); 
sizey(end)=[]; 
yd = prod(sizey);

if length(sizey)>1
   y = reshape(y,yd,yn);
else
   if yn==1
       yn = yd;
       y = reshape(y,1,yn); 
       yd = 1; 
       sizey = yd;
   end
end

nstart = n+length(nanx);
if yn==nstart
   endslopes = [];
elseif nargout==4&&yn==nstart+2
   endslopes = y(:,[1 n+2]); y(:,[1 n+2])=[];
   if any(isnan(endslopes))
      error(message('MATLAB:chckxy:EndslopeNaN'))
   end
   if any(isinf(endslopes))
       error(message('MATLAB:chckxy:EndslopeInf'))
   end
else
   error(message('MATLAB:chckxy:NumSitesMismatchValues',nstart, yn))
end

if ~isempty(nanx)
    y(:,nanx) = [];
end

y=y(:,ind);
nany = find(sum(isnan(y),1));
if ~isempty(nany)
   y(:,nany) = []; x(nany) = [];
   warning(message('MATLAB:chckxy:IgnoreNaN'))
   n = length(x);
   if n<2 
     error(message('MATLAB:chckxy:NotEnoughPts')) 
   end
end
function y0 = sanci(x,y,x0)
[x,y,sizey,endslopes] = chckxy(x,y);
n = length(x); 
yd = prod(sizey);
dd = ones(yd,1); 
dx = diff(x); 
divdif = diff(y,[],2)./dx(dd,:); 
if n==2
   if isempty(endslopes)
      pp=mkpp(x,[divdif y(:,1)],sizey);
   else         
      pp = pwch(x,y,endslopes,dx,divdif);
      pp.dim = sizey;
   end
elseif n==3&&isempty(endslopes)
   y(:,2:3)=divdif;
   y(:,3)=diff(divdif')'/(x(3)-x(1));
   y(:,2)=y(:,2)-y(:,3)*dx(1); 
   pp = mkpp(x([1,3]),y(:,[3 2 1]),sizey);
else
   b=zeros(yd,n);
   b(:,2:n-1)=3*(dx(dd,2:n-1).*divdif(:,1:n-2)+dx(dd,1:n-2).*divdif(:,2:n-1));
   if isempty(endslopes)
      x31=x(3)-x(1);xn=x(n)-x(n-2);
      b(:,1)=((dx(1)+2*x31)*dx(2)*divdif(:,1)+dx(1)^2*divdif(:,2))/x31;
      b(:,n)=...
      (dx(n-1)^2*divdif(:,n-2)+(2*xn+dx(n-1))*dx(n-2)*divdif(:,n-1))/xn;
   else
      x31 = 0; xn = 0; b(:,[1 n]) = dx(dd,[2 n-2]).*endslopes;
   end
   dxt = dx(:);
   c = spdiags([ [x31;dxt(1:n-2);0] ...
        [dxt(2);2*(dxt(2:n-1)+dxt(1:n-2));dxt(n-2)] ...
        [0;dxt(2:n-1);xn] ],[-1 0 1],n,n);

   mmdflag = spparms('autommd');
   spparms('autommd',0);
   s=b/c;
   spparms('autommd',mmdflag);

   pp = pwch(x,y,s,dx,divdif); pp.dim = sizey;

end

if nargin==2
    y0 = pp; 
else
    y0 = ppval(pp,x0); 
end

5.曲线拟合

function [co,error] = squarefit(x,y,phifun,weight)
%x:给定点的横坐标
%y:给定点的纵坐标
%phifun:Ax=b中的x,此处即(1,x,x^2)
%weight:给定点的权重
%co:基函数的系数
%error:拟合函数在各点的误差
if ~exist('weight','var')
    weight=ones(size(x));
end

x=x(:)';
y=y(:)';

N=length(phifun);

%求法方程组
A=zeros(N);
b=zeros(1,N);
for i=1:N
    fi=phifun{i};
    for j=1:N
        fj=phifun{j};
        A(i,j)=sum(weight.*fi(x).*fj(x));
    end
    b(i)=sum(weight.*fi(x).*y);
end

%解法方程组
co=A\b';

%求拟合误差
ys=zeros(1,length(x));
for i=1:N
    ys=ys+co(i)*phifun{i}(x);
end
error=y-ys;
end
clear;
clc;
x=[19,25,31,38,45];
y=[19.5,32.3,49,73,98.7];
f1={@(x)ones(size(x)),@(x)x,@(x)x.^2};
[co,error]=squarefit(x,y,f1)
ys=co(1)*f1{1}(x)+co(2)*f1{2}(x)+co(3)*f1{3}(x);
plot(x,y,'o');
hold on;
plot(x,ys,'r-');
legend('数据点','二次拟合')
title('二次拟合');
set(gcf,'color','w');
figure,plot(x,y-ys)
axis([15,45,-3,3])
grid on;
title('二次拟合的残差');
set(gcf,'color','w');
  • 3
    点赞
  • 40
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

春野与望

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值