拟合与插值|线性最小二乘拟合|非线性最小二乘拟合|一维插值|二维插值

ChoSeitaku

于 2024-08-19 22:15:50 发布

阅读量643

点赞数 7

分类专栏：数学建模文章标签：数学建模

本文链接：https://blog.csdn.net/CoderZzz6310/article/details/141336067

版权

数学建模专栏收录该内容

12 篇文章 0 订阅

订阅专栏

挖掘数据背后的规律是数学建模的重要任务，拟合与插值是常用的分析方法

掌握拟合与插值的基本概念和方法
熟悉Matlab相关程序实现
能够从数据中挖掘数学规律

拟合问题的基本提法

拟合问题的概念

已知一组数据(以二维为例)，即平面上n个点 $(x_{i},y_{i}),i =1,2,\dots,n$ ，寻求一个函数 $y = f (x)$ ，使得 $f (x)$ 在某种准则下与所有数据点最为接近
![[Pasted image 20240819085354.png]]

$\sum_{i=1}^{n}(\delta_{i})^{2}\to min$

确定 $f (x)$ 表达式的形式
估计 $f (x)$ 中待定系数

拟合要解决的两个基本问题

确定 $f (x)$ 表达式的形式

数据可视化，直观判断 $f (x)$ 的形式
通过机理分析来确定 $f (x)$ 的形式
饮酒驾车中酒精浓度拟合问题
$y = ke^{-at}$
SARS的传播规律拟合问题
$i(t)=\frac{1}{1+\left( \frac{1}{i_{0}}-1 \right)e^{-\lambda t}}$
估计 $f (x)$ 中的待定参数
根据 $f (x)$ 表达式中，待定参数具体形式的不同，一般有两种方法
线性最小二乘法：适用于 $f (x)$ 中待定参数为线性形式
非线性最小二乘法：适用于 $f (x)$ 中待定参数为非线性形式

线性最小二乘拟合及其实现

线性最小二乘拟合

基本思路

选定一组基本函数 $r_{1}(x),r_{2}(x),\dots,r_{m}(x);m<n$ ，确定 $f (x)$ 表达式
$f(x)=a_{1}r_{1}(x)+a_{2}r_{2}(x)+\dots+a_{m}r_{m}(x)$
其中， $a_{1},a_{2},\dots,a_{m}$ 为待定参数
按照最小二乘原则确定待定系数
求 $a_{1},a_{2},\dots,a_{m}$ 使n个点 $x_{i},y_{i})$ 与曲线 $f (x)$ 的距离平方和 $J(a_{1},a_{2},\dots,a_{m})$ 最小
$min_{a_{1},a_{2},\dots,a_{m}}J(a_{1},a_{2},\dots,a_{m})=\sum_{i=1}^{n}(\delta_{i})^{2}$
$=\sum_{i=1}^{n}(y_{i}-a_{1}r_{1}x_{i}-a_{1}r_{1}x_{i}-\dots-a_{m}r_{m}x_{i})^{2}$

优化问题的求解
超定方程组，方程个数大于未知量个数的方程组
$\left\{\begin{matrix} r_{11}a_{1}+r_{12}a_{2}+\dots+r_{1m}a_{m}=y_{1} \\ r_{21}a_{1}+r_{22}a_{2}+\dots+r_{2m}a_{m}=y_{2} \\ \dots\dots \\ r_{n1}a_{1}+r_{n2}a_{2}+\dots+r_{nm}a_{m}=y_{n} \end{matrix}\right.$
$R=\begin{pmatrix} r_{11}&&r_{12}&&\dots&&r_{1m} \\ r_{21}&&r_{22}&&\dots&&r_{2m} \\ \dots&&\dots&&\dots&&\dots \\ r_{n1}&&r_{n2}&&\dots&&r_{nm} \end{pmatrix}$
$a=\begin{pmatrix} a_{1} \\ a_{2} \\ \dots \\ a_{m} \end{pmatrix}\qquad y=\begin{pmatrix} y_{1} \\ y_{2} \\ \dots \\ y_{n} \end{pmatrix}$
超定方程组可改写为： $R a = y$
一般地，超定方程组是不存在解的矛盾方程组，不存在准确解
需要在最小二乘的意义下求解
$\sum_{i=1}^{n}(\delta_{i})^{2}=\sum_{i=1}^{n}(y_{i}-a_{1}r_{1}x_{i}-a_{1}r_{1}x_{i}-\dots-a_{m}r_{m}x_{i})^{2}$
最小化的解 $a_{1},a_{2},\dots,a_{m}$ ，称为超定方程组的线性最小二乘解

线性最小二乘的求解

当 $R^{T}R$ 可逆时，超定方程组 $R a = y$ 存在最小二乘解，且为如下正规方程组的解
$R^{T}Ra=R^{T}y$
$a=(R^{T}R)^{-1}R^{T}y$
线性最小二乘的Matlab实现
对超定方程组 $R_{n\times m}a_{m\times_{1}}=y_{n\times_{1}}$

a = R \ y

即可得到相应的最小二乘解

例

![[Pasted image 20240819095322.png]]

求二次多项式
$f(x)=a_{1}x^{2}+a_{2}x+a_{3}$
中的系数 $a_{1},a_{2},a_{3}$
$\begin{pmatrix} x_{1}^{2} &&x_{1}&&1\\ x_{2}^{2} &&x_{2}&&1\\ \dots &&\dots&&\dots\\ x_{11}^{2}&&x_{11}&&1 \end{pmatrix}=\begin{pmatrix} a_{1} \\ a_{2} \\ a_{3} \end{pmatrix}=\begin{pmatrix} y_{1} \\ y_{2} \\ \dots \\ y_{11} \end{pmatrix}$

x = 0 : 0.1 : 1;
y = [-0.447, 1.978, 3.28, 6.16, 7.08, 7.34, 7.66, 9.58, 9.48, 9.30, 11.2];
R = [(x.^2)',x',ones(11,1)];
A=R\y'

计算结果

A = [-9.8108, 20.1293, -0.0317]

$f(x)=-9.8108x^{2}+20.1293x-0.0317$

非线性最小二乘拟合及其实现

基本思路

确定问题的函数表达式的形式
$y=f(a_{1},a_{2},\dots,a_{m};x)$
其中 $a_{1},a_{2},\dots,a_{m}$ 为待定参数
按照最小二乘原则确定待定参数
在最小二乘意义下，求使得
$min_{a_{1},a_{2},\dots,a_{m}}J(a_{1},a_{2},\dots,a_{m})=\sum_{i=1}^{n}(\delta_{i})^{2}$
$=\sum_{i=1}^{n}(y_{i}-f(a_{1},a_{2},\dots,a_{m};x_{i}))^{2}$
最小化的解 $\hat{a_{1}},\hat{a_{2}},\dots,\hat{a_{m}}$ ，称为待定系数 $a_{1},a_{2},\dots,a_{m}$ 的非线性最小二乘解

非线性最小二乘拟合的Matlab实现

两个非线性最小二乘拟合函数

lsqcurvefit
lsqnonlin

非线性最小二乘求解命令lsqcurvefit
已知数据点
$xdata=(xdata_{1},xdata_{2},\dots,xdata_{n})^{T}$
$ydata=(ydata_{1},ydata_{2},\dots,ydata_{n})^{T}$
lsqcurvefit通过求解如下最小二乘问题，得到参数向量x的估计
$|_{2}^{2}=min_{x}\sum_{i}(ydata_{i}-F(x,xdata_{i})^{2}$
其中， $F (x, x d a t a)$ 表示预测值向量
$xdata)=\begin{pmatrix} F(x,xdata_{1}) \\ F(x,xdata_{2}) \\ \dots \\ F(x,xdata_{n}) \end{pmatrix}$
基本格式

x = lsqcurvefit('fun', x0, xdata, ydata);

fun是预先定义的函数，实现待拟合非线性函数 $F (x, x d a t a)$
x0给定参数初始值(lsqcurvefit求解非线性优化问题的初始迭代点)
xdata和ydata是样本数据的输入和输出向量(或矩阵)
x为未知参数的非线性最小二乘估计

例

![[Pasted image 20240819133318.png]]

$xdata)=\begin{pmatrix} F(x,xdata_{1}) \\ F(x,xdata_{2}) \\ \dots \\ F(x,xdata_{n}) \end{pmatrix}$

x，三个待定参数组成的向量
xdata，所有的自变量的取值组成的向量
输出，是待拟合的函数在每一个样本点处对应的拟合值
$=\begin{pmatrix} a+be^{-0.02kt_{1}} \\ a+be^{-0.02kt_{2}} \\ \dots \\ a+be^{-0.02kt_{11}} \end{pmatrix}$

function f = curvefun(x, tdata)
f = x(1) + x(2)*exp(-0.02*x(3)*tdata) %其中x(1)=a;x(2)=b;x(3)=k

%主程序命令
tdata = 100:100:1000
ydata = 1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];
x0 = [0.1; 0.1; 0.1];
x = lsqcurvefit('curvefun', x0, xdata, ydata)
y_pred = curvefun(x, tdata)
plot(tdata, ydata, 'bo', tdata, y_pred, 'r.:') %可视化预测结果

运算结果

x = 0.0070  -0.0030  0.1012
y_pred = 0.0045  0.0050  0.0054  0.0057  0.0059  0.0061  0.0063  0.0064  0.0065  0.0066

结论
a = 0.0070, b=-0.0030, k=1012
$f(t)=0.0070-0.0030e^{-0.02\times 0.1012\times t}$
![[Pasted image 20240819135525.png]]

一维插值及其实现

一维插值的定义

已知n+1个节点 $(x_{i},y_{i})；i= 0,1,\dots,n$ ，其中 $x_{i}$ 互不相同，设
$x_{0}<x_{1}<\dots<x_{n}=b$
构建一个函数 $y = f (x)$ 通过全部节点，即
$f(x_{i})=y_{i};\quad i=0,1,\dots,n$
用函数 $f (x)$ 计算任一被插值点 $x^{*}$ 处的插值，即
$y^{*}=f(x^{*})$
![[Pasted image 20240819140124.png]]

拉格朗日插值

已知n+1个节点 $(x_{i},y_{i})；i= 0,1,\dots,n$ ，构造一个n次多项式 $P_{n}(x)$ ，满足
$P_{n}(x_{i})=y_{i};\quad i=1,2,\dots,n$
被称为拉格朗日插值多项式，构造方法如下
$P_{n}(x)=\sum_{i=1}^{n}L_{i}(x)y_{i}$
其中， $L_{i}(x)$ 为n次多项式
$L_{i}(x)=\frac{(x-x_{0})(x-x_{1})\dots(x-x_{i-1})(x-x_{i+1})\dots(x-x_{n})}{(x_{i}-x_{0})(x_{i}-x_{1})\dots(x_{i}-x_{i-1})(x_{i}-x_{i+1})\dots(x_{i}-x_{n})}$
称为拉格朗日插值基函数
$L_{i}(x_{j})= \left\{\begin{matrix} 0,\ i\ne j \\ 1,\ i=j \end{matrix}\right.$
$P_{n}(x_{j})=y_{j}\quad j=0,1,\dots,n$
两点一次线性插值多项式
$P_{1}(x)=\frac{(x-x_{1})}{(x_{0}-x_{1})}y_{0}+\frac{(x-x_{0})}{(x_{1}-x_{0})}y_{1}$
三点二次抛物插值多项式
$P_{2}(x)=\frac{(x-x_{1})(x-x_{2})}{(x_{0}-x_{1})(x_{0}-x_{2})}y_{0}+\frac{(x-x_{0})(x-x_{2})}{(x_{1}-x_{0})(x_{1}-x_{2})}y_{1}+\frac{(x-x_{0})(x-x_{1})}{(x_{2}-x_{0})(x_{2}-x_{1})}y_{2}$
有可能产生震荡现象

例
![[Pasted image 20240819163240.png]]
![[Pasted image 20240819163348.png]]
![[Pasted image 20240819163404.png]]
![[Pasted image 20240819163423.png]]
![[Pasted image 20240819163435.png]]
![[Pasted image 20240819163450.png]]

随着节点个数的增加，多项式阶数不断增大， $P_{n}(x)$ 产生较大的震荡现象

分段线性插值

基本思想：构造一条一条的直线段进行插值计算
$L_{n}(x)=\sum_{i=1}^{n}l_{i}(x)y_{i}$
其中
$l_{i}(x)= \left\{\begin{matrix} \frac{x-x_{i-1}}{(x_{i}-x_{i-1})},\qquad x_{i-1}\le x<x_{i} \\ \frac{x-x_{i+1}}{(x_{i}-x_{i+1})},\qquad x_{i}\le x<x_{i+1} \\ 0,\qquad otherwise \end{matrix}\right.$

思路简单易实现
计算量小
不会产生震荡现象
光滑性不好

三次样条插值

基本思想
构造一段一段的三次多项式曲线进行插值计算
$s_{i}(x)=a_{i}(x)^{3}+b_{i}(x)^{2}+c_{i}x+d_{i},\ i=1,2,\dots,n$

![[Pasted image 20240819164059.png]]

三次样条插值函数
$S(x)=\left \{ s_{i}(x),x\in [x_{i-1},x_{i}]| i=1,2,\dots,n \right \}$
利用如下三类条件

插值条件n+1
$S(x_{i})=y_{i},\ i=0,1,2,\dots,n$
二阶光滑3n-3
$S(x)\in C^{2}[x_{0},x_{n}]$
自然边界条件2
$S''(x_{0})=S''(x_{n})=0$
可求得三次样条插值函数的所有待定系数

一维插值的Matlab实现

yi = interp1(x, y, xi, 'method')

x，y，插值节点对应的坐标
xi，被插值节点的坐标
yi，与xi对应的插值结果
‘method’，插值方法
nearest：最邻近插值
linear：分段线性插值
spline：三次样条插值
cubic：立方插值
缺省时：分段线性插值
x是单调的
xi不能够超过x的范围(内插)

例

![[Pasted image 20240819165803.png]]

hours = 1:12;
temps =[5 8 9 15 25 29 31 30 22 25 27 24];
h= 1 : 0.1 : 12；
t_nearest = interp1(hours, temps, h, 'nearest');%邻近插值
t_linear = interp1(hours, temps, h, 'linear');%分段线性插值
t_spline = interp1(hours, temps, h, 'spline');%三次样条插值

hours，原始的插值节点，横坐标的取值向量
temps，因变量温度的取值
h，被插值节点对应的向量，每隔1/10小时产生一个时间点对应的等差数列

subplot(2,2,1), plot(hours, temps, 'bo')
subplot(2,2,2), plot(hours, temps, 'bo', h, t_nearest,'r:')%邻近插值
subplot(2,2,3), plot(hours, temps, 'bo', h, t_linear,  'r:')%线性插值
subplot(2,2,4), plot(hours, temps, 'bo', h, t_spline, 'r:')%三次样条插值

![[Pasted image 20240819170440.png]]

二维插值及其实现

二维插值定义

网格节点
已知mxn个节点
$(x_{i},y_{i},z_{ij});\ i=1,\dots,m;j=1,\dots,n$
其中 $x_{i},y_{i}$ 互不相同，设
$x_{1}<x_{2}<\dots<x_{m}=b$
$c=y_{1}<y_{2}<\dots<y_{n}=d$
构造一个二元函数 $z = f (x, y)$ ，通过已知节点，即
$f(x_{i},y_{j})=z_{ij};\ i=1,\dots,m;\ j=1,\dots,n$
并由此计算插值
$z^{*}=f(x^{*},y^{*})$
散乱节点
已知n个节点
$(x_{i},y_{i},z_{i});\ i=1,\dots,n$
其中 $x_{i},y_{i}$ 互不相同
构造一个二元函数 $z = f (x, y)$ ，通过已知节点，即
$f(x_{i},y_{j})=z_{i};\ i=1,\dots,n$
并由此计算插值
$z^{*}=f(x^{*},y^{*})$

二维插值方法

最邻近插值
基本思想
将与被插值点最邻近节点的函数值作为插值结果

思路直观，易理解
计算量小，易于实现
不连续

分片线性插值
基本思想
构造一片一片的空间三角平面进行插值计算

直观，易于实现
连续
光滑性不好

双线性插值
基本思想
构造一片一片的双线性空间二次曲面进行插值计算
双线性插值函数表达式为
$f (x, y) = (a x + b) (cy + d)$
四个待定系数的确定：
利用矩形四个顶点(插值节点)处的函数值，建立四个方程进行求解

易理解实现
函数连续
光滑性不好

三次卷积插值，三次样条插值
保证光滑

二维插值的Matlab实现

网格节点

z = interp2(x0, y0, z0, x, y, 'method')

x0，y0，z0，指定插值节点对应的坐标
x，y，指定被插值节点的信息
z，被插值节点的插值结果
method，插值方法
nearest：最邻近插值
linear：线性插值 $C^{0}$
spline：三次样条插值 $C^{2}$
cubic：三次卷积插值 $C^{1}$
缺省时：线性插值
x0，y0单调
x，y为同型矩阵，或一个行向量另一个列向量
x，y不能够超过x0，y0的范围(内插)

散乱节点

z = gridata(x0, y0, z0, x, y, 'method')

x0，y0，z0，指定插值节点对应的坐标
x，y，指定被插值节点的信息
z，被插值节点的插值结果
method，插值方法
nearest：最邻近插值
linear：线性插值 $C^{0}$
cubic：三次插值 $C^{2}$
v4：双调和样条插值 $C^{2}$
缺省时：线性插值
x，y为同型矩阵，或一个行向量另一个列向量
x，y不能够超过x0，y0的范围(内插)

例1

![[Pasted image 20240819180032.png]]

画出原始数据的温度分布图

x=1:5;
y=1:3;
temps=[82 81 80 82 84;
	79 63 61 65 81;
	84 84 82 85 86];
mesh(x,y, temps)

![[Pasted image 20240819213301.png]]

然后，在x，y方向上每隔0.2个单位处进行插值，
再输入以下命令：

xi = 1: 0.2 : 5;
yi = 1: 0.2 : 3;
zi = interp2(x, y, temps, xi', yi, 'cubic');
mesh(xi, yi, zi)

画出插值后的温度分布曲面图
![[Pasted image 20240819213628.png]]

例2

![[Pasted image 20240819213644.png]]

画出原始海拔高度数据分布图

x=1200: 400: 4000;
y=[1200 1600 2400 2800 3200 3600];
[X,Y] = meshgrid(x,y);
Z=[1130 1250 1280 1230 1040 900 500 700;
1320 1450 1420 1400 1300 700 900 850;
1390 1500 1500 1400 900 1100 1060 950;
1500 1200 1100 1350 1450 1200 1150 1010;
1500 1550 1600 1550 1600 1600 1600 1550;
1480 1500 1550 1510 1430 1300 1200 980];
surfc(X, Y, Z);

meshgrid，对分割向量进行网格交叉，得到网格矩阵
surfc，三维空间曲面绘制

然后，对x，y方向上的取值向量加密，定义被插值数据点

xi = linspace(1200, 4000, 50);
yi = linspace(1200, 3600, 50);
[XI,YI]=meshgrid(xi,yi);

linspace，从1200到4000，产生50个均匀分布的点，构成等差数列
进行插值计算，并可视化

methods = {'nearest', 'linear', 'spline', 'cubic'};
for i = 1 : 4
	ZI = interp2(X, Y, Z, XI, YI, methods{i});
	subplot(2, 2, i);
	surfc(XI, YI, ZI);%这个函数直接就是把等高线画在下面
	title(['方法:',methods{i}]);
	xlabel('x');
	ylabel('y');
end

![[Pasted image 20240819220510.png]]

ChoSeitaku

关注

7
点赞
踩
6

收藏

觉得还不错? 一键收藏
0
评论
拟合与插值|线性最小二乘拟合|非线性最小二乘拟合|一维插值|二维插值

挖掘数据背后的规律是数学建模的重要任务，拟合与插值是常用的分析方法已知一组数据(以二维为例)，即平面上n个点(xi,yi),i=1,2,…,n(x_{i},y_{i}),i =1,2,\dots,n(xi,yi),i=1,2,…,n，寻求一个函数y=f(x)y=f(x)y=f(x)，使得f(x)f(x)f(x)在某种准则下与所有数据点最为接近∑i=1n(δi)2→min\sum_{i=1}^{n}(\delta_{i})^{2}\to mini=1∑n(δi)2→min确定f(x)f(x)f(
复制链接

扫一扫

专栏目录