一、插值
插值：从已知点近似计算未知点的近似计算方法
１.一维插值函数
y=interp1(x0,y0,x,'method')；其中 x0,y0 是已知数据点，x 是插值点，y 是插值点的函数值
mothod默认为线性插值，其值可为：
‘nearest’ 最近项插值
‘linear’ 线性插值
‘spline’ 三次样条插值 （还可直接spline(x0,y0,x)）
‘cubic’ 立方插值/三次Hermite多项式插值(新版本改为pchip)
注：
①所有的插值方法要求x0是单调的
②被插值节点不能超出x0范围
③若y0为矩阵，则对y0的每一列进行插值，x可以为向量
Matlab 中三次样条插值也有现成的函数：
y=spline(x0,y0,x)；
pp=csape(x0,y0,conds)；
pp=csape(x0,y0,conds，valconds)；y=ppval(pp,x)其中 x0,y0 是已知数据点，x 是插值点，y 是插值点的函数值
注：使用csape函数求插值点的函数值，必须调用函数 ppval。
例：机床加工
待加工零件的外形根据工艺要求由一组数据(x,y) 给出（在平面情况下），用程控铣床加工时每一刀只能沿 x方向和 y 方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的坐标(x,y)。
下表中给出的(x,y)数据位于机翼断面的下轮廓线上，假设需要得到 x坐标每改变 0.1 时的 y 坐标。试完成加工所需数据，画出曲线，并求出 0 =x 处的曲线斜率和 13 ≤x≤ 15 范围内 y 的最小值。
要求用分段线性和三次样条三种插值方法计算。
x | 0 | 3 | 5 | 7 | 9 | 11 | 12 | 13 | 14 | 15 |
y | 0 | 1.2 | 1.7 | 2.0 | 2.1 | 2.0 | 1.8 | 1.2 | 1.0 | 1.6 |
解：
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;
y1=interp1(x0,y0,x);%默认的插值方法：线性插值
y2=interp1(x0,y0,x,'spline');%三次样条插值
pp1=csape(x0,y0);%三次样条插值，使用默认边界条件即Lagrange边界条件
y3=fnval(pp1,x);%求被插值点函数值
pp2=csape(x0,y0,'second');%三次样条插值，边界为二阶导数
y4=fnval(pp2,x);
%画图
subplot(1,4,1)
plot(x0,y0,'*',x,y1)
title('分段线性插值')
subplot(1,4,2)
plot(x0,y0,'*',x,y2)
title('三次样条插值')
subplot(1,4,3)
plot(x0,y0,'*',x,y3)
title('三次样条插值(Lagrange)')
subplot(1,4,4)
plot(x0,y0,'*',x,y4)
title('三次样条插值(Second)')
%因三次样条插值曲线光滑，故采用y3计算斜率
%近似计算每一点的斜率
dx=diff(x);
dy=diff(y3);
dy_dx=dy./dx;
%求x=0处斜率(即dy_dx计算的斜率中的第一个元素)
dy_dx0=dy_dx(1)
%求[13,15]范围内最小值
ytemp=y3(131:151);
ymin=min(ytemp);
index=find(y3==ymin);
xmin=x(index);
[xmin,ymin]图像如图
可以看出，分段线性插值的光滑性较差。
得到计算结果为：
dy_dx0 =
0.4972
ans =
13.8000 0.9851
２.二维插值函数
1）插值节点为网格节点
z=interp2(x0,y0,z0,x,y,'method')其中 x0,y0 分别为m 维和n维向量，表示节点，z0 为 m *n 维矩阵（按照表格原样输入就行），表示节点值，x，y 为一维数组，表示插值点。
注：
①x与y应是方向不同的向量，即一个是行向量，另一个是列向量
②z 为矩阵，它的行数为 y 的维数，列数为 x的维数，表示得到的插值
③'method' 的用法同上面的一维插值。
三次样条插值,还可用命令：
pp=csape({x0,y0},z0,conds,valconds)，z=fnval(pp,{x,y})其中 x0,y0 分别为m 维和n维向量，z0 为 n*m 维矩阵（表格数据需要进行转置），z 为矩阵，它的行数为x的维数，列数为y的维数，表示得到的插值，具体使用方法同一维插值。
例：在一丘陵地带测量高程，x和 y 方向每隔100米测一个点，得高程如2表，试插 值一曲面，确定合适的模型，并由此找出最高点和该点的高程。
y x | 100 | 200 | 300 | 400 | 500 |
100 | 636 | 697 | 624 | 478 | 450 |
200 | 698 | 712 | 630 | 478 | 420 |
300 | 680 | 674 | 598 | 412 | 400 |
400 | 662 | 626 | 552 | 334 | 310 |
解：
x0=100:100:500;
y0=100:100:400;
z0=[636 697 624 478 450;698 712 630 478 420;680 674 598 412 400;662 626 552 334 310];
x=100:1:500;
y=100:1:400;
%采用interp2函数进行二维插值时，插值点中x，y应是方向不同的向量，故转置
y=y';
z2=interp2(x0,y0,z0,x,y,'spline');
mesh(x,y,z2)
title('二维三次样条插值')
%求最高点地址
[i j]=find(z2==max(max(z2)));
%求最高点坐标
res_x=x(i)
res_y=y(j)
zmax=z2(i,j)图像如图
结果：
res_x =
175
res_y =
170
zmax =
719.77162)插值节点为散点
ZI = griddata(x,y,z,XI,YI)其中 x、y、z均为 n 维向量。向量XI、 YI 是给定的网格点的横坐标和纵坐标，返回值 ZI 为网格（XI，YI）处的函数值。XI 与 YI 应是方向不同的向量，即一个是行向量，另一个是列向量.
例： 在某海域测得一些点（x,y）处的水深 z 由下表给出，在矩形区域内画出海底曲面的图形。
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 |
解:
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];
%求x和y的最大值和最小值，便于定xi和yi的范围
xmm=minmax(x)
ymm=minmax(y)
xi=xmm(1):xmm(2);
yi=ymm(1):ymm(2);
zi1=griddata(x,y,z,xi,yi','cubic');%立方插值
zi2=griddata(x,y,z,xi,yi','nearest');%最近点插值
zi=zi1;
zi(isnan(zi1))=zi2(isnan(zi1));%将立方插值中的不确定值换成最近点插值的结果
%画图
mesh(xi,yi,zi)图像：
二、拟合
拟合：已知平面上的n个点互不相同，寻求一个函数f(x)在某种准则下，使f(x)所有数据点最为接近。（不要求过已知数据点）
1．多项式拟合
a=polyfit(x,y,m)其中输入参数x,y为要拟合的数据，m 为拟合多项式的次数，输出参数 a 为拟合多项式的系数。
注：当使用polyfit进行拟合时，多项式的阶次最大不超过length(x)-1
求多项式在 x 处的值 y 可用：
y=polyval(a,x)
例： 用最小二乘法求一个形如的公式，使它与下表所示的数据拟合。
x | 19 | 25 | 31 | 38 | 44 |
y | 19.0 | 32.3 | 49.0 | 73.3 | 97.8 |
解：
x=[19 25 31 38 44];
y=[19.0 32.3 49.0 73.3 97.8];
a=polyfit(x,y,2)
%设步长为0.1,进行画图
x0=19:0.1:44;
y0=polyval(a,x0);
plot(x,y,'*',x0,y0,'r')图像如图：
结果为：
2．非线性拟合： 已知输入向量xdata和输出向量ydata，并且知道输入与输出的函数关系为ydata=F(x, xdata)，但不知道系数向量x。今进行曲线拟合，求x使得输出的如下最小二乘表达式成立：
在MATLAB中，可用函数lsqcurvefit和lsqnonlin 解决此类问题
lsqcurvefit 函数：
x= lsqcurvefit (‘fun’,x0,xdata,ydata,options)其中fun是定义函数F(x,data)的M文件，x0为初始解向量，xdata,ydata为满足关系ydata=F（x,xdata）的数据
lsqnonlin 函数：
x= lsqnonlin (‘fun’,x0,options)参数说明同上
例：用下面表中的数据拟合函数中的参数a,b,k。
t | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 |
c | 4.54 | 4.99 | 5.35 | 5.65 | 5.90 | 6.10 | 6.26 | 6.39 | 6.50 | 6.59 |
解：首先编写 M 文件myfun.m 定义函数F(x,xdata):
function f=myfun(x,tdata);
% x(1)=a，x(2)=b，x(3)=k
f=x(1)+x(2)*exp(-0.02*x(3)*tdata);然后调用函数 lsqcurvefit
=100:100:1000;
c=[4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59];
x0=[0.01 0.02 0.01];%初始解向量随便取的
x=lsqcurvefit(@myfun,x0,t,c)
%作图查看拟合效果
scatter(t,c)
hold on
c1=myfun(x,t);
plot(t,c1)
hold off图像：
参考文献：
司守奎，孙兆亮，数学建模算法与应用（第二版），北京.国防工业出版社，2019
刘浩，韩晶，MATLAB R2018a完全自学一本通，中国工信出版集团，2019
