文章详细介绍了多种数学插值方法，包括待定系数法、拉格朗日法、牛顿插值法、分段线性插值、Hermite插值和三次样条插值，以及它们在MATLAB中的实现。同时，文章探讨了拟合技术，如最小二乘法和多项式拟合，展示了如何在实际问题中应用这些方法进行数据处理和分析。
- 插值
【1】待定系数法
% inter_polyval.m
clc,clear
x0=[1:6]';
y0=[16,18,21,17,15,12]';
A=vander(x0),p=A\y0
x=[1.5,2.6];
yh=polyval(p,x)
插值多项式系数
p =
-0.2417
4.3333
-28.9583
87.6667
-115.8000
69.0000
x=1.5和2.6时的y值
yh =
14.9180 20.8846【2】拉格朗日法
% inter_lagrange.m
clc,clear
x0=[1:6]';
y0=[16,18,21,17,15,12]';
x=[1.5,2.6];
yh=lagrange(x0,y0,x)
% lagrange.m
function y = lagrange(x0,y0,x)
n=length(x0);
m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
end
x=1.5和2.6时的y值
yh =
14.9180 20.8846
【3】牛顿插值法
% inter_newton.m
X = [0 pi/6 pi/4 pi/3 pi/2];
Y = [0 0.5 0.7071 0.8660 1];
x = linspace(0,pi,50);
hold on;
grid on;
M = 1;
[y,R,A,C,L] = newton(X, Y, x, M);
plot(X,Y,'r*');
plot(x,y,'b--');
legend('原始数据','牛顿插值');
% 求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
% 输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量，Y是纵坐标向量，
% x是以向量形式输入的m个插值点，M在[a,b]上满足｜f~(n+1)(x)｜≤M
% 注：f~(n+1)(x)表示f(x)的n+1阶导数
% 输出的量：向量y是向量x处的插值，误差限R，n次牛顿插值多项式L及其系数向量C，
% 差商的矩阵A
function[y,R,A,C,L] = newton(X,Y,x,M)
n = length(X);
m = length(x);
for t = 1 : m
z = x(t);
A = zeros(n,n);
A(:,1) = Y';
s = 0.0; p = 1.0; q1 = 1.0; c1 = 1.0;
for j = 2 : n
for i = j : n
A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));
end
q1 = abs(q1*(z-X(j-1)));
c1 = c1 * j;
end
C = A(n, n); q1 = abs(q1*(z-X(n)));
for k = (n-1):-1:1
C = conv(C, poly(X(k)));
d = length(C);
C(d) = C(d) + A(k,k);%在最后一维，也就是常数项加上新的差商
end
y(t) = polyval(C,z);
R(t) = M * q1 / c1;
end
L = poly2sym(C);
end
【4】分段线性插值
可能出现龙格振荡
% inter_piecelinear.m
% 在区间[-5,5]上，用n+1个等距节点做多项式，使得节点处的值与y=1/(1+x.^2)在对应节点处的值相等。
% 考察n=6,8,10时，多项式的次数与逼近误差的关系
clc,clear,close all
global yx x
yx=@(x)1./(1+x.^2);
x=linspace(-5,5,100);
set(gca,'Fontsize',15);
hold on,fun(6,1),fun(8,2),fun(10,3)
fplot(yx,[-5,5],'LineWidth',1.5)
legend({'$n=6$','$n=8$','$n10$','$y=1/(1+x^2)$'},...
'Interpreter','latex','Location','north')
function fun(n,i)
global yx x
s={'--*k','-.k','-pk'};
x0=linspace(-5,5,n+1);
y0=yx(x0);
y=lagrangefun(x0,y0,x);
plot(x,y,s{i})
end
【5】Hermite 插值
如果对插值函数，不仅要求它在 节 点 处 与 函 数 同 值在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值，这就是 Hermite 插值问题。
% inter_hermite.m
%{
设n 个节点的数据以
数组 x0 （已知点的横坐标）
y0 （函数值）
y1（导数值）
输入（注意 Matlat 的数组下标从 1 开始）
m 个插值点以数组 x 输入
输出数组 y 为m个插值
%}
function y=inter_hermite(x0,y0,y1,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)-y1(i))+y0(i));
end
y(k)=yy;
end
【6】三次样条插值
许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续，而且要有连续的曲率，这就导致了样条插值的产生。
% inter_spline3.m
clc,clear,format long g
t0=0.15:0.01:0.18;
v0=[3.5 1.5 2.5 2.8];
pp=csape(t0,v0)
xishu=pp.coefs
s1=integral(@(t)fnval(pp,t),0.15,0.18)
vt=griddedInterpolant(t0,v0,'spline')
s2=integral(@(t)vt(t),0.15,0.18)
format
xishu =
-616666.666666667 33500 -473.333333333334 3.5
-616666.666666667 15000 11.6666666666671 1.5
-616666.666666668 -3499.99999999999 126.666666666667 2.5
s1 =
0.068625
vt =
griddedInterpolant - 属性:
GridVectors: {[0.15 0.16 0.17 0.18]}
Values: [3.5 1.5 2.5 2.8]
Method: 'spline'
ExtrapolationMethod: 'spline'
s2 =
0.068625应用：欧洲国土面积
% inter_spl3_Euro.m
clc,clear,close all
a=load('inter_spl3_Euro.txt');
x0=a(1:3:end,:);x0=x0';x0=x0(:);
y1=a(2:3:end,:);y1=y1';y1=y1(:);
y2=a(3:3:end,:);y2=y2';y2=y2(:);
plot(x0,y1,'*-'),hold on,plot(x0,y2,'.-')
%spline3
pp1=csape(x0,y1);pp2=csape(x0,y2);
dp1=fnder(pp1);dp2=fnder(pp2); %求三次样条插值函数的导数
L1=integral(@(x)sqrt(1+fnval(dp1,x).^2)+...
sqrt(1+fnval(dp2,x).^2),x0(1),x0(end)); %地图弧长
L2=L1/18*40 %边界实际长度
Smap=integral(@(x)fnval(pp2,x)-fnval(pp1,x),x0(1),x0(end));%地图上面积
Strue=Smap/18^2*1600 %实际面积
delta=(Strue-42188)/42188 %相对误差，已知精确面积42188
%直接数值积分
borderlow=trapz(x0,sqrt(1+gradient(y1,x0).^2)); %下边界长度
borderhigh=trapz(x0,sqrt(1+gradient(y2,x0).^2));%上边界长度
bordermap=borderlow+borderhigh; %地图边界长度
bordertrue=bordermap/18*40 %实际边界长度
sfake=trapz(x0,y2-y1); %地图近似面积
strue=sfake/18^2*1600 %换算实际面积
%spline3
L2 =
1.1610e+03
Strue =
4.2476e+04
delta =
0.0068
%数值积分
bordertrue =
968.2873
strue =
4.2414e+04
【7】比较
% inter_cmp.m
clc,clear
x0 = [0 3 5 7 9 11 12 13 14 15];
y0 = [0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x = 0:0.1:15;
% 调用前面编写的Lagrange插值函数
y1 = lagrange(x0,y0,x);
% 分段线性插值
y2 = interp1(x0,y0,x);
% 三次样条插值
y3 = interp1(x0,y0,x,'spline');
% 边界为Largrange的三次样条插值
pp1 = csape(x0,y0); y4=ppval(pp1,x);
% 边界为二阶导数的三次样条插值
pp2 = csape(x0,y0,'second'); y5 = ppval(pp2,x);
fprintf('比较一下不同插值方法和边界条件的结果:\n')
fprintf('x y1 y2 y3 y4 y5\n')
xianshi=[x',y1',y2',y3',y4',y5'];
fprintf('%f\t%f\t%f\t%f\t%f\t%f\n',xianshi')
subplot(2,3,1), plot(x0,y0,'+',x,y1), title('Lagrange')
subplot(2,3,2), plot(x0,y0,'+',x,y2), title('Piecewise linear')
subplot(2,3,3), plot(x0,y0,'+',x,y3), title('Spline1')
subplot(2,3,4), plot(x0,y0,'+',x,y4), title('Spline2')
subplot(2,3,5), plot(x0,y0,'+',x,y5), title('Spline3')
dyx0=ppval(fnder(pp1),x0(1)) %求x=0处的导数
ytemp=y3(131:151);
index=find(ytemp==min(ytemp));
xymin=[x(130+index),ytemp(index)]
【6】二维插值
（1）样条插值
% inter2_demo.m
%{
例2 在一丘陵地带测量高程，x 和 y 方向每隔100米测一个点，得高程如2表，试插
值一曲面，确定合适的模型，并由此找出最高点和该点的高程。
%}
% 二维插值
clear,clc
x=100:100:500;
y=100:100:400;
z=[636 697 624 478 450
698 712 630 478 420
680 674 598 412 400
662 626 552 334 310];
% hold on;
xi = 100:10:500;yi=100:10:400
% 三次样条插值注意得到的cz1需要到转一下
pp = csape({x,y},z')
cz1 = fnval(pp,{xi,yi})
% 二维插值
cz2 = interp2(x,y,z,xi,yi','spline')
surfl(x,y,z-400);
hold on;
surfl(xi,yi,cz1');
surfl(xi,yi,cz2+400);
[i,j] = find(cz1==max(max(cz1)))
x = xi(i),y=yi(j),zmax=cz1(i,j)
i =
8
j =
9
x =
170
y =
180
zmax =
720.6252
（2）网格节点插值
% inter2_netnode.m
clc,clear,close all,z0=load('inter2_netdata.txt');
x0=0:100:1400;y0=1200:-100:0;
x=0:10:1400;y=[1200:-10:0]';
pp=csape({x0,y0},z0');
z=fnval(pp,{x,y});z=z';
max_z=max(max(z)),min_z=min(min(z))
subplot(121),h=contour(x,y,z);clabel(h)
subplot(122),mesh(x,y,z)
m=length(x);n=length(y);s=0;
for i=1:m-1
for j=1:n-1
p1=[x(i),y(j),z(j,i)];
p2=[x(i+1),y(j),z(j,i+1)];
p3=[x(i+1),y(j+1),z(j+1,i+1)];
p4=[x(i),y(j+1),z(j+1,i)];
p12=norm(p1-p2);p23=norm(p3-p2);p13=norm(p3-p1);
p14=norm(p4-p1);p34=norm(p4-p3);
z1=(p12+p23+p13)/2;s1=sqrt(z1*(z1-p12)*(z1-p23)*(z1-p13));
z2=(p13+p14+p34)/2;s2=sqrt(z2*(z2-p13)*(z2-p14)*(z2-p34));
s=s+s1+s2;
end
end
s
（3）散乱数据插值
% inter2_scatter1.m
clc,clear,close all
x0=[129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5];
y0=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];
z0=-[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
xm0=minmax(x0)
ym0=minmax(y0)
xi=xm0(1):xm0(2);
yi=ym0(1):ym0(2);
Fz=scatteredInterpolant(x0',y0',z0','natural','nearest')
[xi,yi]=meshgrid(xi,yi);
zi=Fz(xi,yi);mesh(xi,yi,zi)
% inter2_scatter2.m
x=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5];
y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
z=-[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
xi=75:1:200;
yi=-50:1:150;
zi=griddata(x,y,z,xi,yi','cubic');
subplot(1,2,1), plot(x,y,'*')
subplot(1,2,2), mesh(xi,yi,zi)
2.拟合
【1】解方程法
% fit_equal.m
x=[19 25 31 38 44]';
y=[19.0 32.3 49.0 73.3 97.8]';
r=[ones(5,1),x.^2];
ab=r\y;
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')
【2】最小二乘法
% fit_leastsquare_1.m
clc,clear,t=[0:7]';
y=[27.0,26.8,26.5,26.3,26.1,25.7,25.3,24.8]';
tb=mean(t);yb=mean(y);
ahat=sum((t-tb).*(y-yb))/sum((t-tb).^2)
bhat=yb-ahat*tb
ahat =
-0.3036
bhat =
27.1250% fit_leastsquare_2.m
clc,clear
x0=[5.764 6.286 6.759 7.168 7.408]';
y0=[0.648 1.202 1.823 2.526 3.360]';
a=[x0.^2,x0.*y0,y0.^2,x0,y0];
b=-ones(5,1);c=a\b
fxy=@(x,y)c(1)*x.^2+c(2)*x.*y+c(3)*y.^2+c(4)*x+c(5)*y+1;
h=fimplicit(fxy,[3.5,8,-1,5]),title('')
set(h,'Color','k','LineWidth',1.5)
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex','Rotation',0)
c =
0.0508
-0.0702
0.0381
-0.4531
0.2643
【3】最小二乘优化
（1）应用场景
在最小二乘意义下解约束线性方程组
（2）lsqlin
x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)
% fit_ls_lsqlin.m
clc,clear
x0=[3 5 6 7 4 8 5 9]';
d=[4 9 5 3 8 5 8 5]';
C=[exp(x0),log(x0)];
A=[1 1];b=1;
lb=zeros(2,1);
cs=lsqlin(C,d,A,b,[],[],lb)
cs =
0.0005
0.9995（3）lsqcurvefit
X=LSQCURVEFIT(FUN,X0,XDATA,YDATA,LB,UB,OPTIONS)其中 FUN 是定义函数F(x, xdata)的 M 文件。
td=100:100:1000;
cd=[4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59];
x0=[0.2 0.05 0.05];
x=lsqcurvefit(@fun1,x0,td,cd)
function f=fun1(x,tdata);
f=x(1)+x(2)*exp(-0.02*x(3)*tdata); %其中 x(1)=a，x(2)=b，x(3)=k
end
x =
-85.9771 90.3988 -0.0014clc,clear,close all,rng(1)
x0=linspace(-6,6,50);y0=linspace(-5,5,60);
[X0,Y0]=meshgrid(x0,y0);
XX0=X0(:);YY0=Y0(:);
zxy=@(c,xy)exp(-((xy(:,1)-c(1)).^2+(xy(:,2)-c(2)).^2)/2/c(3)^2);
Z0=zxy([1,4,6],[XX0,YY0])+normrnd(0,1,size(XX0));
c1=lsqcurvefit(zxy,rand(1,3),[XX0,YY0],Z0)
f=fittype('exp(-((x-mu1)^2+(y-mu2)^2)/2/s^2)',...
'independent',{'x','y'},'dependent','z');
[c2,gof]=fit([XX0,YY0],Z0,f,'StartPoint',rand(3,1))
f=@(x,y,c)exp(-((x-c(1)).^2+(y-c(2)).^2)/2/c(3)^2);
subplot(121);fsurf(@(x,y)f(x,y,[1,4,6]))
subplot(122);fsurf(@(x,y)f(x,y,c1))
c1 =
0.7811 4.1812 6.1881
（4）lsqnonlin
X=LSQNONLIN(FUN,X0,LB,UB,OPTIONS)其中 FUN 是定义向量函数F(x)的 M 文件。
x0=[0.2 0.05 0.05]; %初始值是任意取的
x=lsqnonlin(@fun2,x0)
function f=fun2(x);
td=100:100:1000;
cd=[4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59];
f=x(1)+x(2)*exp(-0.02*x(3)*td)-cd;
end
x =
-85.9771 90.3988 -0.0014（5）lsqnonneg
X = LSQNONNEG(C,d,X0,OPTIONS)
c=[0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170];
d=[0.8587;0.1781;0.0747;0.8405];
x=lsqnonneg(c,d)
x =
0
0.6929【4】多项式polyfit
[P,S,mu]=polyfit(X,Y,m)其中输入参数 x0,y0 为要拟合的数据，m 为拟合多项式的次数，输出参数 p 为拟合多项式。根据样本数据X和Y，产生一个m次多项式P及其在采样点误差数据S，m为几次拟合，mu是一个二元向量，mu(1)是mean(X),而mu(2)是std(X).
% fit_ls_polyfit.m
x0=[1990 1991 1992 1993 1994 1995 1996];
y0=[70 122 144 152 174 196 202];
plot(x0,y0,'*')
hold on;
a=polyfit(x0,y0,1)
x = 1990:0.1:2000;
y = polyval(a,x);
plot(x,y,'r--');
y97=polyval(a,1997)
y98=polyval(a,1998)
a =
1.0e+04 *
0.0020 -4.0705
y97 =
233.4286
y98 =
253.9286
【5】曲线拟合与函数逼近
% fit_funapproach.m
clc,clear
syms x,base=[1,x^2,x^4];
y1=base.'*base,y2=cos(x)*base.'
r1=int(y1,-pi/2,pi/2)
r2=int(y2,-pi/2,pi/2)
a=r1\r2
xishu1=double(a)
xishu2=vpa(a,6)
xishu1 =
0.9996
-0.4964
0.0372
【6】fittype和fit函数
% fit_fit1.m
% z=a*exp(b*x)+c* y.^2
% fit_fit1.txt
6 2 6 7 4 2 5 9
4 9 5 3 8 5 8 2
14.2077 39.3622 17.8077 11.8310 32.8618 16.9622 33.0941 11.1737
clc,clear,d=load('fit_fit1.txt');
xy0=d([1,2],:)';z0=d(3,:)';
g=fittype('a*exp(b*x)+c*y^2','dependent','z','independent',{'x','y'})
[d,st]=fit(xy0,z0,g,'StartPoint',rand(1,3))
a = 6.193 (5.043, 7.343)
b = 0.04353 (0.01983, 0.06723)
c = 0.3995 (0.3856, 0.4135)利用函数y=2cos2x+6sin2x产生模拟数据，拟合二阶傅里叶级数
% fit_fit2.m
clc,clear
x0=[1:50]';y0=2*cos(2*x0)+6*sin(2*x0);
lb=[-inf*ones(1,5),1];ub=[inf*ones(1,5),1];
[f,g]=fit(x0,y0,'fourier2','Lower',lb,'Upper',ub)
General model Fourier2:
f(x) = a0 + a1*cos(x*w) + b1*sin(x*w) +
a2*cos(2*x*w) + b2*sin(2*x*w)
Coefficients (with 95% confidence bounds):
a0 = 8.246e-17 (-3.478e-16, 5.127e-16)
a1 = -5.024e-16 (-1.112e-15, 1.075e-16)
b1 = -1.982e-16 (-8.052e-16, 4.089e-16)
a2 = 2 (2, 2)
b2 = 6 (6, 6)
w = 1 (fixed at bound)% fit_fit3.m
clc,clear
g=@(a,b,c,d,k,x)(a+b*x).*(x<k)+(c+d*x).*(x>=k);
g=fittype(g)
x=[0.81;0.91;0.13;0.91;0.63;0.098;0.28;0.55;0.96;0.96;0.16;0.97;0.96];
y=[0.17;0.12;0.16;0.0035;0.37;0.082;0.34;0.56;0.15;-0.046;0.17;-0.091;-0.071];
f=fit(x,y,g,'StartPoint',[1,0,1,0,0.5])
plot(f,x,y),legend({'原始数据点','拟合的曲线'})
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex','Rotation',0)
General model:
f(x) = (a+b*x).*(x<k)+(c+d*x).*(x>=k)
Coefficients (with 95% confidence bounds):
a = -0.03779
b = 1.352
c = 1.237
d = -1.301
k = 0.5
y=a+bx , x<k
c+dx , x>=k
【7】GUI
%拟合5阶傅里叶级数与y=2x-1的交点
x=[0.81;0.91;0.13;0.91;0.63;0.098;0.28;0.55;0.96;0.96;0.16;0.97;0.96];
y=[0.17;0.12;0.16;0.0035;0.37;0.082;0.34;0.56;0.15;-0.046;0.17;-0.091;-0.071];
plot(f,x,y),hold on,fplot(@(x)2*x-1,[0.4,1],'--')
fx=@(x)f(x)-(2*x-1)
x=fsolve(fx,rand)
y=f(x)
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','latex','Rotation',0)
x =
0.6777
y =
0.3555
1092
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
