数学建模——插值拟合

最新推荐文章于 2024-08-26 23:12:24 发布

放牛儿

最新推荐文章于 2024-08-26 23:12:24 发布

阅读量4.2k

点赞数 1

分类专栏：数学建模文章标签：矩阵线性代数算法

本文链接：https://blog.csdn.net/m0_51311105/article/details/125836351

版权

数学建模专栏收录该内容

10 篇文章 21 订阅

订阅专栏

插值：每一个点一定在曲线上；拟合：点不一定在曲线上。

如果构造n次插值多项式，则需要n+1个约束方程。插值多项式P(x)与被插函数f(x)之间的差称为截断误差，用R(x)表示。

此时a0到an是未知数，未知数的系数是一个范德蒙行列式，在matlab中，直接用y除以系数矩阵即可求出未知数。

例子：

对下面的观测点进行插值，并估计x=1.5处的函数值

clc,clear
x0 = [1:6]';
y0 = [16 18 21 17 15 12]';
A  = vander(x0);%生成范德蒙行列式
p  = A\y0
x  = [1.5];
y  = polyval(p,x)%用于多项式计算，可查看帮助文档看具体用法

也可以使用matlab的APP中Curve Fitting生成拟合曲线，选择五次线性拟合就会和插值的结果一样，在文件中，可以生成拟合后的代码。

插值

一维插值

vq = interp1(x,v,xq,method,extrapolation)

%向量 x 包含样本点，v 包含对应值 v(x)。向量 xq 包含查询点的坐标。method指定插值方法，如果希望使用 method 算法进行外插，可将 extrapolation 设置为 'extrap'。

F = griddedInterpolant(X,V,Method,ExtrapolationMethod)

%求插值后某点的值时，直接y = F(x)即可

pp = csape(x0,y0);

%三次样条插值函数，返回值是pp形式。

格式函数功能
pp1=csape(x0,y0) 计算插值函数

pp2=fnder(pp1)

计算pp1对应函数的导数，返回值也是pp格式

pp3=fnint(pp1) 计算pp1对应函数的积分，返回值也是pp格式
y=fnval(pp1,x) 计算pp1对应的函数在x点的取值
pp.coefs 显示每个区间函数的系数

相比于interp1，另外两个函数可以生成匿名函数

例子：

数据点如表：

t	0.15	0.16	0.17	0.18
v(t)	3.5	1.5	2.5	2.8

用三次样条插值求 $s=\int_{0.15}^{0.18}v(t)dt$

clc,clear,format long g
t0 = 0.15:0.01:0.18;
v0 = [3.5 1.5 2.5 2.8];
vt = griddedInterpolant(t0,v0,'spline')%求三次插值样条
s1 = integral(@(t)vt(t),0.15,0.18)%求积分

pp=csape(t0,v0);
xishu=pp.coefs
pp3=fnint(pp)
y=fnval(pp3,0.18)-fnval(pp3,0.15)
%s2=integral(@(t)fnval(pp,t),0.15,0.18)
format%恢复短小数显示格式

二维插值

Vq = interp2(X,Y,V,Xq,Yq,method,extrapval)

%X 和 Y 包含样本点的坐标，X是行向量，Y是列向量。V 包含各样本点处的对应函数值。Xq 和 Yq 包含查询点的坐标。

pp = csape({x0,y0},z0) z = fnval(pp,{x,y})

%x0,y0分别为m维和n维向量，z0为m*n维矩阵。x,y为想插值的点

F = griddedInterpolant(X1,X2,V,Method,ExtrapolationMethod)

%求插值后某点的值时，直接y = F(x1,x2)即可

例如：

x = -5:0.8:5;
y = x';
z = sin(x.^2 + y.^2) ./ (x.^2 + y.^2);
F = griddedInterpolant({x,y},z);
xq = -5:0.1:5;
yq = xq';
vq = F({xq,yq});
surf(xq,yq,vq)

散乱数据插值

vq = griddata(x,y,v,xq,yq) 使 v = f(x,y) 形式的曲面与向量 (x,y,v) 中的散点数据拟合。griddata 函数在 (xq,yq) 指定的查询点对曲面进行插值并返回插入的值 vq。曲面始终穿过 x 和 y 定义的数据点，x,y应该是一个行向量，一个列向量。

F = scatteredInterpolant(x,y,v) 创建一个拟合 v = F(x,y) 形式的曲面的插值。向量 x 和 y 指定样本点的 (x,y) 坐标。v 是一个包含与点 (x,y) 关联的样本值的向量。也相当于一个匿名函数

拟合

线性最小二乘法实现

线性方程组：RA=Y。求取参数向量A方法：A=R\Y

例子：

根据实验数据，建立经验公式y=at+b

t	0	1	2	3	4	5	6	7
y	27.0	26.8	26.5	26.3	26.1	25.7	25.3	24.8

clc,clear
t = [0:7]';
y = [27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8]';
a = [t,ones(8,1)];
cs= a\y

用约束线性最小二乘解

x = lsqlin(C,d,A,b,Aeq,beq,lb,ub) ，C*x = d，线性等式约束 Aeq*x = beq 和边界 lb ≤ x ≤ ub。

例子：

给定拟合函数 $y=ae^{x}+blnx$ ，满足a，b $\geqslant$ 0， $a+b\leqslant 1$

x	3	5	6	7	4	8	5	9
y	4	9	5	3	8	5	8	5

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)

多项式拟合

p = polyfit(x,y,n);%拟合n次多项式，返回值p是多项式对应的系数，次序是从高次幂到低次幂。

y = ployval(p,x);%计算多项式p在x处的函数值。

fittype和fit函数

函数fit需要和函数fittype配合使用，fittype用于定义拟合的函数类，fit用于函数拟合。fit可以拟合一元和二元函数，线性和非线性都可以。

aFittype = fittype(libraryModelName)%利用库模型定义函数类

aFittype = fittype(expression,Name,Value)%利用字符串定义函数类

aFittype = fittype(linearModelTerms,Name,Value)%利用基函数的线性组合定义函数类

aFittype = fittype(anonymousFunction,Name,Value)%利用匿名函数定义函数类

fitobject = fit(x,y,fitType)%x，y分别为自变量和因变量的观测值列向量，返回值是拟合函数的信息。

fitobject = fit([x,y],z,fitType)%[x,y]为自变量的观测值的两列矩阵，z为因变量的观测值列向量，可以拟合二元函数。

[fitobject,gof] = fit(x,y,fitType)%gof可以给出模型的一些检验统计量。

例子：

拟合函数

clc,clear
x1 = [6 2 6 7 4 2 5 9]';
x2 = [4 9 5 3 8 5 8 2]';
z  = [14.2077 39.3622 17.8077 11.8310 32.8618 16.9622 33.0941 11.1737]';
g  = fittype('a*exp(b*x)+c*y^2','dependent','z','independent',{'x','y'});
[f,st] = fit([x1,x2],z,g,'StartPoint',rand(1,3))
plot(f)

拟合分段函数：

g = @(a,b,c,d,k,x)(a+b*x).*(x<k)+(c+d*x).*(x>=k);%注意点乘
g = fittype(g)
...

匿名函数定义中，必须把自变量放在匿名函数输入参数的最后一个。

lsqcurvefit函数

该函数可以拟合多元函数。

x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)

xdata是x观测值，ydata是y观测值， $\theta$ 是参数向量

拟合工具

最好用的还是APP中的拟合工具，可以生成代码，也可以导入到工作空间。如果导入到工作空间，就相当于一个匿名函数，可以直接使用。

其他

①trapz函数可以通过梯形法计算近似积分。

②gradient是求数值梯度的函数。[Fx,Fy]=gradient(F，hx，hy)，其中Fx为其水平方向上的梯度，Fy为其垂直方向上的梯度，Fx的第一列元素为原矩阵第二列与第一列元素之差除以hx，Fx的第二列元素为原矩阵第三列与第一列元素之差除以(2*hx)，以此类推：Fx(i,j)=(F(i,j+1)-F(i,j-1))/(2*hx)。最后一列则为最后两列之差除以hx。同理，可以得到Fy。

③contour函数可以画出矩阵的等高线图

④meshgrid(x,y)基于向量x和y中包含的坐标返回二维网格坐标

⑤quiver(x,y,u,v)该函数使用箭头来直观的显示矢量场。该调用格式表示通过在(x, y)指定的位置绘制小箭头来表示以该点为起点的向量(u,v)。

⑥clabel(h) 为当前等高线图h添加标签

⑦norm(v) 返回向量 v 的欧几里德范数。此范数也称为 2-范数、向量模或欧几里德长度。

⑧fimplicit绘制隐函数图像
f = @(x,y) x.^2 + y.^2 - 3;
fimplicit(f,[-3 0 -2 2])
⑨title、xlabel、ylabel、zlabel、textbox和legend等的Interpreter属性有三个属性：latex 、tex、none

⑩求两个方程的交点（f本身就是一个匿名函数，求f和2x-1的交点，rand赋初值）

格式	函数功能
pp1=csape(x0,y0)	计算插值函数
pp2=fnder(pp1)	计算pp1对应函数的导数，返回值也是pp格式
pp3=fnint(pp1)	计算pp1对应函数的积分，返回值也是pp格式
y=fnval(pp1,x)	计算pp1对应的函数在x点的取值
pp.coefs	显示每个区间函数的系数