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');