《数学建模算法与应用》——插值与拟合
文章目录
- 插值与拟合
- 插值和拟合的区别
- 插值方法
- 分段线段插值
- 拉格朗日插值多项式
- 样条插值
- 三次样条插值
- Matlab插值工具箱
- 一维插值函数
- interp1函数
- 三次样条插值
- 例题1
- 二维插值
- 例题
- 例题2
- 曲线拟合的线性最小二乘法
- 线性最小二乘法
- 公式推导
- 函数rk(x)r_k(x)rk(x)的选取
- 最小二乘法的Matlab实现
- 解方程组法
- 例题5.5
- 多项式拟合法
- 最小二乘优化
- lsqlin函数
- lsqcurvefit函数
- 例题
- lsqnonlin函数
- lsqnonneg函数
- Matlab的曲线拟合用户图形界面解法
- 曲线拟合与函数逼近
- 曲线拟合
- 函数逼近
- 例题
插值与拟合
插值和拟合的区别
图片取自知乎用户yang元祐的回答
插值:函数一定经过原始数据点。
假设f(x)在某区间[a,b]上一系列点上的值
yi=f(xi),i=0,1,…,n。y_i=f(x_i),i=0,1,\dots,n。 yi=f(xi),i=0,1,…,n。
插值就是用较简单、满足一定条件的函数φ(x)\varphi(x)φ(x)去代替f(x)f(x)f(x)。插值函数满足条件
φ(xi)=yi,i=0,1,…,n\varphi(x_i)=y_i,i=0,1,\dots,n φ(xi)=yi,i=0,1,…,n
拟合:用一个函数去近似原函数,不要求过已知数据点,只要求在某种意义下它在这些点上的总偏差最小。
插值方法
分段线段插值
分线段插值就是将每两个相邻的节点用直线连起来,如此形成的一条折线就是就是分段线性插值函数,记作In(x)I_n(x)In(x),它满足In(xi)=yiI_n(x_i)=y_iIn(xi)=yi,且In(x)I_n(x)In(x)在每个小区间[xi,xi+1][x_i,x_{i+1}][xi,xi+1]上是线性函数(i=0,1…,n−1)(i=0,1\dots,n-1)(i=0,1…,n−1)。
In(x)I_n(x)In(x)可以表示为In(x)=∑i=0nyili(x)I_n(x)=\sum_{i=0}^n y_il_i(x)In(x)=∑i=0nyili(x),其中
li(x)={x−xi−1xi−xi−1,x∈[xi−1,xi],i≠0,x−xi+1xi−xi+1,x∈[xi,xi+1],i≠n,0,其他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} li(x)=⎩⎪⎨⎪⎧xi−xi−1x−xi−1,xi−xi+1x−xi+1,0,x∈[xi−1,xi],i=0,x∈[xi,xi+1],i=n,其他
In(x)I_n(x)In(x)有良好的收敛性,即对x∈[a,b]x\in [a,b]x∈[a,b],有
limn→∞In(x)=f(x)\lim _{n \rightarrow \infin}I_n(x)=f(x) n→∞limIn(x)=f(x)
用In(x)I_n(x)In(x)计算x点的插值的时候,只用到x左右的两个点,计算量与节点个数n无关。但是n越大,分段越多,插值误差越小。
拉格朗日插值多项式
朗格朗日(Lagrange)插值的基函数为
li(x)l_i(x)li(x)是xn次多项式,满足
li(xj)={0,j≠i,1,j=i。l_i(x_j)= \begin{cases} 0,&j\neq i,\\ 1,& j = i。 \end{cases} li(xj)={0,1,j=i,j=i。
拉格朗日插值函数函数
Ln(x)=∑i=0nyili(x)=∑i=0nyi(∏j=0j≠inx−xjxi−xj)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}) Ln(x)=i=0∑nyili(x)=i=0∑nyi(j=0j=i∏nxi−xjx−xj)
样条插值
早期工程师制图时,把富有弹性的细长木条(所谓样条)用压铁固定在样点上,在其他地方让它自由弯曲,然后沿木条画下曲线。成为样条曲线。绘图员利用它把一些已知点链接成一条光滑曲线(称为样条曲线),并使连接点处有连续的曲率。三次样条插值就是由此抽象出来的。
数学上将具有一定光滑性的分段的分段多项式称为样条函数。具体地说,给顶区间[a,b]的一个划分。
Δ:a=x0<x1<⋯<xn−1<xn=b。\Delta:a=x_0 < x_1 < \dots < x_{n-1} < x_n = b。 Δ:a=x0<x1<⋯<xn−1<xn=b。
- 在每个小区间[xi,xi=1](i=0,1,…,n−1)[x_i,x_{i=1}](i=0,1,\dots,n-1)[xi,xi=1](i=0,1,…,n−1)上是S(x)是m次多项式。
- S(x)在[a,b]上具有m-1阶连续函数。
则称S(x)为关于划分Δ\DeltaΔ的m次样条函数,其图形为m次样条函数。
三次样条插值
已知函数y=f(x)y=f(x)y=f(x)在区间[a,b]上的n+1个节点
Δ:a=x0<x1<⋯<xn−1<xn=b。\Delta:a=x_0 < x_1 < \dots < x_{n-1} < x_n = b。 Δ:a=x0<x1<⋯<xn−1<xn=b。
的值yi=f(xi)(i0,1,…,n)y_i=f(x_i)(i0,1,\dots,n)yi=f(xi)(i0,1,…,n),求插值函数S(x),使得
- S(xi)=yi(i=0,1,…,n)S(x_i)=y_i(i=0,1,\dots,n)S(xi)=yi(i=0,1,…,n)
- 在每个小区间[xi,xi+1](i=0,1,…,n−1)[x_i,x_{i+1}](i=0,1,\dots,n-1)[xi,xi+1](i=0,1,…,n−1)上S(x)是三次多项式,记为Si(x)S_i(x)Si(x)
- Si(x)S_i(x)Si(x)在[a,b]上二阶连续可微。
由条件2,我们记
S(x)={Si(x),x∈[xi,xi+1],i=0,1,…,n−1}Si(x)=aix3+bix2+cix+di,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, S(x)={Si(x),x∈[xi,xi+1],i=0,1,…,n−1}Si(x)=aix3+bix2+cix+di,
ai,bi,ci,dia_i,b_i,c_i,d_iai,bi,ci,di为待定系数,共4n个
由条件3中得二阶连续可微,有
{Si(xi+1)=Si+1(xi+1),Si′(xi+1)=Si+1′(xi+1),i=0,1,…,n−2。Si′′(xi+1)=Si+1′′(xi+1),\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} ⎩⎪⎨⎪⎧Si(xi+1)=Si+1(xi+1),Si′(xi+1)=Si+1′(xi+1),i=0,1,…,n−2。Si′′(xi+1)=Si+1′′(xi+1),
由上面的式子共确定4n-2方程,为确定S(x)的4n个参数,常用的确定三次样条函数边界条件有3种类型:
S′(a)=y0′,S(b)′=yn′S'(a)=y_0',S(b)'=y_n'S′(a)=y0′,S(b)′=yn′,由这种边界条件建立的样条插值函数称为f(x)的完备三次样条插值函数。
特别的,y0′=yn′=0y_0'=y_n'=0y0′=yn′=0时,样条曲线呈水平状态。如果f′(x)f'(x)f′(x)不知道,我们可以使S′(x)S'(x)S′(x)与f′(x)f'(x)f′(x)在端点处近似相等。这时以x0,x1,x2,x3x_0,x_1,x_2,x_3x0,x1,x2,x3为节点作一个三次Newton插值多项式Na(x)N_a(x)Na(x)。同理,以xn,xn−1,xn−2,xn−3x_n,x_{n-1},x_{n-2},x_{n-3}xn,xn−1,xn−2,xn−3为节点作一个三次Newton插值多项式Nb(x)N_b(x)Nb(x),要求
S′(a)=Na′(a),S′(b)=Nb′(b)。S'(a)=N'_a(a),S'(b)=N'_b(b)。 S′(a)=Na′(a),S′(b)=Nb′(b)。
由这种边界条件建立的三次样条称为f(x)f(x)f(x)的Lagrange三次样条插值函数。S′′(a)=y0′′,S′′(b)=yn′′S''(a)=y''_0,S''(b)=y''_nS′′(a)=y0′′,S′′(b)=yn′′。特别地,y0′′=yn′′=0y''_0=y''_n=0y0′′=yn′′=0时,称为自然边界条件
S′(a+0)=S′(b−0),S′′(a+0)=S′′(b−0)S'(a+0)=S'(b-0),S''(a+0)=S''(b-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)。
interp1官方文档
三次样条插值
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×21 \times 21×2矩阵来表示,conds的取值为0,1,2
例如,conds=[2,1]的意思为,左边界是二阶导数,右边界是一阶导数。对应的值由valconds给出。
csape官方文档
例题1
如下
t | 0.15 | 0.16 | 0.17 | 0.18 |
---|---|---|---|---|
v(t) | 3.5 | 1.5 | 2.5 | 2.8 |
用三次样方插值求位移S=∫0.150.18v(t)dtS=\int_{0.15}^{0.18}v(t)dtS=∫0.150.18v(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})
interp2官方文档
例题
高程数据点
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 |
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});[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>> zmaxzmax =720.6252
例题2
海底水深数据
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 |
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)=a1r1(x)+a2r2(x)+⋯+amrm(x)f(x)=a_1r_1(x)+a_2r_2(x)+\dots+a_mr_m(x) f(x)=a1r1(x)+a2r2(x)+⋯+amrm(x)
rk(x)r_k(x)rk(x)为事先选好的x一组线性无关的函数;aka_kak为待定系数(k=1,2,…,m;m<n)(k=1,2,\dots,m;m<n)(k=1,2,…,m;m<n)。
**定义:**最小二乘法就是yi(k=1,2,…,n)y_i(k=1,2,\dots,n)yi(k=1,2,…,n)与f(xi)f(x_i)f(xi)的距离δi\delta_iδi的平方和最小,因此称为最小二乘法
J(a1,…,am)=∑i=1nδi2=∑i=1n[f(xi)−yi]2J(a_1,\dots,a_m)=\sum_{i=1}^n\delta_i^2=\sum_{i=1}^{n}[f(x_i)-y_i]^2 J(a1,…,am)=i=1∑nδi2=i=1∑n[f(xi)−yi]2
利用取得极值的必要条件∂J∂aj=0\frac{\partial J}{\partial a_j}=0∂aj∂J=0,得到关于a1,…,ama_1,\dots,a_ma1,…,am的线性方程组,即分别对每一个a求偏导。
∑i=1nrj(xi)[∑k=1makrk(xi)−yi]=0,j=1,…,m,\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, i=1∑nrj(xi)[k=1∑makrk(xi)−yi]=0,j=1,…,m,
即,
∑i=1nak[∑k=1mrj(xi)rk(xi)]=∑k=1mrj(xi)yi,j=1,…,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。 i=1∑nak[k=1∑mrj(xi)rk(xi)]=k=1∑mrj(xi)yi,j=1,…,m。
记
R=[r1(x1)…rm(x1)⋮⋮⋮r1(xn)⋯rm(xn)]n×mA=[a1,⋯,am]T,Y=[y1,⋯,yn]T,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=⎣⎢⎡r1(x1)⋮r1(xn)…⋮⋯rm(x1)⋮rm(xn)⎦⎥⎤n×mA=[a1,⋯,am]T,Y=[y1,⋯,yn]T,
方程组可以表示为
RTRA=RTY。R^T RA=R^TY。 RTRA=RTY。
当{r1(x),⋯,rm(x)}\left \{ r_1(x),\cdots,r_m(x) \right \}{r1(x),⋯,rm(x)}线性无关时,R列满秩,RTRR^TRRTR可逆,于是
A=(RTR)−1RTYA=(R^TR)^{-1}R^TY A=(RTR)−1RTY
函数rk(x)r_k(x)rk(x)的选取
常用的曲线有
- 直线y=aix+a2y=a_ix+a_2y=aix+a2
- 多项式y=a1xm+⋯+amx+am+1y=a_1x^m+\cdots+a_mx+a_{m+1}y=a1xm+⋯+amx+am+1(一般m=2,3,不宜太高)
- 双曲线(一支)y=a1x+a2y=\frac{a_1}{x}+a_2y=xa1+a2
- 指数曲线y=a1ea2xy=a_1e^{a_2x}y=a1ea2x,
- 对于指数曲线,拟合前需作变量替换,化为对a1,a2的线性函数
选取时,可在直观判断的基础上,选几种曲线分别拟合,然后比较,选择最小二乘指标J最小的一个。
最小二乘法的Matlab实现
解方程组法
J(a1,…,am)=∑i=1nδi2=∑i=1n[f(xi)−yi]2J(a_1,\dots,a_m)=\sum_{i=1}^n\delta_i^2=\sum_{i=1}^{n}[f(x_i)-y_i]^2 J(a1,…,am)=i=1∑nδi2=i=1∑n[f(xi)−yi]2
记为
J(a1,…,am)=∥RA−Y∥22J(a_1,\dots,a_m)=\Vert RA-Y \Vert_2^2 J(a1,…,am)=∥RA−Y∥22
Matlab中线性最小二乘的标准型为
minA∥RA−Y∥22\min_A \Vert RA-Y \Vert_2^2 Amin∥RA−Y∥22
命令为
A = R\Y
例题5.5
Q:用最小二乘法求一个形如y=a+bx2y=a+bx^2y=a+bx2的经验公式,使其与下列数据表拟合
x | 19 | 25 | 31 | 38 | 44 |
---|---|---|---|---|---|
y | 19.0 | 32.3 | 49.0 | 73.3 | 97.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')
多项式拟合法
如果取{r1(x),⋯,rm+1}={1,x,⋯,xm}\left \{ r_1(x),\cdots,r_{m+1} \right \}=\left \{ 1,x,\cdots,x^m \right \}{r1(x),⋯,rm+1}={1,x,⋯,xm},即用m次多项式来拟合给定数据。
Matlab命令
a = polyfit(x0,y0,m)
其中,x0,y0为要拟合的数据;m为对项式的次数。输出参数a为拟合多项式y=a(1)xm+⋯+a(m)x+a(m+1)y=a(1)x^m+\cdots+a(m)x+a(m+1)y=a(1)xm+⋯+a(m)x+a(m+1)的系数向量a=[a(1,),⋯,a(m),a(m+1)]a=[a(1,),\cdots,a(m),a(m+1)]a=[a(1,),⋯,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]';
a = polyfit(x,y,2);
xi = 19:0.1:44;
yi = polyval(a,xi);
r = [ones(5,1),x.^2];
ab=r\y;
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x0,y0,xi,yi)
legend('最小二乘','多项式拟合')
我们看到其实两者差别不大的,如果我们看一看系数
a(n)a(n)a(n) | a(1) | a(2) | a(3) |
---|---|---|---|
多项式拟合 y1y_1y1 | 0.0497 | 0.0193 | 0.6881 |
最小二乘法 y2y_2y2 | 0.0500 | 0 | 0.9725 |
y1=0.0497x2+0.0193x+0.6881y2=0.9725x2+0.05y_1=0.0497x^2+0.0193x+0.6881\\ y_2=0.9725x^2+0.05 y1=0.0497x2+0.0193x+0.6881y2=0.9725x2+0.05
最小二乘优化
在无约束优化问题中,有些情形,比如目标函数由若干个函数的平方和构成,这类函数一般可以写成
F(x)=∑i=1mfi2(x),x∈Rn,F(x)=\sum_{i=1}^mf_i^2(x),x\in R^n, F(x)=i=1∑mfi2(x),x∈Rn,
式中,x=[x1,⋯,xn]Tx=[x1,\cdots,x_n]^Tx=[x1,⋯,xn]T,一般假设m≥nm\geq nm≥n。
把极小化这类函数的问题
minF(x)=∑i=1mfi2(x)\min F(x)=\sum_{i=1}^mf_i^2(x) minF(x)=i=1∑mfi2(x)
称为最小二乘优化问题。
在Matlab优化工具箱中,有
lsqlin, lsqcurvefit, lsqnonlin, isqnonneg等函数
lsqlin函数
求解
minx12∥C⋅x−d∥22s.t.{A⋅x≤b,Aeq⋅x=beq,lb≤x≤ub,\min _x \frac{1}{2} \Vert C \cdot x -d \Vert_2^2\\ s.t. \begin{cases} A \cdot x \leq b,\\ Aeq \cdot x = beq,\\ lb \leq x \leq ub, \end{cases} xmin21∥C⋅x−d∥22s.t.⎩⎪⎨⎪⎧A⋅x≤b,Aeq⋅x=beq,lb≤x≤ub,
式中,C, Aeq, A为矩阵;d, b, beq, lb, ub, x为向量
Matlab中的函数为
x = lsqlin(C,d,A,b)
x = lsqlin(C,d,A,b,Aeq,beq,lb,ub)
x = lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0,options)
x = lsqlin(problem)
[x,resnorm,residual,exitflag,output,lambda] = lsqlin(___)
%lsqlin命令求解例5.5
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=lsqlin(r,y);
x0=19:0.1:44;
y0=ab(1)+ab(2)*x0.^2;
plot(x,y,'o',x0,y0,'r')
计算结果是一样的
lsqlin官方文档
lsqcurvefit函数
给定输入输出数列xdata,ydata,求参量x,使得
minx∥F(x,xdata)−ydata∥22=∑i(F(x,xdatai)−ydatai)2\min _x \Vert F(x,xdata)-ydata \Vert_2^2 = \sum_i(F(x,xdata_i)-ydata_i)^2 xmin∥F(x,xdata)−ydata∥22=i∑(F(x,xdatai)−ydatai)2
Matlab中的函数为
x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)
% fun为定义函数F(x,xdata)的M文件
注:非线性拟合时,每一次的运行结果可能都是不同的。
例题
Q:用最小二乘法拟合y=12πσe−(x−μ)22σ2y=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}}y=2πσ1e−2σ2(x−μ)2,其中未知参数为σ,μ\sigma,\muσ,μ
clc;
clear;
x0 = -10:0.01:10;
y0 = normpdf(x0,0,1);
mf=@(cs,xdata)1/sqrt(2*pi)/cs(2)*exp(-(xdata-cs(1)).^2/cs(2)^2/2);
% yc = mf([2,1],1);% 测试匿名函数
cs = lsqcurvefit(mf,rand(2,1),x0,y0);% 拟合参数的初始值时任意取的
% 计算出来的估计值 cs(1)=0,cs(2)=1
lsqcurvefit官方文档
lsqnonlin函数
已知函数向量F(x)=[f1(x),⋯,fk(x)]TF(x)=[f_1(x),\cdots,f_k(x)]^TF(x)=[f1(x),⋯,fk(x)]T,使x使得
minx∥F(x)∥22\min _x \Vert F(x) \Vert_2^2 xmin∥F(x)∥22
Matlab中的函数为
x = lsqnonlin(fun,x0,lb,ub,options)
% fun为定义向量函数F(x)的M文件
lsqnonlin官方文档
lsqnonneg函数
求解非负的x,使得,
minx∥Cx−d∥22\min _x \Vert Cx-d \Vert_2^2 xmin∥Cx−d∥22
Matlab中的函数为
x = lsqnonneq(C,d,options)
lsqnonneq官方文档
Matlab的曲线拟合用户图形界面解法
Matlab工具箱提供了命令cftool
,该命令给出了一维数据拟合的交互式环境。
执行步骤:
- 把数据导入到工作空间
- 运行
cftool
,打开用户图形界面窗口 - 选择适当的模型进行拟合
- 生成一些相关的统计量
曲线拟合与函数逼近
曲线拟合
曲线拟合是已知一组离散数据{(xi,yi),i=1,⋯,n}\left \{ (x_i,y_i),i=1,\cdots,n \right \}{(xi,yi),i=1,⋯,n},选择一个较简单的的函数f(x)(如多项式),在一定的准则(如最小二乘法准则)下,最接近这些数据。
函数逼近
如果已知一个较为复杂的连续函数f(x),x∈[a,b]f(x),x\in [a,b]f(x),x∈[a,b],要求选择一个较简单的函数f(x),在一定的准则下最接近f(x),就是所谓的函数逼近
与最小二乘准则相对应,函数逼近常采用的一种准则是最小平方逼近
J=∫ab[f(x)−y(x)]2dxJ=\int_a^b [f(x)-y(x)]^2dx J=∫ab[f(x)−y(x)]2dx
达到最小。与曲线拟合一样,选一组函数{rk(x),k=1,⋯,m}\left \{ r_k(x),k=1,\cdots,m \right \}{rk(x),k=1,⋯,m}构造函数f(x),即令
f(x)=a1r1(x)+⋯+amrm(x),f(x)=a_1r_1(x)+\cdots+a_mr_m(x), f(x)=a1r1(x)+⋯+amrm(x),
带入上式中,求a1,⋯,ama_1,\cdots,a_ma1,⋯,am使J达到最小。利用极值必要条件可得
这里(g,h)=∫abg(x)h(x)dx(g,h)=\int_a^b g(x)h(x)dx(g,h)=∫abg(x)h(x)dx,当方程组的系数矩阵非奇异时,有唯一解。
最简单的是使用多项式逼近,r1(x)=1,r2(x)=x,r3(x)=x2,⋯r_1(x)=1,r_2(x)=x,r_3(x)=x^2,\cdotsr1(x)=1,r2(x)=x,r3(x)=x2,⋯。并且如果能使∫abri(x)rj(x)dx=0,i≠j\int_a^b r_i(x)r_j(x)dx=0,i \neq j∫abri(x)rj(x)dx=0,i=j,方程组的系数矩阵将是对角阵,计算大大简化,满足这种性质的多项式称为正交多项式。
勒让德(Legendre)多项式是在[-1,1]区间上的正交多项式,它的表示式为
P0(x)=1,Pk(x)=12kk!dkdxk(x2−1)k,k=1,2,⋯P_0(x)=1,P_k(x)=\frac{1}{2^kk!}\frac{d^k}{dx^k}(x^2-1)^k,k=1,2,\cdots P0(x)=1,Pk(x)=2kk!1dxkdk(x2−1)k,k=1,2,⋯
可以证明
∫−11Pi(x)Pj(x)dx={0,i≠j,22i+1,i=j,Pk+1(x)=2k+1k+1xPk(x)−kk+1Pk−1(x),k=1,2,⋯。\int_{-1}^1 P_i(x)P_j(x)dx= \begin{cases} 0,&i \neq j,\\ \frac{2}{2i+1},&i=j, \end{cases} \\ P_{k+1}(x)=\frac{2k+1}{k+1}xP_k(x)-\frac{k}{k+1}P_{k-1}(x),k=1,2,\cdots。 ∫−11Pi(x)Pj(x)dx={0,2i+12,i=j,i=j,Pk+1(x)=k+12k+1xPk(x)−k+1kPk−1(x),k=1,2,⋯。
常用的正交多项式还有第一类切比雪夫(Chebyshev)多项式
Tn(x)=cos(narccosx),x∈[−1,1],n=0,1,2,⋯T_n(x)=\cos(narccosx),x\in [-1,1],n=0,1,2,\cdots Tn(x)=cos(narccosx),x∈[−1,1],n=0,1,2,⋯
和拉盖尔(Laguerre)多项式
Ln(x)=exdndxn(xne−x),x∈[0,+∞),n=0,1,2,⋯L_n(x)=e^x \frac{d^n}{dx^n}(x^ne^{-x}),x\in [0,+\infin),n=0,1,2,\cdots Ln(x)=exdxndn(xne−x),x∈[0,+∞),n=0,1,2,⋯
例题
求f(x)=cosx,x∈[−π2,π2]f(x)=\cos x,x\in[-\frac{\pi}{2},\frac{\pi}{2}]f(x)=cosx,x∈[−2π,2π]在H=Span{1,x2,x4}H=Span \left\{ 1,x^2,x^4 \right\}H=Span{1,x2,x4}中最佳平方逼近多项式。
clc;
clear;
syms x
base=[1,x^2,x^3];
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); % 把符号数据转化为保留6位有效数字的符号数据
《数学建模算法与应用》——插值与拟合相关推荐
- c语言埃尔米特插值思路,【数学建模算法】(26)插值和拟合:埃尔米特(Hermite)插值和样条插值...
1.埃尔米特(Hermite)插值 1.1.Hermite插值多项式 如果对插值函数,不仅要求它在节点处与函数同值,而且要求它与函数有相同的一阶.二阶甚至更高阶的导数值,这就是 Hermite 插值问 ...
- matlab中x从0到5不含0,关于MATLAB的数学建模算法学习笔记
关于MATLAB的数学建模算法学习笔记 目录 线性规划中应用: (3) 非线性规划: (3) 指派问题;投资问题:(0-1问题) (3) 1)应用fmincon命令语句 (3) 2)应用指令函数:bi ...
- 数学建模算法学习笔记 已完结
这是为了准备国赛突击学习的模型算法,我在原有的基础上加上自己的理解虽然不知道对不对,就是为了记录下自己学的模型他究竟是个什么东西,语言通俗,但是极不准确,只适合做一个大概的了解,建议大家详细的还是要看 ...
- 数学建模算法学习笔记
数学建模算法学习笔记 作为建模Man学习数学建模时做的笔记 参考文献: <数学建模姜启源第四版> 网上搜罗来的各种资料,侵删 1.线性预测 levinson durbin算法,自相关什么的 ...
- 清风:数学建模算法、编程和写作培训
清风:数学建模算法.编程和写作培训 一.评价模型 1.1 层次分析法 1.2 代码详解 1.3 模型拓展 1.4 课后作业 二.插值与拟合模型 三.相关性模型 四.回归模型 五.图论模型 六.分类问题 ...
- LL1分析构造法_数学建模算法--最优赋权法(含代码)
数学建模算法--最优赋权法(含代码) 作者:郑铿城 本次介绍数学建模和科研写作的方法--最优赋权法最优赋权法经常用于分析评价类问题,从该算法的名称就可以看到,该算法首先要体现"最优" ...
- python dendrogram_【聚类分析】《数学建模算法与应用》第十章 多元分析 第一节 聚类分析 python实现...
第十章 多元分析 第一节 聚类分析 介绍 这里是司守奎教授的<数学建模算法与应用>全书案例代码python实现,欢迎加入此项目将其案例代码用python实现 GitHub项目地址:Math ...
- 数学建模算法:支持向量机_从零开始的算法:支持向量机
数学建模算法:支持向量机 从零开始的算法 (Algorithms From Scratch) A popular algorithm that is capable of performing lin ...
- 数学建模算法与应用:预测模型(3)案例: SARS 疫情对经济指标影响
目录 问题描述: 一.建模思路 二.对模型进行分析预测 2.1.对模型进行假设 三.建立灰色预测模型GM(1,1) 3.1.模型的求解(i)商品零售额 3.2.用MATLAB程序,实现(i)商品零售额 ...
- 数学建模算法与应用 线性规划(cvxpy包)
数学建模算法与应用 线性规划(使用cvxpy包) 说明 使用python中cvxpy库完成<数学建模算法与应用>中课后习题 因为本人也是初学者,若代码有错误还请各位指出 cvxpy库的使用 ...
最新文章
- 【代码片段收集】Python解析AndroidManifest.xml
- 2021年4月28日 深圳头条后台开发实习面试(hr面)
- http statusCode(状态码)
- 博客创办目的——————欢迎相互学习
- 2019山东省赛B - Flipping Game ZOJ - 4114 题解
- 【视觉项目】【day2】8.21号实验记录(手机固定高度15cm拍摄+直方图均衡化+模板匹配,模板12个,测试28个,效果十分差)
- 计算机考研英语词汇书,求助:有知道电脑背考研英语单词的
- 二三星缩水软件手机版_还在抱怨三星手机不好用?用这些软件立马解决
- LeetCode 31. 下一个排列(线性扫描)
- 24种设计模式--命令模式【Command Pattern】
- adboost,随机森林,gbdt,xgboost,lightgbm区别
- springboot+vue3+微信小程序活动报名系统源码
- html使用css居中
- 【七夕节特刊】开源世界里的爱情保卫战
- Android学习笔记--Notification(通知)
- java实现自行车行程
- golang中的dns问题
- Python判断指定日期是不是法定节假日
- 如何在 iPhone 上设置整点报时提醒?
- android type c 耳机检测,USB Type-C 的新音频标准将帮助 Android 设备去掉 3.5mm 耳机孔...
热门文章
- 33 个 2017 年必须了解的 iOS 开源库
- 在Word文档中插入矢量图片
- 元素周期表,列表格式
- 无法登陆计算机报名系统,广东省计算机等级考试无法登录报名系统网页的看这里...
- rocksdb原理_Rocksdb Compaction原理
- linux上mysql忘记密码,linux下mysql忘记密码解决方案
- [BZOJ2428] [HAOI2006]均分数据 模拟退火
- 一度智信:关于拼多多商家经验分享
- 刘强东对京东零售动刀:提醒打工仔,要立新功不吃老本
- php 微信统一下单接口,微信JSAPI支付,统一下单接口