Matlab数学建模学习笔记——插值与拟合
⽬录
插值与拟合
插值和拟合的区别
图⽚取⾃知乎⽤户yang元祐的回答
插值：函数⼀定经过原始数据点。
假设f(x)在某区间[a,b]上⼀系列点上的值
y_i=f(x_i),i=0,1,\dots,n。
插值就是⽤较简单、满⾜⼀定条件的函数\varphi(x)去代替f(x)。插值函数满⾜条件
\varphi(x_i)=y_i,i=0,1,\dots,n
拟合：⽤⼀个函数去近似原函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最⼩。
插值⽅法
分段线段插值
分线段插值就是将每两个相邻的节点⽤直线连起来，如此形成的⼀条折线就是就是分段线性插值函数，记作I_n(x)，它满⾜I_n(x_i)=y_i，且I_n(x)在每个⼩区
间[x_i,x_{i+1}]上是线性函数(i=0,1\dots,n-1)。
I_n(x)可以表⽰为I_n(x)=\sum_{i=0}^n y_il_i(x)，其中
l_i(x)= \begin{cases} \frac{x-x_{i-1}}{x_i-x_{i-1}},&x\in [x_{i-1},x_i],i \neq 0,\\ \frac{x-x_{i+1}}{x_i-x_{i+1}},&x\in [x_i,x_{i+1}],i \neq n,\\ 0,&其他 \end{cases}
I_n(x)有良好的收敛性，即对x\in [a,b]，有
\lim _{n \rightarrow \infin}I_n(x)=f(x)
⽤I_n(x)计算x点的插值的时候，只⽤到x左右的两个点，计算量与节点个数n⽆关。但是n越⼤，分段越多，插值误差越⼩。
拉格朗⽇插值多项式
朗格朗⽇(Lagrange)插值的基函数为
\begin{aligned} l_i(x)&=\frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}\\ &= \prod_{j=0\\j\neq i}^{n} \frac{x-x_j}{x_i -x_j},i=0,1,\cdots,n。 \end{aligned}
l_i(x)是xn次多项式，满⾜
l_i(x_j)= \begin{cases} 0,&j\neq i,\\ 1,& j = i。 \end{cases}
拉格朗⽇插值函数函数
L_n(x)=\sum_{i=0}^{n}y_i l_i(x)=\sum_{i=0}^{n} y_i(\prod_{j=0\\j\neq i}^n \frac{x-x_j}{x_i -x_j})
样条插值
早期⼯程师制图时，把富有弹性的细长⽊条（所谓样条）⽤压铁固定在样点上，在其他地⽅让它⾃由弯曲，然后沿⽊条画下曲线。成为样条曲线。绘图员利⽤它把⼀些已知点链接成⼀条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。三次样条插值就是由此抽象出来的。
数学上将具有⼀定光滑性的分段的分段多项式称为样条函数。具体地说，给顶区间[a,b]的⼀个划分。
\Delta:a=x_0 < x_1 < \dots < x_{n-1} < x_n = b。
在每个⼩区间[x_i,x_{i=1}](i=0,1,\dots,n-1)上是S(x)是m次多项式。
S(x)在[a,b]上具有m-1阶连续函数。
则称S(x)为关于划分\Delta的m次样条函数，其图形为m次样条函数。
三次样条插值
已知函数y=f(x)在区间[a,b]上的n+1个节点
\Delta:a=x_0 < x_1 < \dots < x_{n-1} < x_n = b。
的值y_i=f(x_i)(i0,1,\dots,n)，求插值函数S(x)，使得
S(x_i)=y_i(i=0,1,\dots,n)
在每个⼩区间[x_i,x_{i+1}](i=0,1,\dots,n-1)上S(x)是三次多项式，记为S_i(x)
S_i(x)在[a,b]上⼆阶连续可微。
由条件2，我们记
S(x)=\left \{ S_i(x),x\in[x_i,x_{i+1}],i=0,1,\dots,n-1 \right \}\\ S_i(x)=a_i x^3+b_i x^2+c_i x + d_i，
a_i,b_i,c_i,d_i为待定系数，共4n个
由条件3中得⼆阶连续可微，有
\begin{cases} S_i(x_{i+1})=S_{i+1}(x_{i+1}),\\ S_i^{'}(x_{i+1})=S_{i+1}^{'}(x_{i+1}),i=0,1,\dots,n-2。\\ S_i^{''}(x_{i+1})=S_{i+1}^{''}(x_{i+1}),\\ \end{cases}
由上⾯的式⼦共确定4n-2⽅程，为确定S(x)的4n个参数，常⽤的确定三次样条函数边界条件有3种类型：
S'(a)=y_0',S(b)'=y_n'，由这种边界条件建⽴的样条插值函数称为f(x)的完备三次样条插值函数。
特别的，y_0'=y_n'=0时，样条曲线呈⽔平状态。
如果f'(x)不知道，我们可以使S'(x)与f'(x)在端点处近似相等。这时以x_0,x_1,x_2,x_3为节点作⼀个三次Newton插值多项式N_a(x)。同理，以x_n,x_{n-1},x_{n-2},x_{n-3}为节点作⼀个三次Newton插值多项式N_b(x)，要求
S'(a)=N'_a(a),S'(b)=N'_b(b)。
Processing math: 0%
由这种边界条件建⽴的三次样条称为f(x)的Lagrange三次样条插值函数。
S''(a)=y''_0,S''(b)=y''_n。特别地，y''_0=y''_n=0时，称为⾃然边界条件
S'(a+0)=S'(b-0),S''(a+0)=S''(b-0)。此条件称为周期条件。
Matlab插值⼯具箱
⼀维插值函数
interp1函数
y = interp1(x0,y0,x,'method')
% method 为插值⽅法，默认为线性插值，其值可为
% 'nearest' 最近项插值
% 'linear' 线性插值
% 'spline' ⽴⽅样条插值
% 'cubic' ⽴⽅插值
所有的插值⽅法要求x0是单调的。
当x0为等距时可以使⽤快速插值法，使⽤快速插值法的格式为*nearest,*linear,*spline,*cubic
以下为matlab的官⽅说明
vq = interp1(x,v,xq)
vq = interp1(x,v,xq,method)
vq = interp1(x,v,xq,method,extrapolation)
vq = interp1(v,xq)
vq = interp1(v,xq,method)
vq = interp1(v,xq,method,extrapolation)
pp = interp1(x,v,method,'pp')
说明
vq = interp1(x,v,xq)使⽤线性插值返回⼀维函数在特定查询点的插⼊值。向量 x 包含样本点，v 包含对应值 v(x)。向量 xq 包含查询点的坐标。
如果您有多个在同⼀点坐标采样的数据集，则可以将 v 以数组的形式进⾏传递。数组 v 的每⼀列都包含⼀组不同的⼀维样本值。
vq = interp1(x,v,xq,method) 指定备选插值⽅法：'linear'、'nearest'、'next'、'previous'、'pchip'、'cubic'、'v5cubic'、'makima' 'spline'。默认⽅法为 'linear'。vq = interp1(x,v,xq,method,extrapolation) ⽤于指定外插策略，来计算落在 x 域范围外的点。如果希望使⽤ method 算法进⾏外插，可将 extrapolation 设置为'extrap'。您也可以指定⼀个标量值，这种情况下，interp1 将为所有落在 x 域范围外的点返回该标量值。
vq = interp1(v,xq) 返回插⼊的值，并假定⼀个样本点坐标默认集。默认点是从 1 到 n 的数字序列，其中 n 取决于 v 的形状：
当 v 是向量时，默认点是 1:length(v)。
当 v 是数组时，默认点是 1:size(v,1)。
如果您不在意点之间的绝对距离，则可使⽤此语法。
vq = interp1(v,xq,method) 指定备选插值⽅法中的任意⼀种，并使⽤默认样本点。
vq = interp1(v,xq,method,extrapolation) 指定外插策略，并使⽤默认样本点。
pp = interp1(x,v,method,'pp') 使⽤ method 算法返回分段多项式形式的 v(x)。
三次样条插值
Matlab种数据点称为断点。如果三次样条插值没有边界条件，最常⽤的⽅法，就是采⽤⾮扭结(not - a -kont)条件。这个条件强迫第1个和第2个三次多项式的三阶导数相等。对最后⼀个和倒数第2个多项式也做相同的处理。
% matlab中三次样条插值有以下函数
y = interp1(x0,y0,x,'spline');
y = spline(x0,y0,x);
pp = csape(x0,y0,conds);
pp = csape(x0,y0,conds,valconds);y=fnval(pp,x);
% x0, y0是已知数据点；x是插值点，y是插值点的函数值
对于三次样条插值，推荐使⽤函数csape，csape的返回值是pp形式，要获得插值点的函数值，必须调⽤函数fnval，即为pp =
csape(x0,y0,conds,valconds);y=fnval(pp,x);
pp = csape(x0, y0);% 默认边界条件，Lagrange边界条件
pp = csape(x0, y0, conds, valconds);
% valconds 设置边界的⼆阶导数值为[0,0]
% conds指定插值的边界条件，其值可为
% 'complete' 边界我为⼀阶导数，⼀阶导数的值在valconds参数中给出，若忽略valconds参数，按默认情况处理
% 'not - a - knot' ⾮扭结条件
% 'periodic' 周期条件
% 'second'  边界为⼆阶导数，⼆阶导数的值在valconds参数中给出，若忽略valconds参数，按默认情况处理
对于特殊条件，可以通过conds的⼀个1 \times 2矩阵来表⽰，conds的取值为0,1,2
例如，conds=[2,1]的意思为，左边界是⼆阶导数，右边界是⼀阶导数。对应的值由valconds给出。
例题1
如下
t0.150.160.170.18
v(t)3.5  1.5  2.5  2.8
⽤三次样⽅插值求位移S=\int_{0.15}^{0.18}v(t)dt
clc;
clear;
x0=[0.15,0.16,0.17,0.18];
y0=[3.5,1.5,2.5,2.8];
pp=csape(x0,y0); % 默认的边界条件，Lagrange边界条件
format long g
xinshu = pp.coefs; % 显⽰每个区间上三次多项式的系数
s = quadl(@(t)ppval(pp,t),0.15,0.18); % 求积分
format % 恢复短⼩数的显⽰格式
⼆维插值
若节点是⼆维的，插值函数就是⼆元函数，即曲⾯。
Matlab中计算⼆维插值的命令，如：
z = interp2(x0,y0,z0,x,y,'method')
如果是三次样条插值，可以使⽤命令
pp = csape({x0,y0},z0,conds,valconds),z=fnval(pp,{x,y})
例题
⾼程数据点
y \ x100200300400500
100636697624478450
200698712630478420
300680674598412400
400662626552334310
Q：出最⾼点和该点的⾼程。
clc;
clear;
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];
pp = csape({x,y},z');
xi = 100:10:500;
yi = 100:10:400;
cz = fnval(pp,{xi,yi});matlab拟合数据
[i,j]= find(cz==max(max(cz)));
% 要⽤两层max，因为max(cz)为y=180时，和x=100:10:500的⼀系列值，max(max(cz))才是z的最⼤值。 x = xi(i);
y = yi(j);
zmax = cz(i,j);
>> [x,y]
ans =
170  180
>> zmax
zmax =
720.6252
例题2
海底⽔深数据
x129140103.588185.5195105157.5107.57781162162117.5
y7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5
z48686889988949
Q：绘制海底曲⾯的图形
clc;
clear;
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];
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));% 把⽴⽅插值中不确定值换成最近点插值的结果
subplot(1,2,1),plot(x,y,'*');
subplot(1,2,2),mesh(xi,yi,zi);% 绘制三维图形
注：Matlab插值时外插值是不确定的，这⾥使⽤了混合插值，把不确定的插值换成了最近点插值的结果。
曲线拟合的线性最⼩⼆乘法
线性最⼩⼆乘法
公式推导
f(x)=a_1r_1(x)+a_2r_2(x)+\dots+a_mr_m(x)
r_k(x)为事先选好的x⼀组线性⽆关的函数;a_k为待定系数(k=1,2,\dots,m;m<n)。
定义：最⼩⼆乘法就是y_i(k=1,2,\dots,n)与f(x_i)的距离\delta_i的平⽅和最⼩，因此称为最⼩⼆乘法
J(a_1,\dots,a_m)=\sum_{i=1}^n\delta_i^2=\sum_{i=1}^{n}[f(x_i)-y_i]^2
利⽤取得极值的必要条件\frac{\partial J}{\partial a_j}=0，得到关于a_1,\dots,a_m的线性⽅程组，即分别对每⼀个a求偏导。
\sum_{i=1}^n r_j(x_i)\left[ \sum_{k=1}^{m}a_kr_k(x_i)-y_i \right]=0,j=1,\dots,m,
即，
\sum_{i=1}^n a_k\left[ \sum_{k=1}^{m}r_j(x_i)r_k(x_i)\right]= \sum_{k=1}^{m}r_j(x_i)y_i,j=1,\dots,m。
记
R= \left[ \begin{matrix} r_1(x1) & \dots & r_m(x_1)\\ \vdots & \vdots & \vdots\\ r_1(x_n) & \cdots & r_m(x_n)\\ \end{matrix} \right]_{n\times m}\\ A= [a_1,\cdots,a_m]^T,Y=[y_1,\cdots,y_n]^T,
⽅程组可以表⽰为
R^T RA=R^TY。
当\left \{ r_1(x),\cdots,r_m(x) \right \}线性⽆关时，R列满秩，R^TR可逆，于是
A=(R^TR)^{-1}R^TY
函数r_k(x)的选取
常⽤的曲线有
直线y=a_ix+a_2
多项式y=a_1x^m+\cdots+a_mx+a_{m+1}（⼀般m=2,3，不宜太⾼）
双曲线(⼀⽀)y=\frac{a_1}{x}+a_2
指数曲线y=a_1e^{a_2x},
对于指数曲线，拟合前需作变量替换，化为对a1,a2的线性函数
选取时，可在直观判断的基础上，选⼏种曲线分别拟合，然后⽐较，选择最⼩⼆乘指标J最⼩的⼀个。
最⼩⼆乘法的Matlab实现
解⽅程组法
J(a_1,\dots,a_m)=\sum_{i=1}^n\delta_i^2=\sum_{i=1}^{n}[f(x_i)-y_i]^2
记为
J(a_1,\dots,a_m)=\Vert RA-Y \Vert_2^2
Matlab中线性最⼩⼆乘的标准型为
\min_A \Vert RA-Y \Vert_2^2
命令为
A = R\Y
例题5.5
Q:⽤最⼩⼆乘法求⼀个形如y=a+bx^2的经验公式，使其与下列数据表拟合
x1925313844
y19.032.349.073.397.8
clc;
clear;
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')
多项式拟合法
如果取\left \{ r_1(x),\cdots,r_{m+1} \right \}=\left \{ 1,x,\cdots,x^m \right \}，即⽤m次多项式来拟合给定数据。
Matlab命令
a = polyfit(x0,y0,m)
其中，x0,y0为要拟合的数据；m为对项式的次数。输出参数a为拟合多项式y=a(1)x^m+\cdots+a(m)x+a(m+1)的系数向量a=[a(1,),\cdots,a(m),a(m+1)]求多项式在x处的值y可⽤以下命令
y = polyval(a,x)
我们⽤多项式拟合来拟合上⾯的例题
clc;
clear;
x = [19,25,31,38,44]';
y = [19.0,32.3,49.0,73.3,97.8]';
a = polyfit(x,y,2);
xi = 19:0.1:44;
yi = polyval(a,xi);
plot(x,y,'o',xi,yi,'r')
如果我们⽐较⼀下两者的区别
clc;
clear;
x = [19,25,31,38,44]';
y = [19.0,32.3,49.0,73.3,97.8]';

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。