插值与拟合.docx
《插值与拟合.docx》由会员分享,可在线阅读,更多相关《插值与拟合.docx(13页珍藏版)》请在冰豆网上搜索。
插值与拟合
插值与拟合
实验内容:
1.三种插值方法
2.用Matlab计算插值
3.拟合的基本原理
4.用Matlab拟合曲线
实验目的:
掌握插值与拟合方法
一、概念的引入
1.插值与拟合在现实生活中的应用
●机械制造:
汽车外观设计
●采样数据的重新建构:
电脑游戏中场景的显示,地质勘探,医学领域(CT)
2.概念的定义
●插值:
基于[a,b]区间上的n个互异点,给定函数f(x),寻找某个函数去逼近f(x)。
若要求φ(x)在xi处与f(xi)相等,这类的函数逼近问题称为插值问题,xi即是插值点
●逼近:
当取值点过多时,构造通过所有点的难度非常大。
此时选择一个次数较低的函数最佳逼近这些点,一般采用最小二乘法
●光顾:
曲线的拐点不能太多,条件:
二阶几何连续
不存在多余拐点
曲率变化较小
●拟合:
曲线设计过程中用插值或通过逼近方法是生成的曲线光滑(切变量连续)光顾
二、插值理论
设函数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)=yklk(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),使之满足
●P(x)在[a,b]上连续
●P(xk)=yk
●P(x)在[xi,xi+1]上是线性函数,P(x)=
2)两点带导数插值---避免尖点、一阶连续
区间[a,b]上两个互异节点xi,xi+1,已知实数yi,yi+1,mi,mi+1,为了构造次数不大于3的多项式
满足条件
引入
使之满足
可以求出
此时
=
+
,其中
4.三次样条插值------二阶可导
对于给定n+1个不同节点x0,x1,…,xn及函数值y0,y1,…,yn,其中a=x0构造三次样条插值函数S(x)。
S(x)称为三次样条函数时需满足:
●S(x)在[a,b]上二阶导数连续
●S(xk)=yk(k=0,1,…,n)
●每个子区间[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插值
1)方法回顾
对给定的n个节点x1,x2,…,xn及对应的函数值y1,y2,…,yn,利用n次Lagrange插值多项式公式,插值区间内任意x的函数值y可以通过下式求出:
2)Matlab实现函数Lagrange.m
functiony=lagrange(x0,y0,x)
n=length(x0);m=length(x);
fori=1:
m
z=x(i);
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
ifj~=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=[1234];
y=[0–5–63];
lagrange(x,y,2.5)
2.Runge现象及分段线性插值
1)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')
holdon
plot(x0,y1,'-b')
2)Matlab实现分段插值------一维插值interp1
●yi=interp1(x,y,xi)对(x,y)进行插值,计算插值点xi的函数值
●yi=interp1(y,xi)默认x=1:
n,n是向量y的元素个数
●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
functiony=hermite(x0,y0,y1,x)
n=length(x0);m=length(x);
fork=1:
m
yy=0.0;
fori=1:
n
h=1.0;
a=0.0;
forj=1:
n
ifj~=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.30.320.35];
y0=[0.295520.314570.34290];
y1=[0.955340.949240.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);
holdon
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=[012345];
y=[021627077110];
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,'*')
holdon
plot(x,ybest)
例2
如下给出从二阶到十阶多项式拟合曲线的比较程序,并给出拟合曲线