数值分析课程设计三次样条插值Word文件下载.docx
《数值分析课程设计三次样条插值Word文件下载.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计三次样条插值Word文件下载.docx（20页珍藏版）》请在冰豆网上搜索。
总成绩：
教研室主任签名：
年月日
课程设计任务书
2008—2009学年第二学期
专业班级：
09普本信计1班学号：
姓名：
课程设计名称：
数值分析
设计题目：
三次样条插值
完成期限：
自2012年6月8日至2012年6月13日共1周
设计依据、要求及主要内容：
一、设计目的
熟练掌握三次样条插值算法的原理和推导过程，并且能够应用Matlab软件编写相应的程序和使用Matlab软件函数库软件。
二、设计要求
（1）用Matlab函数库中相应函数对选定的问题，求出具有一定精度的结果。
（2）使用所用的方法编写Matlab程序求解,对数值结果进行分析。
（3）对于使用多个方法解同一问题的，在界面上设计成菜单形式。
三、设计内容
首先构造三次样条插值函数的定义和一般特征,并对实例问题进行实例分析,并总结
四、参考文献
[1]黄明游,冯果忱.数值分析[M].北京:
高等教育出版社,2008.
[2]马东升,雷勇军.数值计算方法[M].北京:
机械工业出版社,2006.[3]石博强,赵金.MATLAB数学计算与工程分析范例教程[M].北京:
中国铁道出版社.2005.
[4]郝红伟，MATLAB6，北京，中国电力出版社，2001
[5]姜健飞，胡良剑，数值分析及其MATLAB实验，科学出版社，2004
[6]薛毅，数值分析实验，北京工业大学出版社，2005
计划答辩时间：
2012年6月18日
指导教师（签字）：
教研室主任（签字）：
批准日期：
年月
摘要
分段低次样条插值虽然计算简单、稳定性好、收敛性有保证且易在电子计算机上实现，但只能保证各小段曲线在连接处的连续性，不能保证整件曲线的光滑性。
利用样条插值，既可保持分段低次插值多项式，又可提高插值函数光滑性。
故给出分段三次样条插值的构造过程、算法步骤，利用MATLAB软件编写三次样条插值函数通用程序，并通过数值算例证明程序的正确性。
关键字:
三次样条插值函数MATLAB编程收敛性算法步骤
一三次样条函数定义及特征
定义1：
若函数,且在每个小区间上上是三次多项式，其中是给定节点，则称是节点上的三次样条函数。
若节点上给定函数值,且
（1.1）
成立，则称为三次样条差值函数。
从定义知，要求出，在每个应小区间上确定4个待定系数，共有n个小区间，故应确定4n个参数，根据在上二阶导数连续，在节点处应满足连续性条件
（1.2）
共有3n-3个条件，再加上满足插值条件（1.1），共有4n-2个条件，因此还需要2个条件才能确定。
通常可在区间端点上各加一个条件（称边界条件），边界条件可根据实际的问题要求给定。
常见的三种：
（1）已知两端的一节导数值，即
（1.3）
（2）两端的二阶导数已知，即
（1.4）
特殊情况下的边界条件
（1.4）’
称为自然边界条件
（3）当是以为周期函数时，则要求也是周期函数，这时边界条件应满足
而此时式中,这样确定的样条函数称为周期函数。
二函数推导原理及构造
我们采用待定一阶导数的方法即设S（Xj）=Mj,j=0,1,...,n,因为分段三次Hermite插值多项式已经至少是一阶连续可导了，为了让它成为三次样条函数只需确定节点处的一阶导数使这些节点处的二阶导数连续即可！
由于在内部节点处二阶导数连续条件：
整理化简后得：
第一类三次样条插值问题方程组由于已知：
基本方程组化为n-1阶方程组
化为矩阵形式
\
这是一个严格对角占优的三对角方程组,用追赶法可以求解！
第二类三次样条插值问题的方程组,由于已知：
故得：
稍加整理得
联合基本方程组得一个n+1阶三对角方程组，化成矩阵形式为：
仍然是严格对角占优
第三类样条插值问题的方程组，由于：
立即可得下式：
其中:
联合基本方程得一个广义三对角或周期三对角方程组：
求解这些不同类型的样条插值问题的方程组，我们可得所要待定的一阶导数：
再代入S（x）的每一段表达式，就求得三次样条函数的表达式!
利用插值（即求过已知有限个数据点的近似函数）的基本原理，用多项式作为研究插值的工具，进行代数插值。
其基本问题是：
已知函数f（x）在区间[a,b]上n+1个不同点x0，…，xn处的函数值（i=0,1,…,n），求一个至多n次多项式ψn（x）使其在给定点处与f（x）同值，即满足插值条件:
ψn（x）==.许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续，而且要有连续的曲率，这就导致了样条插值的产生。
数学上将具有一定光滑性的分段多项式称为样条函数。
具体地说，给定区间[a,b]的一个分划Δ如果函数s（x）满足：
（i）在每个小区间[]（i=0,1,…,n）上s（x）是k次多项式；
（ii）s（x）在[a,b]上具有k−1阶连续导数。
则称s（x）为关于分划Δ的k次样条函数，其图形称为k次样条曲线。
由于三次样条插值我、函数s（x）的插值节点处的二阶导数存在，因此令各节点处的二阶导数为
（1.01）
根据样条插值函数的定义，三次样条插值函数是s（x）在每一个小区间上市不超过三次的多项式。
在每一个小区间上，其二阶导数为线性函数，即
（1.02）
对式（1.02）积分两次，则得到
（1.03）
其中为任意常数。
又根据样条插值函数定义中的条件（3），即
可以确定与为
=
（1.04）
将式（1.04）中与的值代入表达式（1.03后，就可以得到样条插值函数在区间上的表达式为
（1.05）
其中与分别为区间两端点处的二阶导数值。
由此可以看处，只要能确定各点处的二阶导数值，则子渠道间上的三次样条插值函数也确定了。
在区间[a,b]上的一阶导数连续，在各节点的左右两子区间上的s（x）虽然不同，但在连接点处的导数存在，即在连接点处的左右导数相等，有
（1.06）
为了利用条件（2.18），在x属于时，县求为
（1.07）
当x属于时，
（1.08）
整理得：
（1.09）
其中
三问题的提出
上面讨论的分段低次插值函数都有一致收敛性，但光滑性较差，对于像高速飞机的机翼形线，船体放样等型值线往往要求有二阶光滑度，即有二阶连续导数，早期工程师制图时，把富有弹性的细长木条（所谓样条）用压铁固定在样点上，在其他地方让它自由弯曲，然后画下长条的曲线，称为样条曲线。
它实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续，从数学上加以概括就得到数学样条这一概念。
下面我们讨论最常用的三次样条函数。
四实例应分析函数的MATLAB的程序设计
例1：
已知一组数据点，编写一程序求解三次样条插值函数满足
并针对下面一组具体实验数据
0.25
0.3
0.39
0.45
0.53
0.5000
0.5477
0.6245
0.6708
0.7280
求解，其中边界条件为.
1）三次样条插值自然边界条件源程序：
functions=spline3（x,y,dy1,dyn）
%x为节点，y为节点函数值，dy1，dyn分别为x=0.25,0.53处的二阶导
m=length（x）;
n=length（y）;
ifm~=n
error（'
xory输入有误'
）
return
end
h=zeros（1,n-1）;
h（n-1）=x（n）-x（n-1）;
fork=1:
n-2
h（k）=x（k+1）-x（k）;
v（k）=h（k+1）/（h（k+1）+h（k））;
u（k）=1-v（k）;
g
（1）=3*（y
（2）-y
（1））/h
（1）-h
（1）/2*dy1;
g（n）=3*（y（n）-y（n-1））/h（n-1）+h（n-1）/2*dyn;
fori=2:
n-1
g（i）=3*（u（i-1）*（y（i+1）-y（i））/h（i）+v（i-1）*（y（i）-y（i-1））/h（i-1））;
n-1;
A（i,i-1）=v（i-1）;
A（i,i+1）=u（i-1）;
A（n,n-1）=1;
A（1,2）=1;
A=A+2*eye（n）;
M=zhuigf（A,g）;
%调用函数，追赶法求M
fprintf（'
三次样条（三对角）插值的函数表达式\n'
）;
symsX;
fprintf（'
S%d--%d:
\n'
k,k+1）;
s（k）=（h（k）+2*（X-x（k）））./h（k）.^3.*（X-x（k+1））.^2.*y（k）...
+（h（k）-2*（X-x（k+1）））./h（k）.^3.*（X-x（k））.^2.*y（k+1）...
+（X-x（k））.*（X-x（k+1））.^2./h（k）.^2*M（k）+（X-x（k+1））.*...
（X-x（k））.^2./h（k）.^2*M（k+1）;
s=s.'
;
s=vpa（s,4）;
%画三次样条插值函数图像
fori=1:
X=x（i）:
0.01:
x（i+1）;
st=（h（i）+2*（X-x（i）））./（h（i）^3）.*（X-x（i+1））.^2.*y（i）...
+（h（i）-2.*（X-x（i+1）））./（h（i）^3）.*（X-x（i））.^2.*y（i+1）...
+（X-x（i））.*（X-x（i+1））.^2./h（i）^2*M（i）+（X-x（i+1））.*...
（X-x（i））.^2./h（i）^2*M（i+1）;
plot（x,y,'
o'
X,st）;
holdon
End
plot（x,y）;
gridon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%调用的函数：
%追赶法
functionM=zhuigf（A,g）
n=length（A）;
L=eye（n）;
U=zeros（n）;
U（i,i+1）=A（i,i+1）;
U（1,1）=A（1,1）;
n
L（i,i-1）=A（i,i-1）/U（i-1,i-1）;
U（i,i）=A（i,i）-L（i,i-1）*A（i-1,i）;
Y
（1）=g
（1）;
Y（i）=g（i）-L（i,i-1）*Y（i-1）;
M（n）=Y（n）/U（n
升级会员 