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

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

  • 掌握拟合与插值的基本概念和方法
  • 熟悉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]]

∑ i = 1 n ( δ i ) 2 → m i n \sum_{i=1}^{n}(\delta_{i})^{2}\to min i=1n(δi)2min

  1. 确定 f ( x ) f(x) f(x)表达式的形式
  2. 估计 f ( x ) f(x) f(x)中待定系数
拟合要解决的两个基本问题

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

  1. 数据可视化,直观判断 f ( x ) f(x) f(x)的形式
    ![[Pasted image 20240819090140.png]]
    ![[Pasted image 20240819090155.png]]

  2. 通过机理分析来确定 f ( x ) f(x) f(x)的形式
    饮酒驾车中酒精浓度拟合问题
    y = k e − a t y = ke^{-at} y=keat
    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+(i011)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)中待定参数为非线性形式

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

线性最小二乘拟合

基本思路

  1. 选定一组基本函数 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为待定参数
  2. 按照最小二乘原则确定待定系数
    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=1n(δ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=1n(yia1r1xia1r1xiamrmxi)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= r11r21rn1r12r22rn2r1mr2mrnm
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= a1a2am y= y1y2yn
超定方程组可改写为: 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=1n(δi)2=i=1n(yia1r1xia1r1xiamrmxi)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]]

求二次多项式
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} x12x22x112x1x2x11111 = a1a2a3 = y1y2y11

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.1293x0.0317

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

基本思路
  1. 确定问题的函数表达式的形式
    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为待定参数
  2. 按照最小二乘原则确定待定参数
    在最小二乘意义下,求使得
    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=1n(δ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=1n(yif(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∣∣ydataF(x,xdata)22=minxi(ydataiF(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]]

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+be0.02kt1a+be0.02kt2a+be0.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.00700.0030e0.02×0.1012×t
![[Pasted image 20240819135525.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]]

拉格朗日插值

已知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=1nLi(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)=(xix0)(xix1)(xixi1)(xixi+1)(xixn)(xx0)(xx1)(xxi1)(xxi+1)(xxn)
称为拉格朗日插值基函数
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)=(x0x1)(xx1)y0+(x1x0)(xx0)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)=(x0x1)(x0x2)(xx1)(xx2)y0+(x1x0)(x1x2)(xx0)(xx2)y1+(x2x0)(x2x1)(xx0)(xx1)y2
有可能产生震荡现象


![[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 ) 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=1nli(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)= (xixi1)xxi1,xi1x<xi(xixi+1)xxi+1,xix<xi+10,otherwise

  • 思路简单易实现
  • 计算量小
  • 不会产生震荡现象
  • 光滑性不好
    ![[Pasted image 20240819163641.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]]

三次样条插值函数
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[xi1,xi]i=1,2,,n}
利用如下三类条件

  1. 插值条件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
  2. 二阶光滑3n-3
    S ( x ) ∈ C 2 [ x 0 , x n ] S(x)\in C^{2}[x_{0},x_{n}] S(x)C2[x0,xn]
  3. 自然边界条件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]]

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]]

二维插值及其实现

二维插值定义
  1. 网格节点
    已知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]]

  2. 散乱节点
    已知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]]

二维插值方法
  1. 最邻近插值
    基本思想
    将与被插值点最邻近节点的函数值作为插值结果
    ![[Pasted image 20240819174727.png]]
  • 思路直观,易理解
  • 计算量小,易于实现
  • 不连续
  1. 分片线性插值
    基本思想
    构造一片一片的空间三角平面进行插值计算
    ![[Pasted image 20240819174905.png]]
  • 直观,易于实现
  • 连续
  • 光滑性不好
  1. 双线性插值
    基本思想
    构造一片一片的双线性空间二次曲面进行插值计算
    双线性插值函数表达式为
    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]]
  • 易理解实现
  • 函数连续
  • 光滑性不好
  1. 三次卷积插值,三次样条插值
    保证光滑
二维插值的Matlab实现
  1. 网格节点
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的范围(内插)
  1. 散乱节点
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]]

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

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,三维空间曲面绘制
    ![[Pasted image 20240819215809.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 : 4ZI = 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]]

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/877489.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

C语言指针详解-上

C语言指针详解-上 前言1.指针的基本概念1.1指针是什么1.2指针的声明与初始化1.3取地址符&和解引用符*& 运算符用于**获取变量的地址*** 运算符用于访问指针指向的值 2.指针的类型常见数据类型的指针指针与数组、字符串数组指针结构体指针函数指针二级指针void指针 3.指…

对零基础想转行网络安全同学的一点建议

最近有同学在后台留言&#xff0c;0基础怎么学网络安全&#xff1f;0基础可以转行做网络安全吗&#xff1f;以前也碰到过类似的问题&#xff0c;想了想&#xff0c;今天简单写一下。 我的回答是先了解&#xff0c;再入行。 具体怎么做呢&#xff1f; 首先&#xff0c;你要确…

滑动变阻器的未来发展趋势和前景如何?是否有替代品出现?

滑动变阻器是常见的电子元件&#xff0c;主要用于调节电路中的电阻值。随着科技的不断发展&#xff0c;滑动变阻器的未来发展趋势和前景也引起了广泛关注。 滑动变阻器的未来发展将更加注重智能化&#xff0c;随着物联网、人工智能等技术的快速发展&#xff0c;滑动变阻器也将与…

C语言 | Leetcode C语言题解之第347题前K个高频元素

题目&#xff1a; 题解&#xff1a; struct hash_table {int key;int val;// 查看 https://troydhanson.github.io/uthash/ 了解更多UT_hash_handle hh; };typedef struct hash_table* hash_ptr;struct pair {int first;int second; };void swap(struct pair* a, struct pair*…

YUM和NFS

文章目录 yum软件仓库的提供方式RPM软件包的来源Linux系统各家厂商用的安装源命令---yum 配置本地yum源具体操作 搭建ftp yum仓库环境具体操作实操环境服务端一、安装 vsftpd服务二、创建一个文件&#xff0c;并且挂载三、开启服务四、查看挂载 客户端五、备份六、搭建ftp yum仓…

【联想电脑】:使用拓展坞后转接HDMI,无法识别显示屏

项目场景&#xff1a; 作为一个嵌入式软件开发者&#xff0c;有两个外接屏幕&#xff0c;不足为奇。 但是在今天的使用电脑过程中&#xff0c;出现了接了一个拓展坞上面有HDMI接口&#xff0c;但是HDMI接口接上外接显示屏的时候电脑无法识别到&#xff0c;导致只有电脑直连的HD…

使用Docker-compose一键部署Wordpress平台

一、Docker-compose概述&#xff1a; docker-compose&#xff1a;单机容器编排 Dockerfile&#xff1a;先配置好文件&#xff0c;然后build&#xff0c;镜像——>容器。 docker-compose&#xff1a;即可基于Dockerfile&#xff0c;也可以基于镜像&#xff0c;可以一键式拉…

安全基础学习-SM4加密算法

SM4 是一种中国国家密码标准(GB/T 32907-2016)中定义的分组加密算法,又称为“中国商用密码算法SM4”。它是由中国国家密码管理局发布的,并广泛应用于金融、电子商务和其他需要数据加密的场景。 1、SM4 算法概述 SM4 是一种对称加密算法,意味着加密和解密使用相同的密钥。…

使用WINUI3 编写一个小软件1 C#

本篇主要是记录安装和运行的问题。 先说安装 因为我是WIN11&#xff0c;所以勾了&#xff0c;如果你是WIN10就不用勾选11那个&#xff0c;但是我不确定用11要不要10那个&#xff0c;所以就勾了&#xff0c;按安装手册来的。 2、创建项目 照着选就完事了&#xff0c;别选错 这…

【数据结构与算法】快速排序

快速排序目录 一.快速排序的原理二.快速排序的图解三.快速排序的实现1.基准两边分2.分而治之 四.完整代码 一.快速排序的原理 每次选取第一个数为基准数.然后使用乾坤大挪移将大于或者小于基准的元素分别放置于基准数两边.继续分别对基准数两侧未排序的数据使用分治法进行处理…

springboot项目配置https安装ssl证书教程

1.将下载的ssl证书文件中的jks后缀文件放在/src/main/resource文件夹里面 2.在配置文件中&#xff08;yml后缀配置文件的格式不同&#xff09;添加如下配置即可

Postman接口测试项目实战

第 1 章 什么是接口测试 1.1、为什么要进行接口测试 目前除了特别Low的公司外&#xff0c;开发都是前后端分离的&#xff0c;就是说前端有前端的工程师进行编码&#xff0c;后端有后端的工程师进行编码&#xff0c;前后端进行数据基本都是通过接口进行交互的。 1.2、接口测…

IDEA:如何在idea中设置自动导包

这里使用的是idea2020版本,但是不同版本操作不会有较大的差别. 在Editer中展开General之后,选中Auto Import,最后勾选中Add unambiguous imports on the fly.

Java Web|day5.MyBatis

MyBatis 定义 它是一款半自动的ORM持久层框架&#xff0c;具有较高的SQL灵活性&#xff0c;支持高级映射(一对一&#xff0c;一对多)&#xff0c;动态SQL&#xff0c;延迟加载和缓存等特性&#xff0c;但它的数据库无关性较低 **ORM: **Object Relation Mapping&#xff0c;…

培训学校课程管理系统-计算机毕设Java|springboot实战项目

&#x1f34a;作者&#xff1a;计算机毕设匠心工作室 &#x1f34a;简介&#xff1a;毕业后就一直专业从事计算机软件程序开发&#xff0c;至今也有8年工作经验。擅长Java、Python、微信小程序、安卓、大数据、PHP、.NET|C#、Golang等。 擅长&#xff1a;按照需求定制化开发项目…

Java面试--框架--Spring MVC

Spring MVC 目录 Spring MVC1.spring mvc简介2.spring mvc实现原理2.1核心组件2.2工作流程 3.RESTful 风格4.Cookie&#xff0c;Session4.1 会话4.2 保存会话的两种技术 5.拦截器5.1过滤器、监听器、拦截器的对比5.2 过滤器的实现5.3 拦截器基本概念5.4 拦截器的实现 1.spring …

如何使用unittest和pytest进行python脚本的单元测试

1. 关于unittest和pytest unittest是python内置的支持单元测试的模块&#xff0c;他提供了核心类&#xff0c;TestCase&#xff0c;让单元测试 代码的编写不再是从0开始&#xff0c;不再是作坊式&#xff0c;而是标准化&#xff0c;模板化&#xff0c;工厂化。 pytest是第三方…

【深度解析】WRF-LES与PALM微尺度气象大涡模拟

查看原文>>>【深度解析】WRF-LES与PALM微尺度气象大涡模拟 针对微尺度气象的复杂性&#xff0c;大涡模拟&#xff08;LES&#xff09;提供了一种无可比拟的解决方案。微尺度气象学涉及对小范围内的大气过程进行精确模拟&#xff0c;这些过程往往与天气模式、地形影响和…

Linux - 模拟实现 shell 命令行解释器

目录 简介 shell 的重要性 解释为什么学习 shell 的工作原理很重要 模拟实现一个简单的 shell 循环过程 1. 获取命令行 2. 解析命令行 3. 建立一个子进程&#xff08;fork&#xff09; 4. 替换子进程&#xff08;execvp&#xff09; 5. 父进程等待子进程退出&#xff08;wai…

合宙LuatOS AIR700 IPV6 TCP 客户端向NodeRed发送数据

为了验证 AIR700 IPV6 &#xff0c;特别新建向NodeRed Tcp发送的工程。 Air700发送TCP数据源码如下&#xff1a; --[[ IPv6客户端演示, 仅EC618系列支持, 例如Air780E/Air600E/Air780UG/Air700E ]]-- LuaTools需要PROJECT和VERSION这两个信息 PROJECT "IPV6_SendDate_N…