java的拟合与插值库,插值和拟合区别

插值和拟合都是函数逼近或者数值逼近的重要组成部分。他们的共同点都是通过已知一些离散点集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]上有互异点x0,x1,…,xn处取值y0,y1,…,yn

。如果函数φ(x)在点xi上满足φ(xi)=yi

(i=0,1,2,…,n)，则称φ(x)是函数y=f(x)的插值函数，x0,x1,…,xn是插值节点。若此时φ(x)是代数多项式P(x)，则称P(x)为插值多项式。显然

f(x)≈φ(x)，x∈[a,b]

1. 拉格朗日插值

构造n次多项式Pn

(x)= yk lk (x)=y0l0 (x)+y1l1

(x)+…+ynln (x)，这是不超过n次的多项式，其中基函数lk(x)=

显然lk (x)满足lk (xi)=

此时 Pn(x)≈f(x),误差Rn(x)=f(x)-Pn(x)=

其中 ∈(a,b)且依赖于x，

=(x-x0)(x-x1)…(x-xn)

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

P1(x)=yklk(x)+yk+1lk+1(x)

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

2. 牛顿插值

构造n次多项式Nn(x)=f(x0)+f(x0,x1)(x-x0)+f(x0,x1,x2)(x-x0)(x-x1)+…

+f(x0,x1,x2,…,xn)(x-x0)(x-x1)…(x-xn)

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

(二个节点，一阶差商)

(三个节点，二阶差商)

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

注意：由于插值多项式的唯一性，有时为了避免拉格朗日余项Rn(x)中n+1阶导数的运算，用牛顿插值公式Rn

(x)=f(x)-Nn(x)=f(x,x0,…,xn)ωn+1(x),

其中ωn+1(x)=(x-x0)(x-x1)…(x-xn)

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

1) 线性插值

已知n+1个不同节点x0,x1,…,xn

,构造分段一次线性多项式P(x)，使之满足

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

l P(xk)=yk

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

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

区间[a,b]上两个互异节点xi,xi+1，已知实数y i,y i+1,m

i,m i+1，为了构造次数不大于3的多项式 满足条件

引入 , 使之满足

可以求出

此时 = + ，其中

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

对于给定n+1个不同节点x0,x1,…,xn及函数值y0,y1,…,yn，其中a=x0

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

l S(xk)=yk (k=0,1,…,n)

l 每个子区间[xk,xk+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)拉格朗日插值

P3(x)=y0l0(x)+y1l1(x)+y2l2(x)+y3l3(x)

=(-5)

+ (-6)  +3

=x3-4x2+3

φ(x)≈P3(x)= x3-4x2+3 φ(2.5)=2.53-4*2.52+3=-6.375

(2)牛顿插值

xi

f(xi)

一阶差商

二阶差商

三阶差商

1

0

2

–5

-5

3

–6

-1

2

4

3

9

5

1

N3(x)=f(x0)+f(x0,x1)(x-x0)+f(x0,x1,x2)(x-x0)(x-x1)+f(x0,x1,x2,x3)(x-x0)(x-x1)(x-x2)

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

=x3-4x2+3

三、 Matlab在插值中的应用

1. Lagrange插值

１)方法回顾

对给定的n个节点x1,x2,…,xn及对应的函数值y1,y2,…,yn，利用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个等距结点作插值多项式Pn(x)，使其在个结点的值与函数y(x)=1/(1+25x2)在结点的值相等。但在n→∞时，插值多项式Pn(x)在区间中部趋于y(x)。但对于0.726≤∣x∣≤1的x，Pn(x)严重发散。

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

)

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个插值节点x1,x2,…,xn及起对应的函数值y1,y2,…,yn和一阶导数值y1’,y2’,…,yn’，则计算插值区域内任意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. 方法介绍

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

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

1)线性拟合

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

求出此时的a1,a0

2)超定方程的解法

求解方程Ax=b，其中A为m*n阶矩阵，x(1:n)，b(1:m)均为列向量，且m>n。由于该超定方程个数多于未知数个数，当增广矩阵的秩大于系数矩阵的秩时无解。现在求其最小二乘解，它就是使余向量rx=b-Ax的谱范数║rx║2=(rxTrx)1/2

最小的n维向量。具体解法可以通过求解该方程组的法方程组ATAx=ATb获得。

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=kxn，求参数k与n。

提示:

y=kxnàlny=lnk+nlnx

LOG(x)

EXP(x)

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

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')