matlab实现插值法和曲线拟合

插值法和曲线拟合
电子科技大学
摘要：理解拉格朗日多项式插值、分段线性插值、牛顿前插，曲线拟合，用matlab编程求解函数，用插值法和分段线性插值求解同一函数，比较插值余项；用牛顿前插公式计算函数，计算函数值；对于曲线拟合，用不同曲线拟合数据。

关键字：拉格朗日插值多项式；分段线性插值；牛顿前插；曲线拟合
引言：
在数学物理方程中，当给定数据是不同散点时，无法确定函数表达式，求解函数就需要很大的计算量，我们有多种方法对给定的表格函数进行求解，我们这里，利用插值法和曲线拟合对函数进行求解，进一步了解函数性质，两种方法各有利弊，适合我们进行不同的散点函数求解。

正文：
一、插值法和分段线性插值
1拉格朗日多项式原理
对某个多项式函数，已知有给定的k + 1个取值点：
其中对应着自变量的位置，而对应着函数在这个位置的取值。

假设任意两个不同的x j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：
其中每个为拉格朗日基本多项式（或称插值基函数），其表达式为：
[3]
拉格朗日基本多项式的特点是在上取值为1，在其它的点上取值为0。

2分段线性插值原理
给定区间[a,b], 将其分割成a=x
0<x
1
<…<x
n
=b, 已知函数y= f(x) 在这些插值结点的函数
值为
y k =f(x
k
)(k=0,1,…,n)求一个分段函数I
h
(x), 使其满足：
(1) I h (x k )=y k ，(k=0,1,…,n) ;
(2) 在每个区间[x k ,x k+1 ] 上,I h (x)是个一次函数。

易知,I h (x)是个折线函数, 在每个区间[x k ,x k+1 ]上,(k=0,1,…,n) k
1k k 1k 1
k k 1k k k ,1)
()
()(x x x x x f x x x x x f x L --+--=++++，
于是, I h (x)在[a,b]上是连续的，但其一阶导数是不连续的。

3拉格朗日插值多项式算法
○1输入,(0,1,2,,)i i x y i n = ，令0)(=x L n 。

○
2对0,1,2,,i n = ，计算 0,()()/()
n
i j i j j j i
l x x x x x -≠=
--∏
()()()n n i i
L x L x l x y ←−−+
4分段线性插值算法
○1输入（x k ,y k ）,k=0,1,…,n;
○2计算k
1k k 1k 1
k k 1k k k ,1)
()
()(x x x x x f x x x x x f x L --+--=++++
5插值法和分段线性插值程序
按下列数据分别作五次插值和分段线性插值，画出两条插值曲线以及给定数据点。

求x 1=0.32,
functionlagrint
xi=[0.32,0.55,0.68];
%xi=[0.2:0.001:0.8];
x=[0.3,0.42,0.50,0.58,0.66,0.72];
y=[1.04403,1.08462,1.11803,1.15603,1.19817,1.23223]; L=zeros(size(y)); m=length(xi); fori=1:m
dxi=xi(i)-x;
L(1)=prod(dxi(2:6))/prod(x(1)-x(2:6));
L(6)=prod(dxi(1:6-1))/prod(x(6)-x(1:6-1)); for j=2:6-1
num=prod(dxi(1:j-1))*prod(dxi(j+1:6));
den=prod(x(j)-x(1:j-1))*prod(x(j)-x(j+1:6));
L(j)=num/den;
end
yi(i)=sum(y.*L);
fprintf('x=%f,y=%f\n',xi(i),yi(i));
end
plot(xi,yi,'r');
axis([0.2 0.8 1.03 1.24]);
hold on
plot(x,y,'b.','markersize',20)
grid on
分段线性插值算法程序：
function [y]=div
%xi=[0.3:0.001:0.72];
x0=[0.3,0.42,0.50,0.58,0.66,0.72];
y0=[1.04403,1.08462,1.11803,1.15603,1.19817,1.23223]; k=1;
xi=[0.32,0.55,0.68];
for j=1:3
fori=1:5
if xi(j)>=x0(i) && xi(j)<=x0(i+1) && k<=3
lx(1)=(xi(j)-x0(i+1))/(x0(i)-x0(i+1));
lx(2)=(xi(j)-x0(i))/(x0(i+1)-x0(i));
y(k)=lx(1)*y0(i)+lx(2)*y0(i+1);
k=k+1;
end
end
end
plot(xi,y,'r');
axis([0.2 0.8 1.03 1.24]);
hold on
plot(x0,y0,'b.','markersize',20)
grid on
6运算结果
拉格朗日插值结果
x=0.320000,y=1.049958 x=0.550000,y=1.141271 x=0.680000,y=1.209300 拉格朗日插值余项：
)
72.0)(66.0)(58.0)(5.0)(42.0)(3.0(!
6)()
()()(6
55------=-=x x x x x x f x L x f x R ξ
分段插值结果
ans =
1.0508 1.1418 1.2095 分段线性插值余项：
)
(max
8
)
)((2
)(max
)()(max
1
1
12
1,1x f h x x x x f x L x f R i i i i i i x x x i
i i i x x x i x x x ''≤--''≤-=+++≤≤+≤≤≤≤ξ
由于拉格朗日插值的余项比分段线性插值的余项要求更为严格，点少、区间小的时候，拉格朗日插值要更好。

但在区间较大、节点较多的时候，分段线性插值要更好。

二、牛顿前插 1牛顿前插原理
n 次牛顿前插公式：
)
1()1(!
)1(!
2!
1)(!
)(002
001
00+--∆+
+-∆+
∆+
=-∏∆=
+-==∑
n t t t n f t t f t f f j t k f th x N n
k j n
k k
n 插值余项：
)
()!
1()
()1()()
1(1
0ξ+++--=
+n n n f
h
n n t t t th x R ，
)
,(0n x x ∈ξ
m 阶差分记作i m i m i m
f f f 1
11
-+-∆
-∆=∆。

[]i f x k 阶差商是[]∑
∏=≠=-=
k
m k
m
i i i m
m k x x
x f x x x f 0
010)()
(,,,
差分和差商之间的关系是[]k
i k
k i i i h
k f x x x f !,,,1∆=++
2牛顿前插算法
○1输入),2,1,0(,,n i y x n i
i =。

○
2对n k ,,3,2,1 =，k i ,,2,1 =计算各阶差分i m i m i m f f f 111-+-∆-∆=∆ ○
3计算函数值)
1()1(!
)1(!
2!
1)(!
)(002
001
00+--∆+
+-∆+
∆+
=-∏∆=
+-==∑
n t t t n f t t f t f f j t k f th x N n
k j n
k k
n 3牛顿前插程序：
编写一个用牛顿前插公式计算函数值的程序，要求先输出差分表，再计算x 点的函数值
分别求x =0.158和x =0.636的三次插值的值，并比较二者的插值余项。

这里以x=0.636为例 function [P]=newtoncha x0=0.636;
X=[0.125 0.250 0.375 0.500 0.625 0.750]; Y=[0.796 0.773 0.744 0.704 0.656 0.602]; h=abs(X(2)-X(1)); n=find(abs(x0-X)<3*h); X=X(n(1):n(end)); Y=Y(n(1):n(end));
w=length(X);
R=zeros(w,w);
R(:,1)=Y(:);
for k=2:w
for j=k:w
R(j,k)=R(j,k-1)-R(j-1,k-1);
end
end
t=(x0-X(1))/h;
T=1;
for m=1:w-1
T=T*(t-m+1);
N(m)=R(m+1,m+1)*T/factorial(m); end
P=R(1,1)+sum(N);
4运行结果：
差分表
X=0.636时ans =0.651459661824000 x=0,158时ans =0.790229818880000 三、曲线拟合 1曲线拟合原理： 给定数据
(,),1,2,,j j x y j n
= 。

记拟合函数的形式为
0011()()()()
m m p x a x a x a x ϕϕϕ=+++ （1.1），
其中{}0()m
k
k x ϕ=为已知的线性无关函数。

求系数
*
*
*
01,,,m
a a a 使得
2
2
011
1
(,,,)[()][()]
n
n
m
m j
j k
k
j j j j k a a a p x
y a x y ϕϕ
====
-=
-∑∑∑ （1.2）
取最小值。

称
*
*0
()()
m
k
k k p x a
x ϕ==
∑（1.3）为拟合函数或经验公式。

如果
)
,,1,0()(m k x x k
k ==ϕ，则（1.3）为m 次最小二乘拟合多项式
2曲线拟合算法：
已知数据对
)
,,2,1)(,(n j y x i i =，求多项式
)
()(0
n m x
a x P m
i i
i
<=
∑=，使得
∑∑==-=
Φn
j j m
i i j
i
m y x
a a a a 1
2
010)
(),,,( 为最小。

注意到此时k
k x
x =)(ϕ，
多项式系数m
a a a ,,,10 满足下面的线性方程组： ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡++m m m m m
m m T T T a a a S S S S S S S S S
101021
1
2110
其中
)
2,,2,1,0(1
m k x
S n
j k j
k ==
∑=，)
,,2,1,0(1
m k x y T n
j k
j
j k ==
∑
=
然后只要调用线性方程组的函数程序即可 3曲线拟合程序：
试分别用抛物线y =a +bx 2
和指数曲线y =ae bx
拟合下列数据
画出数据点和两条拟合曲线，并通过计算2个拟合函数残差向量的2范数来比较拟合优劣。

用抛物线y=a+bx 拟合程序： function ZXE
x=[1 2.5^2 3.5^2 4^2];
y=[3.8 1.50 26.0 33.0];
m=1;
S=zeros(1,2*m+1);T=zeros(m+1,1); for k=1:2*m+1
S(k)=sum(x.^(k-1));
end
for k=1:m+1
T(k)=sum(x.^(k-1).*y);
end
A=zeros(m+1,m+1);a=zeros(m+1,1); fori=1:m+1
for j=1:m+1
A(i,j)=S(i+j-1);
end
end
a=A\T;
for k=1:m+1
fprintf('a[%d]=%f\n',k,a(k)); end
p=polyfit(x,y,1);
u=polyval(p,x);
plot(sqrt(x),u,'b')
hold on
plot(sqrt(x),y,'b.')
grid on
指数曲线y=aebx拟合程序：function ZXE2
x=[1 2.5 3.5 4];
y=[3.8 1.50 26.0 33.0];
y=log(y);
m=1;
S=zeros(1,2*m+1);T=zeros(m+1,1); for k=1:2*m+1
S(k)=sum(x.^(k-1));
end
for k=1:m+1
T(k)=sum(x.^(k-1).*y);
end
A=zeros(m+1,m+1);a=zeros(m+1,1); fori=1:m+1
for j=1:m+1
A(i,j)=S(i+j-1);
end
end
a=A\T;
for k=1:m+1
fprintf('a[%d]=%f\n',k,a(k)); end
p=polyfit(x,y,1);
u=polyval(p,x);
plot(x,exp(u),'r')
hold on
plot(x,exp(y),'b.')
grid on
4运行结果：
拟合曲线：
a[1]=-0.168731
a[2]=0.833636
结论：
求解散点函数的时候，点少、区间小的时候，拉格朗日插值要更好。

但在区间较大，节点较多的时候，分段线性插值要更好。

对函数的求解和线性拟合如上。

参考文献：
[1] 孙志忠等《计算方法与实习》第五版东南大学出版社。