（十）插值与拟合

来源：互联网发布：福特野马轮毂数据编辑：程序博客网时间：2024/06/05 09:04

插值与拟合

插值：求过已知有限个数据点的近似函数
拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小

插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显

插值方法

拉格朗日多项式插值

用多项式作为研究插值的工具，称为代数插值。其基本问题是：已知函数f(x)在区间[a,b]上n+1个不同点x0,x1,⋯,xn处的函数值yi=f(xi)(i=0,1,⋯,n)，求一个至多n 次多项式φ(x)使其在给定点处与f(x)同值，即满足插值条件。约束方程组的系数矩阵为A，则是范德蒙特(Vandermonde)行列式，只要n+1个节点互不相同，满足插值要求的插值多项式是唯一的。插值多项式与被插函数之间的差Rn(x)=f(x)−φn(x)称为截断误差，又称为插值余项
拉格朗日插值多项式
实际上比较方便的作法不是解方程组求待定系数，而是先构造一组基函数
$l i (x) = \prod j = 0 j \neq i n x - x j x i - x j$
li(x)是n次多项式，满足
$l i (x) = {0 j \neq i 1 j = i$
令
$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 x - x j x i - x j ⎞ ⎠ ⎟ ⎟ ⎟$
上式称为n次 Lagrange 插值多项式，由方程解的唯一性，
Lagrange 插值多项式存在唯一
用Matlab 作Lagrange 插值
Matlab 中没有现成的Lagrange 插值函数，必须编写一个M文件实现Lagrange插值。设n个节点数据以数组x0,y0输入（注意 Matlab的数组下标从 1 开始），m个插值点以数组x输入，输出数组y为m个插值。
```
function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end
```

牛顿（Newton）插值

导出Newton 公式前，先介绍公式表示中所需要用到的差商、差分的概念及性质

差商
称f[xi,xj]=f(xi)−f(xj)xi−xj为关于点xi,xj的一阶差商。一阶差商的差商为二阶差商。一般的，称f[x0,x1,⋯,xk−1]−f[x1,x1,⋯,xk]x0−xk为k阶差商。差商满足的性质：
- f[xi,xj]=f[xj,xi]
- f[xi,xj,xk]=f[xi,xk,xj]=f[xj,xi,xk]
Newton 插值公式
Newton 插值的优点是：每增加一个节点，插值多项式只增加一项,即Nn+1(x)=Nn(x)+(x−x0)⋯(x−xn)f[x0,x1,⋯,xn+1]因而便于递推运算。而且 Newton 插值的计算量小Lagrange 插值。
差商遇导数关系
f[x0,x1,⋯,xn+1]=f(n)(ζ)n!ζ∈(α,β)
其中α为自变量最小值，β为最大值
差分
当节点等距时，即相邻两个节点之差（称为步长）为常数，Newton 插值公式的形式会更简单。此时关于节点间函数的平均变化率（差商）可用函数值之差（差分）来表示
差分具有以下性质：
- 各阶差分均可表成函数值的线性组合
- 各种差分之间可以互化
等距节点插值公式
$N n (x 0 + t h) = f 0 + t △ f 0 + \dots + t ( t - 1 ) \dots ( t - n + 1 ) n ! △ n f 0$
式称为 Newton 向前插值公式

分段线性插值

插值多项式的振荡
用 Lagrange插值多项式Ln(x)近似f(x)(a≤x≤b)，虽然随着节点个数的增加，Ln(x)的次数n变大，多数情况下误差|Rn(x)|会变小。但是n增大时，Ln(x)的光滑性变坏，有时会出现很大的振荡。理论上，当n→∞，在[a,b]内并不能保证处处收敛于f(x).如
$f (x) = 1 1 + x 2, x \in [- 5, 5]$
高次插值多项式的这些缺陷，促使人们转而寻求简单的低次多项式插值
分段线性插值
简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数。用In(x)计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n无关。但n越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。

用Matlab 实现分段线性插值

Matlab 中有现成的一维插值函数interp1

y=interp1(x0,y0,x,'method')%method 指定插值的方法，默认为线性插值。其值可为：%'nearest' 最近项插值%'linear' 线性插值%'spline' 逐段3次样条插值%'cubic' 保凹凸性3次插值%所有的插值方法要求 x0 是单调的。当 x0 为等距时可以用快速插值法，使用快速插值法的格式为'*nearest'、'*linear'、'*spline'、'*cubic'。

埃尔米特(Hermite)插值

Hermite 插值多项式
如果对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值，这就是Hermite 插值问题。本节主要讨论在节点处插值函数与函数的值及一阶导数值均相等的Hermite 插值
Hermite 插值多项式为：
$H (x) = \sum i = 0 n h i [(x i - x (2 a i y i - y' i) + y i)]$
其中hi=∏j=0j≠i(x−xjxi−xj)2,ai=∑j=0j≠in1xi−xj
用Matlab 实现Hermite 插值
Matlab 中没有现成的Hermite 插值函数，必须编写一个M 文件实现插值.设n个节点的数据以数组x0（已知点的横坐标）， y0（函数值）， y1（导数值）输入（注意Matlab 的数组下标从1 开始），m 个插值点以数组x 输入，输出数组y 为m个插值。编写一个名为hermite.m 的M 文件：
```
function y=hermite(x0,y0,y1,x);n=length(x0);m=length(x);for k=1:myy=0.0;for i=1:nh=1.0;a=0.0;for j=1:nif j~=ih=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;a=1/(x0(i)-x0(j))+a;endendyy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));endy(k)=yy;end
```

样条插值

许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续，而且要有连续的曲率，这就导致了样条插值的产生

样条函数的概念
所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。
利用样条函数进行插值，即取插值函数为样条函数，称为样条插值。例如分段线性插值是一次样条插值。我们只介绍二次、三次样条插值见^196页^
三次样条插值在Matlab 中的实现
在 Matlab 中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件
```
y=interp1(x0,y0,x,'spline')；y=spline(x0,y0,x)；pp=csape(x0,y0,conds)，y=ppval(pp,x)
```
其中 x0,y0 是已知数据点，x 是插值点，y 是插值点的函数值。对于三次样条插值，我们提倡使用函数 csape，csape 的返回值是pp 形式，要求插值点的函数值，必须调用函数ppval。pp=csape(x0,y0)：使用默认的边界条件，即Lagrange 边界条件。pp=csape(x0,y0,conds)中的conds 指定插值的边界条件，其值可为：
- ’complete’ 边界为一阶导数，即默认的边界条件
- ‘not-a-knot’ 非扭结条件
- ‘periodic’ 周期条件
- ‘second’ 边界为二阶导数，二阶导数的值[0, 0]
- ‘variational’ 设置边界的二阶导数值为[0,0]
对于一些特殊的边界条件，可以通过 conds 的一个1×2矩阵来表示，conds 元素的取值为1，2。此时，使用命令
pp=csape(x0,y0_ext,conds)。其中y0_ext=[left, y0, right]，这里left 表示左边界的取值，right 表示右边界的取值。conds(i)=j 的含义是给定端点i 的j 阶导数，即conds 的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶
导数，对应的值由left 和right 给出。help csape解决

案例比较：197页

clc,clearx0=[0 3 5 7 9 11 12 13 14 15];y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];x=0:0.1:15;y1=lagrange(x0,y0,x); %调用前面编写的Lagrange插值函数y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,'spline');pp1=csape(x0,y0); y4=ppval(pp1,x);pp2=csape(x0,y0,'second'); y5=ppval(pp2,x);fprintf('比较一下不同插值方法和边界条件的结果:\n')fprintf('x y1 y2 y3 y4 y5\n')xianshi=[x',y1',y2',y3',y4',y5'];fprintf('%f\t%f\t%f\t%f\t%f\t%f\n',xianshi')subplot(2,2,1), plot(x0,y0,'+',x,y1), title('Lagrange')subplot(2,2,2), plot(x0,y0,'+',x,y2), title('Piecewise linear')subplot(2,2,3), plot(x0,y0,'+',x,y3), title('Spline1')subplot(2,2,4), plot(x0,y0,'+',x,y4), title('Spline2')dyx0=ppval(fnder(pp1),x0(1)) %求x=0处的导数ytemp=y3(131:151);index=find(ytemp==min(ytemp));xymin=[x(130+index),ytemp(index)]

B 样条函数插值方法

磨光函数

实际中的许多问题，往往是既要求近似函数（曲线或曲面）有足够的光滑性，又要求与实际函数有相同的凹凸性，一般插值函数和样条函数都不具有这种性质。如果对于一个特殊函数进行磨光处理生成磨光函数（多项式），则用磨光函数构造出样条函数作
为插值函数，既有足够的光滑性，而且也具有较好的保凹凸性，因此磨光函数在一维插值（曲线）和二维插值（曲面）问题中有着广泛的应用

由积分理论可知，对于可积函数通过积分会提高函数的光滑度，因此，我们可以利用积分方法对函数进行磨光处理。

定义：f(x)的k次磨光函数为

f k, h (x) = 1 h \int x + h 2 x - h 2 f k - 1, h (t) d t (k > 1)

事实上，磨光函数

fk,h(x)比

f(x)的光滑程度要高，且当磨光宽度

h很少时 ( ) ,

fk,h(x)很接近于

f(x)

等距B 样条函数

一维等距B 样条函数插值

二维等距B 样条函数插值

这是一种具有良好保凸性的光滑曲面（函数），在工程设计中是常用的，但只能使用于均匀划分或近似均匀划分的情况

二维插值

前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。若节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的高程（节点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这些点的高程（插值）

网格节点

已知m×n个节点：求点(x,y)处的插值z

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

其中 x0,y0 分别为m维和n维向量，表示网格节点。z0为n × m维矩阵，表示节点值，x,y为一维数组，表示插值点，x 与y 应是方向不同的向量，即一个是行向量，另一个是列向量，z 为矩阵，它的行数为y 的维数，列数为x 的维数，表示得到的插值，’method’
的用法同上面的一维插值。如果是三次样条插值，可以使用命令
pp=csape({x0,y0},z0,conds,valconds)，z=fnval(pp,{x,y})
其中 x0,y0 分别为m 维和n维向量，z0 为m × n 维矩阵，z 为矩阵，它的行数为x的维数，列数为y 的维数，表示得到的插值，具体使用方法同一维插值

clear,clcx=100:100:500;y=100:100:400;z=[636 697 624 478 450698 712 630 478 420680 674 598 412 400662 626 552 334 310];pp=csape({x,y},z')xi=100:10:500;yi=100:10:400cz1=fnval(pp,{xi,yi})cz2=interp2(x,y,z,xi,yi','spline')[i,j]=find(cz1==max(max(cz1)))x=xi(i),y=yi(j),zmax=cz1(i,j)

插值节点为散乱节点

ZI = GRIDDATA(X,Y,Z,XI,YI)

其中X、Y、Z 均为n 维向量，指明所给数据点的横坐标、纵坐标和竖坐标。向量XI、YI 是给定的网格点的横坐标和纵坐标，返回值ZI 为网格（XI，YI）处的函数值。XI与YI 应是方向不同的向量，即一个是行向量，另一个是列向量。

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];xi=75:1:200;yi=-50:1:150;zi=griddata(x,y,z,xi,yi','cubic')subplot(1,2,1), plot(x,y,'*')subplot(1,2,2), mesh(xi,yi,zi)

曲线拟合的线性最小二乘法

线性最小二乘法

无法知道 y 与x之间的关系，通常可以将数据(xi,yi)作图，直观地判断应该用什么样的曲线去作拟合。人们常用的曲线有

直线
多项式
双曲线
指数函数

对于指数曲线，拟合前需作变量代换，化为对a1,a2的线性函数

已知一组数据，用什么样的曲线拟合最好，可以在直观判断的基础上，选几种曲线分别拟合，然后比较，看哪条曲线的最小二乘指标J最小

最小二乘法的Matlab 实现

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\yx0=19:0.1:44;y0=ab(1)+ab(2)*x0.^2;plot(x,y,'o',x0,y0,'r')

多项式拟合方法

a=polyfit(x0,y0,m)

用m 次多项式拟合给定数据,多项式在x处的值y可用下面的函数计算y=polyval(a,x)

最小二乘优化

min F (x) = \sum i = 1 m f 2 i (x)

最小二乘优化是一类比较特殊的优化问题，在处理这类问题时，Matlab 也提供了一些强大的函数。在Matlab 优化工具箱中，用于求解最小二乘优化问题的函数有：lsqlin、lsqcurvefit、lsqnonlin、lsqnonneg，用法介绍如下

未完待续

曲线拟合的用户图形界面求法

Matlab 工具箱提供了命令cftool，该命令给出了一维数据拟合的交互式环境。具体执行步骤如下：

把数据导入到工作空间
运行cftool，打开用户图形界面窗口
对数据进行预处理
选择适当的模型进行拟合
生成一些相关的统计量，并进行预测

可以通过帮助（运行 doc cftool）熟悉该命令的使用细节

曲线拟合与函数逼近

如果已知一个较为复杂的连续函数y(x),x∈[a,b]，要求选择一个较简单的函数f(x)，在一定准则下最接近y(x)，就是所谓函数逼近。与曲线拟合的最小二乘准则相对应，函数逼近常用的一种准则是最小平方逼近。最简单的当然是用多项式逼近函数。

正交多项式

基底构成正交系

勒让得（Legendre）多项式
第一类切比雪夫（Chebyshev）多项式

案例：

f(x)=cos(x),x∈[−π2,π2]在H=Span{1,x2,x4}中的最佳平方逼近多项式

syms xbase=[1,x^2,x^4];y1=base.'*basey2=cos(x)*base.'r1=int(y1,-pi/2,pi/2)r2=int(y2,-pi/2,pi/2)a=r1\r2xishu1=double(a)digits(8),xishu2=vpa(a)

黄河小浪底调水调沙问题

未完待续

阅读全文

0 0