挖掘数据背后的规律是数学建模的重要任务，拟合与插值是常用的分析方法
- 掌握拟合与插值的基本概念和方法
- 熟悉Matlab相关程序实现
- 能够从数据中挖掘数学规律
拟合问题的基本提法
拟合问题的概念
已知一组数据(以二维为例)，即平面上n个点
(
x
i
,
y
i
)
,
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)在某种准则下与所有数据点最为接近
![![[Pasted image 20240819085354.png]]](https://i-blog.csdnimg.cn/direct/44a5f49ab5ce451ab9ad3d6d9a47dfbc.png)
∑ i = 1 n ( δ i ) 2 → m i n \sum_{i=1}^{n}(\delta_{i})^{2}\to min i=1∑n(δi)2→min
- 确定 f ( x ) f(x) f(x)表达式的形式
- 估计 f ( x ) f(x) f(x)中待定系数
拟合要解决的两个基本问题
确定 f ( x ) f(x) f(x)表达式的形式
-
数据可视化，直观判断 f ( x ) f(x) f(x)的形式
![![[Pasted image 20240819090140.png]]](https://i-blog.csdnimg.cn/direct/1b83f22283174f30b2bbd4de6bdeec08.png)
![![[Pasted image 20240819090155.png]]](https://i-blog.csdnimg.cn/direct/f1873640888a4290af427aa163321b3f.png)
-
通过机理分析来确定 f ( x ) f(x) f(x)的形式
饮酒驾车中酒精浓度拟合问题
y = k e − a t y = ke^{-at} y=ke−at
SARS的传播规律拟合问题
i ( t ) = 1 1 + ( 1 i 0 − 1 ) e − λ t i(t)=\frac{1}{1+\left( \frac{1}{i_{0}}-1 \right)e^{-\lambda t}} i(t)=1+(i01−1)e−λt1
估计 f ( x ) f(x) f(x)中的待定参数
根据 f ( x ) f(x) f(x)表达式中，待定参数具体形式的不同，一般有两种方法
线性最小二乘法：适用于 f ( x ) f(x) f(x)中待定参数为线性形式
非线性最小二乘法：适用于 f ( x ) f(x) f(x)中待定参数为非线性形式
线性最小二乘拟合及其实现
线性最小二乘拟合
基本思路
- 选定一组基本函数
r
1
(
x
)
,
r
2
(
x
)
,
…
,
r
m
(
x
)
;
m
<
n
r_{1}(x),r_{2}(x),\dots,r_{m}(x);m<n
r1(x),r2(x),…,rm(x);m<n，确定
f
(
x
)
f(x)
f(x)表达式
f ( x ) = a 1 r 1 ( x ) + a 2 r 2 ( x ) + ⋯ + a m r m ( x ) f(x)=a_{1}r_{1}(x)+a_{2}r_{2}(x)+\dots+a_{m}r_{m}(x) f(x)=a1r1(x)+a2r2(x)+⋯+amrm(x)
其中， a 1 , a 2 , … , a m a_{1},a_{2},\dots,a_{m} a1,a2,…,am为待定参数 - 按照最小二乘原则确定待定系数
求 a 1 , a 2 , … , a m a_{1},a_{2},\dots,a_{m} a1,a2,…,am使n个点 ( x i , y i ) (x_{i},y_{i}) (xi,yi)与曲线 f ( x ) f(x) f(x)的距离平方和 J ( a 1 , a 2 , … , a m ) J(a_{1},a_{2},\dots,a_{m}) J(a1,a2,…,am)最小
m i n a 1 , a 2 , … , a m J ( a 1 , a 2 , … , a m ) = ∑ i = 1 n ( δ i ) 2 min_{a_{1},a_{2},\dots,a_{m}}J(a_{1},a_{2},\dots,a_{m})=\sum_{i=1}^{n}(\delta_{i})^{2} mina1,a2,…,amJ(a1,a2,…,am)=i=1∑n(δi)2
= ∑ i = 1 n ( y i − a 1 r 1 x i − a 1 r 1 x i − ⋯ − a m r m x 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} =i=1∑n(yi−a1r1xi−a1r1xi−⋯−amrmxi)2
优化问题的求解
超定方程组，方程个数大于未知量个数的方程组
{
r
11
a
1
+
r
12
a
2
+
⋯
+
r
1
m
a
m
=
y
1
r
21
a
1
+
r
22
a
2
+
⋯
+
r
2
m
a
m
=
y
2
…
…
r
n
1
a
1
+
r
n
2
a
2
+
⋯
+
r
n
m
a
m
=
y
n
\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.
⎩
⎨
⎧r11a1+r12a2+⋯+r1mam=y1r21a1+r22a2+⋯+r2mam=y2……rn1a1+rn2a2+⋯+rnmam=yn
R
=
(
r
11
r
12
…
r
1
m
r
21
r
22
…
r
2
m
…
…
…
…
r
n
1
r
n
2
…
r
n
m
)
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}
R=
r11r21…rn1r12r22…rn2…………r1mr2m…rnm
a
=
(
a
1
a
2
…
a
m
)
y
=
(
y
1
y
2
…
y
n
)
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}
a=
a1a2…am
y=
y1y2…yn
超定方程组可改写为：
R
a
=
y
Ra=y
Ra=y
一般地，超定方程组是不存在解的矛盾方程组，不存在准确解
需要在最小二乘的意义下求解
∑
i
=
1
n
(
δ
i
)
2
=
∑
i
=
1
n
(
y
i
−
a
1
r
1
x
i
−
a
1
r
1
x
i
−
⋯
−
a
m
r
m
x
i
)
2
\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}
i=1∑n(δi)2=i=1∑n(yi−a1r1xi−a1r1xi−⋯−amrmxi)2
最小化的解
a
1
,
a
2
,
…
,
a
m
a_{1},a_{2},\dots,a_{m}
a1,a2,…,am，称为超定方程组的线性最小二乘解
线性最小二乘的求解
当
R
T
R
R^{T}R
RTR可逆时，超定方程组
R
a
=
y
Ra=y
Ra=y存在最小二乘解，且为如下正规方程组的解
R
T
R
a
=
R
T
y
R^{T}Ra=R^{T}y
RTRa=RTy
a
=
(
R
T
R
)
−
1
R
T
y
a=(R^{T}R)^{-1}R^{T}y
a=(RTR)−1RTy
线性最小二乘的Matlab实现
对超定方程组
R
n
×
m
a
m
×
1
=
y
n
×
1
R_{n\times m}a_{m\times_{1}}=y_{n\times_{1}}
Rn×mam×1=yn×1
a = R \ y
即可得到相应的最小二乘解
例
![![[Pasted image 20240819095322.png]]](https://i-blog.csdnimg.cn/direct/19f8053548ac4d4fa32a7aa93d3c05ab.png)
求二次多项式
f
(
x
)
=
a
1
x
2
+
a
2
x
+
a
3
f(x)=a_{1}x^{2}+a_{2}x+a_{3}
f(x)=a1x2+a2x+a3
中的系数
a
1
,
a
2
,
a
3
a_{1},a_{2},a_{3}
a1,a2,a3
(
x
1
2
x
1
1
x
2
2
x
2
1
…
…
…
x
11
2
x
11
1
)
=
(
a
1
a
2
a
3
)
=
(
y
1
y
2
…
y
11
)
\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}
x12x22…x112x1x2…x1111…1
=
a1a2a3
=
y1y2…y11
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.8108 x 2 + 20.1293 x − 0.0317 f(x)=-9.8108x^{2}+20.1293x-0.0317 f(x)=−9.8108x2+20.1293x−0.0317
非线性最小二乘拟合及其实现
基本思路
- 确定问题的函数表达式的形式
y = f ( a 1 , a 2 , … , a m ; x ) y=f(a_{1},a_{2},\dots,a_{m};x) y=f(a1,a2,…,am;x)
其中 a 1 , a 2 , … , a m a_{1},a_{2},\dots,a_{m} a1,a2,…,am为待定参数 - 按照最小二乘原则确定待定参数
在最小二乘意义下，求使得
m i n a 1 , a 2 , … , a m J ( a 1 , a 2 , … , a m ) = ∑ i = 1 n ( δ i ) 2 min_{a_{1},a_{2},\dots,a_{m}}J(a_{1},a_{2},\dots,a_{m})=\sum_{i=1}^{n}(\delta_{i})^{2} mina1,a2,…,amJ(a1,a2,…,am)=i=1∑n(δi)2
= ∑ i = 1 n ( y i − f ( a 1 , a 2 , … , a m ; x i ) ) 2 =\sum_{i=1}^{n}(y_{i}-f(a_{1},a_{2},\dots,a_{m};x_{i}))^{2} =i=1∑n(yi−f(a1,a2,…,am;xi))2
最小化的解 a 1 ^ , a 2 ^ , … , a m ^ \hat{a_{1}},\hat{a_{2}},\dots,\hat{a_{m}} a1^,a2^,…,am^，称为待定系数 a 1 , a 2 , … , a m a_{1},a_{2},\dots,a_{m} a1,a2,…,am的非线性最小二乘解
非线性最小二乘拟合的Matlab实现
两个非线性最小二乘拟合函数
- lsqcurvefit
- lsqnonlin
非线性最小二乘求解命令lsqcurvefit
已知数据点
x
d
a
t
a
=
(
x
d
a
t
a
1
,
x
d
a
t
a
2
,
…
,
x
d
a
t
a
n
)
T
xdata=(xdata_{1},xdata_{2},\dots,xdata_{n})^{T}
xdata=(xdata1,xdata2,…,xdatan)T
y
d
a
t
a
=
(
y
d
a
t
a
1
,
y
d
a
t
a
2
,
…
,
y
d
a
t
a
n
)
T
ydata=(ydata_{1},ydata_{2},\dots,ydata_{n})^{T}
ydata=(ydata1,ydata2,…,ydatan)T
lsqcurvefit通过求解如下最小二乘问题，得到参数向量x的估计
m
i
n
∣
∣
y
d
a
t
a
−
F
(
x
,
x
d
a
t
a
)
∣
∣
2
2
=
m
i
n
x
∑
i
(
y
d
a
t
a
i
−
F
(
x
,
x
d
a
t
a
i
)
2
min| |ydata-F(x,xdata)| |_{2}^{2}=min_{x}\sum_{i}(ydata_{i}-F(x,xdata_{i})^{2}
min∣∣ydata−F(x,xdata)∣∣22=minxi∑(ydatai−F(x,xdatai)2
其中，
F
(
x
,
x
d
a
t
a
)
F(x,xdata)
F(x,xdata)表示预测值向量
F
(
x
,
x
d
a
t
a
)
=
(
F
(
x
,
x
d
a
t
a
1
)
F
(
x
,
x
d
a
t
a
2
)
…
F
(
x
,
x
d
a
t
a
n
)
)
F(x, xdata)=\begin{pmatrix} F(x,xdata_{1}) \\ F(x,xdata_{2}) \\ \dots \\ F(x,xdata_{n}) \end{pmatrix}
F(x,xdata)=
F(x,xdata1)F(x,xdata2)…F(x,xdatan)
基本格式
x = lsqcurvefit('fun', x0, xdata, ydata);
- fun是预先定义的函数，实现待拟合非线性函数 F ( x , x d a t a ) F(x,xdata) F(x,xdata)
- x0给定参数初始值(lsqcurvefit求解非线性优化问题的初始迭代点)
- xdata和ydata是样本数据的输入和输出向量(或矩阵)
- x为未知参数的非线性最小二乘估计
例
![![[Pasted image 20240819133318.png]]](https://i-blog.csdnimg.cn/direct/0b6ab2059cb54df2852230468bac9a4e.png)
F ( x , x d a t a ) = ( F ( x , x d a t a 1 ) F ( x , x d a t a 2 ) … F ( x , x d a t a n ) ) F(x, xdata)=\begin{pmatrix} F(x,xdata_{1}) \\ F(x,xdata_{2}) \\ \dots \\ F(x,xdata_{n}) \end{pmatrix} F(x,xdata)= F(x,xdata1)F(x,xdata2)…F(x,xdatan)
- x，三个待定参数组成的向量
- xdata，所有的自变量的取值组成的向量
- 输出，是待拟合的函数在每一个样本点处对应的拟合值
= ( a + b e − 0.02 k t 1 a + b e − 0.02 k t 2 … a + b e − 0.02 k t 11 ) =\begin{pmatrix} a+be^{-0.02kt_{1}} \\ a+be^{-0.02kt_{2}} \\ \dots \\ a+be^{-0.02kt_{11}} \end{pmatrix} = a+be−0.02kt1a+be−0.02kt2…a+be−0.02kt11
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.0030
e
−
0.02
×
0.1012
×
t
f(t)=0.0070-0.0030e^{-0.02\times 0.1012\times t}
f(t)=0.0070−0.0030e−0.02×0.1012×t
![![[Pasted image 20240819135525.png]]](https://i-blog.csdnimg.cn/direct/cfb264ad513d41e8975bb1c5796528c0.png)
一维插值及其实现
一维插值的定义
已知n+1个节点
(
x
i
,
y
i
)
；
i
=
0
,
1
,
…
,
n
(x_{i},y_{i})；i= 0,1,\dots,n
(xi,yi)；i=0,1,…,n，其中
x
i
x_{i}
xi互不相同，设
a
=
x
0
<
x
1
<
⋯
<
x
n
=
b
a = x_{0}<x_{1}<\dots<x_{n}=b
a=x0<x1<⋯<xn=b
构建一个函数
y
=
f
(
x
)
y=f(x)
y=f(x)通过全部节点，即
f
(
x
i
)
=
y
i
;
i
=
0
,
1
,
…
,
n
f(x_{i})=y_{i};\quad i=0,1,\dots,n
f(xi)=yi;i=0,1,…,n
用函数
f
(
x
)
f(x)
f(x)计算任一被插值点
x
∗
x^{*}
x∗处的插值，即
y
∗
=
f
(
x
∗
)
y^{*}=f(x^{*})
y∗=f(x∗)
![![[Pasted image 20240819140124.png]]](https://i-blog.csdnimg.cn/direct/1cd6dba9308647e7aed7bf51c8dc7e05.png)
拉格朗日插值
已知n+1个节点
(
x
i
,
y
i
)
；
i
=
0
,
1
,
…
,
n
(x_{i},y_{i})；i= 0,1,\dots,n
(xi,yi)；i=0,1,…,n，构造一个n次多项式
P
n
(
x
)
P_{n}(x)
Pn(x)，满足
P
n
(
x
i
)
=
y
i
;
i
=
1
,
2
,
…
,
n
P_{n}(x_{i})=y_{i};\quad i=1,2,\dots,n
Pn(xi)=yi;i=1,2,…,n
被称为拉格朗日插值多项式，构造方法如下
P
n
(
x
)
=
∑
i
=
1
n
L
i
(
x
)
y
i
P_{n}(x)=\sum_{i=1}^{n}L_{i}(x)y_{i}
Pn(x)=i=1∑nLi(x)yi
其中，
L
i
(
x
)
L_{i}(x)
Li(x)为n次多项式
L
i
(
x
)
=
(
x
−
x
0
)
(
x
−
x
1
)
…
(
x
−
x
i
−
1
)
(
x
−
x
i
+
1
)
…
(
x
−
x
n
)
(
x
i
−
x
0
)
(
x
i
−
x
1
)
…
(
x
i
−
x
i
−
1
)
(
x
i
−
x
i
+
1
)
…
(
x
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})}
Li(x)=(xi−x0)(xi−x1)…(xi−xi−1)(xi−xi+1)…(xi−xn)(x−x0)(x−x1)…(x−xi−1)(x−xi+1)…(x−xn)
称为拉格朗日插值基函数
L
i
(
x
j
)
=
{
0
,
i
≠
j
1
,
i
=
j
L_{i}(x_{j})= \left\{\begin{matrix} 0,\ i\ne j \\ 1,\ i=j \end{matrix}\right.
Li(xj)={0, i=j1, i=j
P
n
(
x
j
)
=
y
j
j
=
0
,
1
,
…
,
n
P_{n}(x_{j})=y_{j}\quad j=0,1,\dots,n
Pn(xj)=yjj=0,1,…,n
两点一次线性插值多项式
P
1
(
x
)
=
(
x
−
x
1
)
(
x
0
−
x
1
)
y
0
+
(
x
−
x
0
)
(
x
1
−
x
0
)
y
1
P_{1}(x)=\frac{(x-x_{1})}{(x_{0}-x_{1})}y_{0}+\frac{(x-x_{0})}{(x_{1}-x_{0})}y_{1}
P1(x)=(x0−x1)(x−x1)y0+(x1−x0)(x−x0)y1
三点二次抛物插值多项式
P
2
(
x
)
=
(
x
−
x
1
)
(
x
−
x
2
)
(
x
0
−
x
1
)
(
x
0
−
x
2
)
y
0
+
(
x
−
x
0
)
(
x
−
x
2
)
(
x
1
−
x
0
)
(
x
1
−
x
2
)
y
1
+
(
x
−
x
0
)
(
x
−
x
1
)
(
x
2
−
x
0
)
(
x
2
−
x
1
)
y
2
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}
P2(x)=(x0−x1)(x0−x2)(x−x1)(x−x2)y0+(x1−x0)(x1−x2)(x−x0)(x−x2)y1+(x2−x0)(x2−x1)(x−x0)(x−x1)y2
有可能产生震荡现象
例
![![[Pasted image 20240819163240.png]]](https://i-blog.csdnimg.cn/direct/de8e3eef0d994764b5dd4460a441ef4d.png)
![![[Pasted image 20240819163348.png]]](https://i-blog.csdnimg.cn/direct/d7a0b4a8ee2d4ea8aeac83486cb63772.png)
![![[Pasted image 20240819163404.png]]](https://i-blog.csdnimg.cn/direct/b7dc82efa32d43d89c6c4048f5432420.png)
![![[Pasted image 20240819163423.png]]](https://i-blog.csdnimg.cn/direct/b55828712e2e461a98ab6023d575888a.png)
![![[Pasted image 20240819163435.png]]](https://i-blog.csdnimg.cn/direct/3de5fd6270664ad1af0818b0df0517bf.png)
![![[Pasted image 20240819163450.png]]](https://i-blog.csdnimg.cn/direct/92609c407f384a308a799d16f3876033.png)
随着节点个数的增加，多项式阶数不断增大， P n ( x ) P_{n}(x) Pn(x)产生较大的震荡现象
分段线性插值
基本思想：构造一条一条的直线段进行插值计算
L
n
(
x
)
=
∑
i
=
1
n
l
i
(
x
)
y
i
L_{n}(x)=\sum_{i=1}^{n}l_{i}(x)y_{i}
Ln(x)=i=1∑nli(x)yi
其中
l
i
(
x
)
=
{
x
−
x
i
−
1
(
x
i
−
x
i
−
1
)
,
x
i
−
1
≤
x
<
x
i
x
−
x
i
+
1
(
x
i
−
x
i
+
1
)
,
x
i
≤
x
<
x
i
+
1
0
,
o
t
h
e
r
w
i
s
e
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.
li(x)=⎩
⎨
⎧(xi−xi−1)x−xi−1,xi−1≤x<xi(xi−xi+1)x−xi+1,xi≤x<xi+10,otherwise
- 思路简单易实现
- 计算量小
- 不会产生震荡现象
- 光滑性不好
![![[Pasted image 20240819163641.png]]](https://i-blog.csdnimg.cn/direct/e77e0c6691f14629b2679157dfb4577c.png)
三次样条插值
基本思想
构造一段一段的三次多项式曲线进行插值计算
s
i
(
x
)
=
a
i
(
x
)
3
+
b
i
(
x
)
2
+
c
i
x
+
d
i
,
i
=
1
,
2
,
…
,
n
s_{i}(x)=a_{i}(x)^{3}+b_{i}(x)^{2}+c_{i}x+d_{i},\ i=1,2,\dots,n
si(x)=ai(x)3+bi(x)2+cix+di, i=1,2,…,n
![![[Pasted image 20240819164059.png]]](https://i-blog.csdnimg.cn/direct/1da0b97cc05d4488af1bc80d37bf1eec.png)
三次样条插值函数
S
(
x
)
=
{
s
i
(
x
)
,
x
∈
[
x
i
−
1
,
x
i
]
∣
i
=
1
,
2
,
…
,
n
}
S(x)=\left \{ s_{i}(x),x\in [x_{i-1},x_{i}]| i=1,2,\dots,n \right \}
S(x)={si(x),x∈[xi−1,xi]∣i=1,2,…,n}
利用如下三类条件
- 插值条件n+1
S ( x i ) = y i , i = 0 , 1 , 2 , … , n S(x_{i})=y_{i},\ i=0,1,2,\dots,n S(xi)=yi, i=0,1,2,…,n - 二阶光滑3n-3
S ( x ) ∈ C 2 [ x 0 , x n ] S(x)\in C^{2}[x_{0},x_{n}] S(x)∈C2[x0,xn] - 自然边界条件2
S ′ ′ ( x 0 ) = S ′ ′ ( x n ) = 0 S''(x_{0})=S''(x_{n})=0 S′′(x0)=S′′(xn)=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]]](https://i-blog.csdnimg.cn/direct/8b0a263353164cc595d95a359459966c.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]]](https://i-blog.csdnimg.cn/direct/b47c9df082d342389b2aaf38993cee7d.png)
二维插值及其实现
二维插值定义
-
网格节点
已知mxn个节点
( x i , y i , z i j ) ; i = 1 , … , m ; j = 1 , … , n (x_{i},y_{i},z_{ij});\ i=1,\dots,m;j=1,\dots,n (xi,yi,zij); i=1,…,m;j=1,…,n
其中 x i , y i x_{i},y_{i} xi,yi互不相同，设
a = x 1 < x 2 < ⋯ < x m = b a = x_{1}<x_{2}<\dots<x_{m}=b a=x1<x2<⋯<xm=b
c = y 1 < y 2 < ⋯ < y n = d c=y_{1}<y_{2}<\dots<y_{n}=d c=y1<y2<⋯<yn=d
构造一个二元函数 z = f ( x , y ) z=f(x,y) z=f(x,y)，通过已知节点，即
f ( x i , y j ) = z i j ; i = 1 , … , m ; j = 1 , … , n f(x_{i},y_{j})=z_{ij};\ i=1,\dots,m;\ j=1,\dots,n f(xi,yj)=zij; i=1,…,m; j=1,…,n
并由此计算插值
z ∗ = f ( x ∗ , y ∗ ) z^{*}=f(x^{*},y^{*}) z∗=f(x∗,y∗)
![![[Pasted image 20240819171122.png]]](https://i-blog.csdnimg.cn/direct/ec8b0111802c456bba676041a035e7bf.png)
-
散乱节点
已知n个节点
( x i , y i , z i ) ; i = 1 , … , n (x_{i},y_{i},z_{i});\ i=1,\dots,n (xi,yi,zi); i=1,…,n
其中 x i , y i x_{i},y_{i} xi,yi互不相同
构造一个二元函数 z = f ( x , y ) z=f(x,y) z=f(x,y)，通过已知节点，即
f ( x i , y j ) = z i ; i = 1 , … , n f(x_{i},y_{j})=z_{i};\ i=1,\dots,n f(xi,yj)=zi; i=1,…,n
并由此计算插值
z ∗ = f ( x ∗ , y ∗ ) z^{*}=f(x^{*},y^{*}) z∗=f(x∗,y∗)
![![[Pasted image 20240819174355.png]]](https://i-blog.csdnimg.cn/direct/b15621bfb7df4ae288d72170df1e9fac.png)
二维插值方法
- 最邻近插值
基本思想
将与被插值点最邻近节点的函数值作为插值结果
![![[Pasted image 20240819174727.png]]](https://i-blog.csdnimg.cn/direct/fee5dbcfd689478b862b510b1f4cb83d.png)
- 思路直观，易理解
- 计算量小，易于实现
- 不连续
- 分片线性插值
基本思想
构造一片一片的空间三角平面进行插值计算
![![[Pasted image 20240819174905.png]]](https://i-blog.csdnimg.cn/direct/8fd69e7d60464a4982cba14070f18572.png)
- 直观，易于实现
- 连续
- 光滑性不好
- 双线性插值
基本思想
构造一片一片的双线性空间二次曲面进行插值计算
双线性插值函数表达式为
f ( x , y ) = ( a x + b ) ( c y + d ) f(x,y)=(ax+b)(cy+d) f(x,y)=(ax+b)(cy+d)
四个待定系数的确定：
利用矩形四个顶点(插值节点)处的函数值，建立四个方程进行求解
![![[Pasted image 20240819175043.png]]](https://i-blog.csdnimg.cn/direct/37719b13d9704b259870022d42ba142c.png)
- 易理解实现
- 函数连续
- 光滑性不好
- 三次卷积插值，三次样条插值
保证光滑
二维插值的Matlab实现
- 网格节点
z = interp2(x0, y0, z0, x, y, 'method')
- x0，y0，z0，指定插值节点对应的坐标
- x，y，指定被插值节点的信息
- z，被插值节点的插值结果
- method，插值方法
nearest：最邻近插值
linear：线性插值 C 0 C^{0} C0
spline：三次样条插值 C 2 C^{2} C2
cubic：三次卷积插值 C 1 C^{1} C1
缺省时：线性插值 - 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 C^{0} C0
cubic：三次插值 C 2 C^{2} C2
v4：双调和样条插值 C 2 C^{2} C2
缺省时：线性插值 - x，y为同型矩阵，或一个行向量另一个列向量
- x，y不能够超过x0，y0的范围(内插)
例1
![![[Pasted image 20240819180032.png]]](https://i-blog.csdnimg.cn/direct/be28a2179f3946f4b9cd3d9e0ca308f0.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]]](https://i-blog.csdnimg.cn/direct/7c8934a68a5e47fab392e60a52061f7e.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]]](https://i-blog.csdnimg.cn/direct/cc89478fb94b4a4782660c85545355b4.png)
例2
![![[Pasted image 20240819213644.png]]](https://i-blog.csdnimg.cn/direct/6a64dfef3c354b36a8e52bf2d0774a6b.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，三维空间曲面绘制
![![[Pasted image 20240819215809.png]]](https://i-blog.csdnimg.cn/direct/78d7e8fd17584fcb8843258863c33519.png)
然后，对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]]](https://i-blog.csdnimg.cn/direct/a35923c1fefa4839afc99a206003c90e.png)
