插值和拟合

插值和拟合都是函数逼近或者数值逼近的重要组成部分

         他们的共同点都是通过已知一些离散点集M上的约束，求取一个定义
         在连续集合S(M包含于S)的未知连续函数，从而达到获取整体规律的
         目的，即通过"窥几斑"来达到"知全豹"。

         简单的讲，所谓拟合是指已知某函数的若干离散函数值{f1,f2,…,fn}，通
         过调整该函数中若干待定系数f(λ1, λ2,…,λ3), 使得该函数与已知点集的
         差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或者
         线性回归(主要在统计中)，否则叫作非线性拟合或者非线性回归。表
         达式也可以是分段函数，这种情况下叫作样条拟合。


         而插值是指已知某函数的在若干离散点上的函数值或者导数信息，通
         过求解该函数中待定形式的插值函数以及待定系数，使得该函数在给
         定离散点上满足约束。插值函数又叫作基函数，如果该基函数定义在
         整个定义域上，叫作全域基，否则叫作分域基。如果约束条件中只有
         函数值的约束，叫作Lagrange插值，否则叫作Hermite插值。

         从几何意义上将，拟合是给定了空间中的一些点，找到一个已知形式
         未知参数的连续曲面来最大限度地逼近这些点；而插值是找到一个(
         或几个分片光滑的)连续曲面来穿过这些点。

  一、           概念的引入

1.   插值与拟合在现实生活中的应用

l         机械制造：汽车外观设计

l         采样数据的重新建构：电脑游戏中场景的显示，地质勘探，医学领域（CT）

2.   概念的定义

l         插值： 基于[a,b]区间上的n个互异点，给定函数f(x)，寻找某个函数去逼近f(x)。若要求φ(x)在xi处与f(xi)相等，这类的函数逼近问题称为插值问题，xi即是插值点

l         逼近： 当取值点过多时，构造通过所有点的难度非常大。此时选择一个次数较低的函数最佳逼近这些点，一般采用最小二乘法

l         光顾： 曲线的拐点不能太多，条件：

①二阶几何连续②不存在多余拐点③曲率变化较小

l         拟合：曲线设计过程中用插值或通过逼近方法是生成的曲线光滑(切变量连续)光顾

二、        插值理论

设函数y=f(x)在区间[a,b]上连续，在[a,b]上有互异点x₀,x₁,…,x_n处取值y₀,y₁,…,y_n 。如果函数φ(x)在点x_i上满足φ(x_i)=y_i (i=0,1,2,…,n)，则称φ(x)是函数y=f(x)的插值函数，x₀,x₁,…,x_n是插值节点。若此时φ(x)是代数多项式P(x)，则称P(x)为插值多项式。显然 f(x)≈φ(x)，x∈[a,b]

 

1.   拉格朗日插值

构造n次多项式P_n (x)=   y_k l_k (x)=y₀l₀ (x)+y₁l₁ (x)+…+y_nl_n (x)，这是不超过n次的多项式，其中基函数l_k(x)=

显然l_k (x)满足l_k (x_i)=

此时 P_n(x)≈f(x),误差R_n(x)=f(x)-P_n(x)=  

其中 ∈(a,b)且依赖于x， =（x-x₀）(x-x₁)…(x-x_n)

 

很显然，当n=1、插值节点只有两个x_k,x_k+1时

P₁(x)=y_kl_k(x)+y_k+1l_k+1(x)

其中基函数l_k(x)=   l_k+1(x)=

2.  牛顿插值

构造n次多项式N_n(x)=f(x₀)+f(x₀,x₁)(x-x₀)+f(x₀,x₁,x₂)(x-x₀)(x-x₁)+…

+f(x₀,x₁,x₂,…,x_n)(x-x₀)(x-x₁)…(x-x_n)

称为牛顿插值多项式，其中

               (二个节点，一阶差商)

          (三个节点，二阶差商)

  (n+1个节点，n阶差商)

注意：由于插值多项式的唯一性，有时为了避免拉格朗日余项R_n(x)中n+1阶导数的运算，用牛顿插值公式R_n (x)=f(x)-N_n(x)=f(x,x₀,…,x_n)ω_n+1(x),

其中ω_n+1(x)=(x-x₀)(x-x₁)…(x-x_n)

 

3.  分段插值------子区间内，避免函数在某些区间失真

1）  线性插值

已知n+1个不同节点x₀,x₁,…,x_n ,构造分段一次线性多项式P(x)，使之满足

l         P(x)在[a,b]上连续

l         P(x_k)=y_k

l         P(x)在[x_i,x_i+1]上是线性函数，P(x)=

2）  两点带导数插值---避免尖点、一阶连续

 

区间[a,b]上两个互异节点x_i,x_i+1，已知实数y _i,y _i+1,m _i,m _i+1，为了构造次数不大于3的多项式 满足条件  

引入 , 使之满足  

可以求出

此时 = + ，其中

 

4.        三次样条插值------二阶可导

对于给定n+1个不同节点x₀,x₁,…,x_n及函数值y₀,y₁,…,y_n，其中a=x₀<x₁<…<x_n=b。构造三次样条插值函数S(x)。S(x)称为三次样条函数时需满足：

l              S(x)在[a,b]上二阶导数连续

l              S(x_k)=y_k (k=0,1,…,n)

l              每个子区间[x_k,x_k+1]上S(x)是三次多项式(k=0,1,…,n)

 

5.   例题

已知函数y=f(x)的观测值

X	1	2	3	4
Y	0	-5	-6	3

求三次插值多项式φ(x)及φ(2.5).

解：

(1)拉格朗日插值

P₃(x)=y₀l₀(x)+y₁l₁(x)+y₂l₂(x)+y₃l₃(x)

=(-5)   + (-6)  +3

=x³-4x²+3

φ(x)≈P₃(x)= x³-4x²+3  φ(2.5)=2.5³-4*2.5²+3=-6.375

 

(2)牛顿插值

xi	f(xi)	一阶差商	二阶差商	三阶差商
1	0
2	–5	-5
3	–6	-1	2
4	3	9	5	1

 

N₃(x)=f(x₀)+f(x₀,x₁)(x-x₀)+f(x₀,x₁,x₂)(x-x₀)(x-x₁)+f(x₀,x₁,x₂,x₃)(x-x₀)(x-x₁)(x-x₂)

        =0+(-5)(x-1)+2*(x-1)(x-2)+1*(x-1)(x-2)(x-3)

=x³-4x²+3

 

 

三、        Matlab在插值中的应用

1.   Lagrange插值

 

１)方法回顾

 

对给定的n个节点x₁,x₂,…,x_n及对应的函数值y₁,y₂,…,y_n，利用n次Lagrange插值多项式公式，插值区间内任意x的函数值y可以通过下式求出：

2)  Matlab实现函数    Lagrange.m

function y=lagrange(x0,y0,x)

n=length(x0);m=length(x);

for i=1:m

z=x(i);

s=0.0;

for k=1:n

      p=1.0;

      for j=1:n

   if j~=k

        p=p*(z-x0(j))/(x0(k)-x0(j));

   end

      end

      s=p*y0(k)+s;

end

y(i)=s;

end

 

3) 例题的Matlab实现

x=[1 2 3 4];

y=[0 –5 –6 3];

lagrange(x,y,2.5)

 

2.  Runge现象及分段线性插值

１）Runge现象

Runge在本世纪初发现：在[-1,1]上用n+1个等距结点作插值多项式P_n(x)，使其在个结点的值与函数y(x)=1/(1+25x²)在结点的值相等。但在n→∞时，插值多项式P_n(x)在区间中部趋于y(x)。但对于0.726≤∣x∣≤1的x，P_n(x)严重发散。

通过下面的例子，以图形的方式体会Runge现象（ f(x)=1/(1+x²) ）

x=[-5:1:5];

y=1./(1+x.^2);

x0=[-5:0.1:5];

y0=lagrange(x,y,x0);

y1=1./(1+x0.^2);

%绘制图形

plot(x0,y0,'--r')

hold on

plot(x0,y1,'-b')

 

 

 

 

 

2）Matlab实现分段插值------一维插值interp1

l              yi=interp1(x,y,xi)          对(x,y)进行插值，计算插值点xi的函数值

l              yi=interp1(y,xi)                    默认x=1:n，n是向量y的元素个数

l              yi=interp1(x,y,xi,’method’)   指定特定算法插值，method可以是如下字符串

Ø    linear       线性插值

Ø    spline      三次样条插值

Ø    cubic      三次插值

要求：x是单调，但不要求连续等距。

如果x连续等距，可以选用快速插值法。调用函数时只需在method前加”*”，如”*spline”

 

3） 例题

①用一维线性插值解决Runge 现象

y2=interp1(x,y,x0);

plot(x0,y2,'*m')

②正弦曲线的插值示例

x=0:0.1:10;

y=sin(x);

xi=0:0.25:10;

yi=interp1(x,y,xi);

plot(x,y,’o’,xi,yi)

 

3.   Hermite插值

1) 方法介绍

已知n个插值节点x₁,x₂,…,x_n及起对应的函数值y₁,y₂,…,y_n和一阶导数值y₁’,y₂’,…,y_n’，则计算插值区域内任意x的函数值y的Hermite插值公式

其中：  , 

2) Matlab实现       Hermite.m

 

function y=hermite(x0,y0,y1,x)

n=length(x0);m=length(x);

for k=1:m

   yy=0.0;

   for i=1:n

      h=1.0;

      a=0.0;

      for j=1:n

         if j~=i

            h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2;

            a=1/(x0(i)-x0(j))+a;

         end

      end

      yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i));

   end

   y(k)=yy;

end

 

3） 例题

如下数据表，构造Hermite多项式，并求出Sin0.34 的近似值

x	0.30	0.32	0.35
Sin(x)	0.29552	0.31457	0.34290
Cos(x)	0.95534	0.94924	0.93937

x0=[0.3 0.32 0.35];

y0=[0.29552 0.31457 0.34290];

y1=[0.95534 0.94924 0.93937];

x=[0.3:0.005:0.35];

y=hermite(x0,y0,y1,x);

plot(x,y)

y=hermite(x0,y0,y1,0.34);

y

sin(0.34)

y2=sin(x);

hold on

plot(x,y2,'--r')

 

4.  三次样条插值

y3=interp1(x,y,x0,'*spline');

y3=spline(x,y,x0);

plot(x0,y3,'-g')

四、        数据拟合

1.   方法介绍

在实际生活中，往往需要从一组实验数据（x_i,y_i）中寻找出变量x,y之间的函数关系。由于观测数据不可避免出现误差，因此并不需要y=f(x)一定要经过所有的点，而只要求在给定点x_i上误差Δi=f(x_i)-y_i按某种标准达到最小。通常用欧式范数║Δ║²作为误差量度的标准。这就是所谓的最小二乘法。

 

注意：数据拟合与插值的最大区别在于拟合需要给出一个曲线方程的具体解析形式，而插值只需求出该点的内插数值。

 

1)线性拟合

线性拟合以最简单的一次线性方程f(x)=a₁x+a₀拟合数据。按最小二乘法，a₁,a₀需满足 最小，因此可以通过  求出此时的a1,a0

2)超定方程的解法

求解方程Ax=b，其中A为m*n阶矩阵，x(1:n)，b(1:m)均为列向量，且m>n。由于该超定方程个数多于未知数个数，当增广矩阵的秩大于系数矩阵的秩时无解。现在求其最小二乘解，它就是使余向量r_x=b-Ax的谱范数║r_x║₂=(r_x^Tr_x)^1/2 最小的n维向量。具体解法可以通过求解该方程组的法方程组A^TAx=A^Tb获得。

 

2.  Matlab的实现

1)线性拟合及多项式拟合

ployfit(x,y,i)   以最高次为i的多项式拟合数据点(x,y)

 

例1 

x=[0 1 2 3 4 5];

y=[0 21 62 70 77 110];

coef=polyfit(x,y,1);

a1=coef(1),a0=coef(2);

ybest=a1*x+a0;

s=sum((y-ybest).^2);

axis([-1,6,-20,120]);

plot(x,y, '*')

hold on

plot(x,ybest)

 

例2

如下给出从二阶到十阶多项式拟合曲线的比较程序，并给出拟合曲线

 

 

x=[0 1 2 3 4 5];

y=[0 21 62 70 77 110];

xi=0:0.2:5;

for n=2:10

    bb=polyfit(x,y,n);

    yi=polyval(bb,xi);

    plot(xi,yi,x,y, '* ')

    title([int2str(n), '次多项式拟合曲线'])

    grid on

    pause

end

 

例3

在某个实验中得到如下一组数据：

x	1	2	3	4	5	6	7
y	0.3101	0.4900	0.6400	0.8000	0.9200	1.0500	1.2000

已知x,y满足y=kxⁿ，求参数k与n。

 

 提示:

   y=kxⁿàlny=lnk+nlnx     

   LOG(x)

           EXP(x)

 

** 可线性化的非线性模型

模型形式	变换后形式	变量和参数的变化
		Y	X	a1	a2

 

2) 超定方程的解法

例：用最小二乘法求一个形如y=a+bx2的经验公式，是其拟合下表数据：

xi	19	25	31	38	44
yi	19.0	32.3	49.0	73.3	97.8

 

x=[19 25 31 38 44];

y=[19.0 32.3 49.0 73.3 97.8];

x1=x.^2;

x1=[ones(5,1),x1'];

ab=x1/y';

ab

x0=[19:0.2:44];

y0=ab(1)+ab(2)*x0.^2;

plot(x,y, 'o')

hold on

plot(x0,y0, '-r')

插值与拟合的区别在于是否通过所有已知点来逼近函数，而在图象放大中，假定 原图象是无噪音的即图象的每一点都是有用信息，所以大多数图形放大都采用 8 插值法。

最近邻域法，线性插值法和三次样条法是常用的几种算法 所谓 “ 最近邻域法 ” 就是把这个非整数坐标作一个四舍五入，取最近的整 数点坐标处的点的颜色。

而 “ 线性插值法 ” 就是根据周围最接近的几个点（对于 平面图像来说，共有四点）颜色作线性插值计算（对于平面图像来说就是二维 线性插值）来估计这点的颜色，大多数情况下，准确度要高于最近邻域法， 当然效果也要好得多，最明显的就是放大时，图像边缘的锯齿比最近邻域法小 非常多。

线性插值：对于一个目的像素，设置坐标通过反向变换得到浮点坐标为 i+u,j+v 其中 i j 均为非负整数， u v 为 [0,1 区间的浮点数，则这个像素 得值 f i+u,j+v 可由原图像中坐标为 ( i,j i+1,j i,j+1 i+1,j+1 所对应的 周围四个像素的值决定，即： f i+u,j+v = 1-u 1-v f i,j + 1-u vf i,j+1 +u 1-v f i+1,j +uvf i+1,j+1 其中 f i,j 表示源图像 ( i,j 处的像素值，以此类推，这就是线性插值法。

线性 插值法计算量大，但缩放后图像质量高，不会出现像素值不连续的情况。由于 线性插值具有低通滤波器的性质，使高频分量受损，所以可能会使图像轮廓在一 定程度上变得模糊。