三次样条拟合范例.docx
《三次样条拟合范例.docx》由会员分享，可在线阅读，更多相关《三次样条拟合范例.docx（10页珍藏版）》请在冰豆网上搜索。
三次样条拟合范例
1设计目的、要求
对龙格函数
在区间[-1,1]上取
的等距节点，分别作多项式插值、三次样条插值和三次曲线拟合，画出
及各逼近函数的图形，比较各结果。
2设计原理
（1）多项式插值：
利用拉格朗日多项式插值的方法，其主要原理是拉格朗日多项式，即：
表示待插值函数的
个节点，
，其中
；
（2）三次样条插值：
三次样条插值有三种方法，在本例中，我们选择第一边界条件下的样条插值，即两端一阶导数已知的插值方法：
（3）三次曲线拟合：
本题中采用最小二乘法的三次多项式拟合。
最小二乘拟合是利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。
在本题中，n=10，故有11个点，以这11个点的
和
值为已知数据，进行三次多项式拟合，设该多项式为
，该拟合曲线只需
的值最小即可。
3采用软件、设备
计算机、matlab软件
4设计内容
1、多项式插值：
在区间
上取
的等距节点，带入拉格朗日插值多项式中，求出各个节点的插值，并利用matlab软件建立m函数，画出其图形。
在matlab中建立一个lagrange.m文件，里面代码如下：
%lagrange函数
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
建立一个polynomial.m文件，用于多项式插值的实现，代码如下：
%lagrange插值
x=[-1:
0.2:
1];
y=1./（1+25*x.^2）;
x0=[-1:
0.02:
1];
y0=lagrange（x,y,x0）;
y1=1./（1+25*x0.^2）;
plot（x0,y0,'--r'）
%插值曲线
holdon
%原曲线
plot（x0,y1,'-b'）
运行duoxiangshi.m文件，得到如下图形：
2、三次样条插值：
所谓三次样条插值多项式
是一种分段函数，它在节点
分成的每个小区间
上是3次多项式，其在此区间上的表达式如下：
因此，只要确定了
的值，就确定了整个表达式，
的计算方法如下：
令：
则
满足如下
个方程：
对于第一种边界条件下有
如果令
那么解就可以为
求函数的二阶导数：
>>symsx
>>f=sym（1/（1+25*x^2））
f=
1/（1+25*x^2）
>>diff（f）
ans=
-（50*x）/（25*x^2+1）^2
将函数的两个端点，代入上面的式子中：
f’（-1）=0.0740
f’
（1）=-0.0740
求出从-1到1的n=10的等距节点，对应的x，y值
对应m文件代码如下：
forx=-1:
0.2:
1
y=1/（1+25*x^2）
end
y=
得出
x=-1-0.8-0.6-0.4-0.200.20.40.60.81
y=0.03850.05880.10.20.510.50.20.10.05880.0385
编写m文件Three_Spline.m
x=linspace（-1,1,11）;
y=1./（1+25*x.^2）;
[m,p]=scyt1（x,y,0.0740,-0.0740）;
holdon
x0=-1:
0.01:
1;
y0=1./（1+25*x0.^2）;
plot（x0,y0,'--b'）
得到如下图像：
.
其中蓝色曲线为原图，红色曲线为拟合后的图像。
3、三次曲线拟合：
这里我们使用最小二乘法的3次拟合
建立一个Three_fitting.m文件，代码如下：
%主要代码
x=[-1-0.8-0.6-0.4-0.200.20.40.60.81];
y=[0.03850.05880.10.20.510.50.20.10.05880.0385];
a=polyfit（x,y,3）;
x1=[-1:
0.01:
1];
y1=a（4）+a（3）*x1+a
（2）*x1.^2+a
（1）*x1.^3;
x0=[-1:
0.01:
1];
y0=1./（1+25*x0.^2）
%原曲线
plot（x0,y0,'-r'）
holdon
%三次拟合曲线
plot（x1,y1,'-b'）
上图中，蓝色部分为三次拟合曲线，红色部分为原曲线
6结果分析
拉格朗日插值的优点是对于某一区域，不限于被估计点周围，公式简单易实施。
一般认为n的次数越高，逼近
的精度就越好，但在本题中，对龙格函数
，中间部分插值效果比较好，而对于两端，插值结果是非常不理想的，即龙格现象。
样条函数可以给出光滑的插值曲线，从本题中就能体现出来。
从以上图形可以看出，三次样条插值的图形是比较逼近于原图的，收敛性相对而言是非常好的，但在本题中，仅将原区间分成10个等距区间，因此，逼近效果还不是特别理想，当我们将n增大时，插值后的曲线越逼近于原曲线。
总的来说，三次样条插值的稳定性比较好，收敛性比较强。
在这三种方法中，三次曲线拟合的效果是最差的，所得的图形与原曲线差距甚远。
最小二乘法中，并不要求拟合后的曲线经过所有已知的点，只需要拟合多项式上的点在某种标准上与定点之间的差距最小即可，因此与原曲线的逼近程度是最差的。
最小二乘法的多项式拟合只适用于多项式，而本题中的函数并不是一个多项式，因此，不建议使用最小二乘法拟合。
参考文献:
[1]李庆扬王能超等.数值分析[M].清华大学出版社
[2]吴振远.科学计算实验指导书基于MATLAB数值分析[M].中国地质大学出版社
[3]宋叶志.MATLAB数值分析与应用[M].机械工业出版社,2009.07
附录
三次样条插值主要代码：
function[m,p]=scyt1（x,y,df0,dfn）
n=length（x）;
r=ones（n-1,1）;
u=ones（n-1,1）;
d=ones（n,1）;
r
（1）=1;
d
（1）=6*（（y
（2）-y
（1））/（x
（2）-x
（1））-df0）/（x
（2）-x
（1））;
u（n-1）=1;
d（n）=6*（dfn-（y（n）-y（n-1））/（x（n）-x（n-1）））/（x（n）-x（n-1））;
fork=2:
n-1
u（k-1）=（x（k）-x（k-1））/（x（k+1）-x（k-1））;r（k）=（x（k+1）-x（k））/（x（k+1）-x（k-1））;
d（k）=6*（（y（k+1）-y（k））/（x（k+1）-x（k））-（y（k）-y（k-1））/（x（k）-x（k-1）））/（x（k+1）-x（k-1））;
end
A=eye（n,n）*2;
fork=1:
n-1
A（k,k+1）=r（k）;
A（k+1,k）=u（k）;
end
m=A\d;
ft=d
（1）;
symst
fork=1:
n-1%求s（x）即插值多项式
p（k,1）=m（k）/（6*（x（k+1）-x（k）））;
p（k,2）=m（k+1）/（6*（x（k+1）-x（k）））;
p（k,3）=（y（k）-m（k）*（x（k+1）-x（k））^2/6）/（x（k+1）-x（k））;
p（k,4）=（y（k+1）-m（k+1）*（x（k+1）-x（k））^2/6）/（x（k+1）-x（k））;
sx（k）=p（k,1）*（x（k+1）-t）^3+p（k,2）*（t-x（k））^3+p（k,3）*（x（k+1）-t）+p（k,4）*（t-x（k））;
end
kmax=1000;
xt=linspace（x
（1）,x（n）,kmax）;
fori=1:
n-1%出点xt对应的y值
fork=1:
kmax
ifx（i）<=xt（k）&xt（k）<=x（i+1）
fx（k）=subs（sx（i）,xt（k））;
end
end
end
plot（xt,fx,'r'）;xlabel（'x'）;
ylabel（'y'）;title（'f'）;
text（x（fix（n/2））,y（fix（n/2））,'f'）
holdon
plot（x,y,'*'）
holdoff
升级会员 