插值与拟合算法 - 知乎 - 程序猿·D·安宇雨 DeepMind

插值与拟合算法 - 知乎

安宇雨 - 随手采集
2023-08-22 17:56:07
随手采集
0000-未整理-等待研究

一：基本概述

适用的问题：已知函数在某区间（域）内若干点处的值，求函数在该区间（域）内其它点处的值

解决方法：先利用已知点拟合出曲线方程，一般是多项式，然后求插值点的值

二：解决方法

可以用范德蒙行列式和克莱姆法则证明出函数的多项式存在且唯一，即插值问题的解唯一存在。

常用的插值方法是Lagrange插值法和Newton插值法

2.1Lagrange插值法

拉格朗日公式：

指的是在节点上给出节点基函数，然后做基函数的线性组合，组合系数为节点函数值的一种插值多项式。

比如对上图进行插值处理，已知图中标注的三点

法1：可以看出这是一个二次函数，所以可以先设一个多项式 a+bx+x2a+bx+x^{2} ,然后代入x1、x2得到两个式子，最后三式联立解未知数即可，但是这种方法计算较复杂而且多项式的次数越高越难解方程

法2：

在图二中画一个图，使之过（x1，1）、（x2,0）、（x3,0）同理图三和图四也是过了如图所示的三个其他的点

由图二可得，用图二的函数在x1处的函数值乘以y1一定会经过（x1，y1）点同理可由图三和图四得

综上的一下式子

就简化了运算，从而得到需要插值函数的方程

2.2matlab插值法

Maple和Matlab都可以进行插值计算，Maple的一维插值计算较为便捷，而MIatlab的二维插值功能较强,还能进行散乱点插值。

本次主要介绍MIatlab的一维和二维插值命令

2.2.1一维插值

一维插值命令是interp1，其基本格式为yi= interp1(x,y,xi,'method')。

x,y为插值点，xi,yi为被插值点和插值结果，x,y和xi,yi通常为向量;‘method'表示插值方法: nearest'-最邻近插值'linear'线性插值'spline'-三次样条插值，'cubic’立方插值，若是不写method则默认是线性插值。

例：

x=0:2:24;
y=[12 9 9 10 18 24 28 27 25 20 18 15 13];
x1=13;
yl=interp1(x,y,x1,'spline')%在x1出插入，方法是三次样条插样
xi=0:1/3600:24;
yi=interp1(x,y,xi,'spline');%在xi处插入，此时xi是一个向量，方法是三次样条插样
plot(x,y,'*',xi,yi)%绘图，在（x，y）点上标志*

生成的图片如下图所示：

用*标志的位置是已知点坐标

例2：比较拉格朗日插值法、线性插值法和三次插值法

function plane
x0=[0 3 5 7 9 11 12 13 14 15 1];
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;
yl=lagrange(x0,y0,x);%拉格朗日法
y2=interp1(x0,y0,x);%线性法
y3=interp1(x0,y0,x,'spline');%三次插值法
subplot(3,1,1);%在一个平面内同时绘制三个图形而且分三行一列摆放
plot(x0,y0,'k+',x,y,'r');
grid;
title('lagrange');
subplot(3,1,2);
plot(x0,y0,'k+',y2,'r');
grid;
title('piecewise linear');
subplot(3,1,3);
plot(x0,y0,'k+',x,y3,'r')
grid;
title('spline');
%拉格朗日函数，因为matlab中没有lagrange告辞插值功能，所以要单独编写这个函数
function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=l:n
p=1.0;
for j=l:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(i));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end

最终得到上图，有图可得，拉格朗日法的波动较大，所以一般不用，而且用三次插值发得到的图形更加圆滑，所以此方法较好

2.2.2二维插值

二维插值命令是interp2，基本格式为zi=interp2(x,y,z,xi,yi,'method').二维插值命令的使用较复杂。

x,y,z为插值点，z可以理解为被插值函数在(x,y)处的值;xi,yi为被插值点,zi为输出的插值结果，可理解为插值函数在(xi,yi）处的值;x,y为向量，xi,yi为向量或矩阵，而z和zi则为矩阵。

‘method'表示插值方法:‘nearest'一最邻近插值，linear-双线性插值，spline'-双三次样条插值，'cubi c'双立方插值，默认是双线性插值。为什么这里都含有“双”，因为每一个z都对应着x和y，所以进行插值时x和y都要进行

例1：

x=1:5;
y=1:3;
temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];%一共有3*5个数字，而且排列顺序为y的最大数是行数，x的最大数是列数，这个一定不能弄倒
figure(1);
mesh(x,y,temps);%绘空间曲面，有很多线条组成，是空间网格图
xi=1:0.2:5;
yi=1:0.2:3;
zi=interp2(x,y,temps,xi,yi','cubic');%按照cubic的方法在三维曲面上插入值，记得对yi进行转置
figure(2);
mesh(xi,yi,zi);%绘制插入的值的曲面
figure(3);
contour(xi,yi,zi,20,'r');%绘制等高线，就是将所有z相等的点连起来，等高线互不相交，20表示绘制多少条等高线，r表示红色
[i,j]=find(zi==min(min(zi)));%最里面的min是每一列的最小值，外层的的min是求矩阵中最小的值
x=xi(j),y=yi(i),zmin=zi(i,j)%在最小值的位置，注意是xi（j），yi（i），行是y，列是x
[i,j]=find(zi==max(max(zi)));%与找最小值同理，这个是找最大值
x=xi(i),y=yi(i),zmax=zi(i,i)

上述程序较复杂，说明如下

(1) interp2中的xi为行向量，而yi'为列向量，其实xi和yi'行列不同即可（别忘了对yi进行转置）

(2) plot3(空间曲线),mesh(空间曲面)surf (空间曲面)，contour(等高线)是三维作图中的常用命令。mesh和surf的区别是mesh画的是曲面网格图，而surf画的是曲面表面图。

contour(x,y,z,n)的功能是作出由点(x,y,z)插值而成曲面的n条等高线

用meshc和surfc可在曲面下方画等高线，meshz和surfz是画垂帘图

(3) 程序的最后部分为求最高(低)点，注意品味find的功能

例2：在某山区测得一些地点的高程如下表。平面区域为0<x<5600，0<y<4800

试用MIatlab中的最邻近插值、双线性插值和双三次插值3种方法作出该山区的地貌图（三维空间图）和等高线图，并求出最高和最低点

x=0:400:5600;
y=0:400:4800;
z=[370 470 550 600 670 690 670 620 580 450 400 300 100 150 250;...510 620 730 800 850 870 850 780 720 650 500 200 300 350 320;...650 760 880 970 1020 1050 1020 830 900 700 300 500 550 480 350;...740 880 1080 1130 1250 1280 1230 1040 900 500 700 780 750 650 550;...830 980 1180 1320 1450 1420 1400 1300 700 900 850 840 380 780 750;..880 1060 1230 1390 1500 1500 1400 900 1100 1060 950 870 900 930 950;...910 1090 1270 1500 1200 1100 1350 1450 1200 1150 1010 880 1000 1050 1100;...950 1190 1370 1500 1200 1100 1550 1600 1550 1380 1070 900 1050 1150 1200;...1430 1430 1460 1500 1550 1600 1550 1600 1600 1600 1550 1500 1500 1550 1550;...1420 1430 1450 1480 1500 1550 1510 1430 1300 1200 980 850 750 550 500;...1380 1410 1430 1450 1470 1320 1280 1200 1080 940 780 620 460 370 350;...1370 1390 1410 1430 1440 1140 1110 1050 950 820 690 540 380 300 210;...1350 1370 1390 1400 1410 960 940 880 800 690 570 430 290 210 150];
figure(1);
meshz(x,y,z);
xlabel('X'),ylabel('y'),zlabel('z');
[xi,yi]=meshgrid(0:50:5600,0:50:4800);%与x程序用“xi,yi=meshgrid(0:50:560生成网格点(xi,yi)，作用0,0:50:4800);”相当于“xi=0:50:5600;yi=0:50:4800;但meshgrid(x,y)生成的xi,yi为同维矩阵xi的行均为x，而yi的列均为y。
figure(2);
zli=interp2(x,y,z,xi,yi,'nearest');%用最近值的方法插入
surfc(xi,yi,zli);
xlabel('X'),ylabel('y'),zlabel('z');
figure(3);
z2i=interp2(x,y,z,xi,yi);%用线性法插入
surfc(xi,yi,z2i);
xlabel('x'),ylabel('y'),zlabel('z');
figure(4);
z3i=interp2(x,y,z,xi,yi,'cubic');%按照立方插入
[i,j]=find(z3i==min(min(z3i)));%最里面的min是每一列的最小值，外层的的min是求矩阵中最小的值
x=xi(j),y=yi(i),zmin=zi(i,j)%在最小值的位置，注意是xi（j），yi（i），行是y，列是x
[i,j]=find(zi==max(max(z3i)));%与找最小值同理，这个是找最大值
x=xi(j),y=yi(i),zmax=zi(i,i)
surfc(xi,yi,z3i);
xlabel('X'),ylabel('y'),zlabel('z');
%在绘制三个等高线在同一个平面内
figure(5);
subplot(1,3,1),contour(xi,yi,zli,10,'r');
subplot(1,3,2),contour(xi,yi,z2i,10,'r');
subplot(1,3,3),contour(xi,yi,z3i,10,'r');

程序用“[xi,yi]=meshgrid(0:50:5600,0:50:4800);"生成网格点（xi，yi），作用相当于“xi=0:50:5600;yi'=0:50:4800;”但meshgrid(x,y)生成的xi,yi为同维矩阵,xi的行均为x，而yi的列均为y。

生成的图形如下所示：

可得用cubic的方法得到的图形最光滑，最接近事实

2.3散乱点插值

前面讨论的插值问题的插值点(x,y)均为网格点。当(x,y)为散乱点时，可用griddata(x,y,z,xi,yi,'method') 命令进行二维插值

例1：在某海域测得一些点 (x,y)处的水深z如下表，船的吃水深度为5英尺，在矩形区域(75,200)*(-50,150)内的哪些地方船要避免进入。

水深是负值，负值会影响插值的准确率所以最好变成正值，但是不改也行

clear
x=[129 140 103.5 88 185.5 195 105.5 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 356.5 -66.5 84 -33.5];
z=[-4 -8 -6-8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9]
[xi,yi]=meshgrid(75:0.5:200,-70:0.5:150);
zi=griddata(x,y,z,xi,yi,'cubic');%按照cubic的方法插入值
figure(1);
meshz(xi,yi,zi);
xlabel('X'),ylabel'y'),zlabel('z');
figure(2),contour(xi,yi,zi,[-5 -5], 'b');%注意这里的[5,5]表示绘制zi为-5的等高线，且格式只能是[k,k]
grid;
hold on;%在一个平面内再画一个图
plot(x,y,'+');%在等高线的图中再绘制已知点的坐标，用+标注
xlabel('X'),ylabel('y');

得到图形：

在行船时最好在圈外行驶，否则可能会搁浅

三：高次插值的龙格现象

很多人认为在进行插值是，将曲线拟合成越高次的多项式越好，但实际上并不是这样，多项式的次数越高曲线的波动性更大。研究发现当次数超过七时，插值多项式就会出现严重的震荡现象，就称之为龙格现象。

避免龙格现象的常用方法是: 将插值区间分成若干小区间，在小区间内用低次(二次，三次)插值，即分段低次插值，如样条函数插值

总结：一维插值用样条（三次），二维插值用立方

四：拟合问题

1.拟合问题

对于超过已知点范围的点的值预测既不能用高次插值也不能用样条插值，前者是因为龙格现象后者是因为在超过已知点之后的值再仅进行插值波动很大，误差很大。所以这时就要用拟合的方法，找出大致符合已知点的方程，然后进行值的预测。

或者是对函数进行求导的问题，也需要先用拟合的方法求出函数方程然后再求导

这种根据离散数据求数据间近似函数关系的问题称为曲线拟合问题

拟合问题与插值问题的区别在于:

(1)插值函数过已知点，而拟合函数不一定过已知点;

(2) 插值主要用于求函数值，而拟合的主要目的是求函数关系，从而进行预测等进步的分析。

当然，某些特定问题既可以用插值也可以用拟合。

2.拟合的计算

曲线拟合需解决如下两个问题

(1) 线型的选择

(2) 线型中参数的计算线型的选择是拟合计算的关键和难点。

通常主要根据专业知识和散点图确定线型，线性拟合中参数的计算可采用最小二乘法，而非线性拟合参数的计算则要应用Gauss-Newton迭代法。

3.Matlab拟合

(1) 多项式拟合

Matlab多项式拟合命令格式为

[a,S]=polyfit(x,y,n)

其中，x和y是被拟合数据的自变量和因变量(已知点); n为拟合多项式的最高次数; a为拟合多项式系数构成的向量; S为分析拟合效果所需的指标(可省略)

例：

x=1:12;
y=[5, 8, 9, 15, 25, 29, 31, 30, 22, 25, 27, 24];
a=polyfit(x,y,9)%多项式的最高次为9
xp=1:0.1:12;
yp=polyval(a,xp);%利用新生成的多项式求yp
plot(x,y,'.k',xp,yp,'r');%已知点是黑色，利用多项式预测的点是红色的

得到如下的图形：

（2）非线性拟合

Matlab非线性拟合命令格式为:[b,r]=polyfit(x,y,fun,b0,option)

其中，x和y是被拟合数据的自变量和因变量(已知点); fun为拟合函数; b0为拟合参数的初始迭代值（拟合函数的参数的初始值，这个值是假设的）; option为拟合选项; b为拟合参数（b是求出来的参数值）;r为拟合残差

例：

x=1:16;
y=[4.00 6.40 8.00 8.80 9.22 9.50 9.70 9.86 10.0010.20 10.32 10.42 10.50 10.55 10.8 10.60];
yl=@(b,t)b(1)*exp(-t/b(2))+b(3)*exp(-t/b(4))+b(5);%是刚开始根据已知点设的函数，也是上述的fun
b0=[-1 1 -1 1 1];%函数中的参数初始值，随意设置
a=nlinfit(x,y,yl,b0)%求函数的参数值，nlinfit与plotfit的作用一致
xp=1:0.1:16;
yp=y1(a,xp);%利用函数预测未知的点
plot(x,y,'.k',xp,yp,'r');%绘图

图形如下图所示;