沈阳航空航天大学
数学软件课程设计
（设计程序）
题目三次样条插值函数
班级 / 学号
学生姓名
指导教师
沈阳航空航天大学
课程设计任务书
课程名称数学软件课程设计
院（系）理学院专业信息与计算科学
班级学号姓名
课程设计题目三次样条插值函数
课程设计时间: 2010 年12月20日至2010 年12月31日
课程设计的内容及要求：
1.三次样条插值函数
给出函数在互异点处的值分别为。
(1)掌握求三次样条插值函数的基本原理；
(2)编写程序求在第一边界条件下函数的三次样条插值函数；
(3)在区间上取n=10，20，分别用等距节点对函数
作三次样条插值函数，利用(1)的结果画出插值函数的图形，并在该图形界面中同时画出的图形。
[要求]
1.学习态度要认真，要积极参与课程设计，锻炼独立思考能力；
2.严格遵守上机时间安排；
3.按照MATLAB编程训练的任务要求来编写程序；
4.根据任务书来完成课程设计论文；
5.报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”；
6.报告上交时间：课程设计结束时上交报告；
7.严谨抄袭行为。
指导教师年月日负责教师年月日学生签字年月日
沈阳航空航天大学
课程设计成绩评定单
课程名称数学软件课程设计
院（系）理学院专业信息与计算科学课程设计题目三次样条插值函数
学号姓名
指导教师评语：
课程设计成绩
指导教师签字
年月日
目录
一正文 (1)
1问题分析 (1)
1.1 题目 (1)
1.2 分析 (1)
2 研究方法原理 (1)
2.1 求三次样条插值多项式，算法组织 (1)
3 算例结果 (3)
二总结 (7)
参考文献 (8)
附录 (9)
源程序: (9)
程序1 (9)
程序2 (10)
程序3 (12)
程序 4 (12)
一 正文
1问题分析 1.1 题目
三次样条插值函数
给出函数在互异点处的值分别为
。
（1）掌握求三次样条插值函数的基本原理；
（2）编写程序求在第一边界条件下函数的三次样条插值函数； （3）在区间
上取n=10，20，分别用等距节点对函数
作三次样条插值函数，利用(1)的结果画出插值函数图形，并在该图形界面中同时画出
的图形。
1.2 分析
一般认为插值次数n 越高，)(x f 的精度就越高，但实际并非如此，20世纪初龙格(Runge)就发现了这一现象，因此就提出了分段低次插值分段线性插值有一致收敛性，但光滑性差，而三次样条插值具有二介光滑度，三次样条插值首先要给定n 个点和对应的函数值，还要给出边界条件如第一边界条件)(')('),0(')0('xn f xn s x f x s ==，第二边界条件)('')(''),0('')0(''xn f xn s x f x s ==，而题目要求是在给定第一边界条件下的三次样条插值。
2 研究方法原理
2.1 求三次样条插值多项式，算法组织
所谓三次样条插值多项式)(x S n 是一种分段函数，它在节点
i x 011()n n a x x x x b -=<
2233
1111111()[()()]()()666[,]1,2,,.
i i i i i i i i i i i i i i i
i i h x x h x x S x x x M x x M y M y M h h h x x x i n --------=-+-+-+-∈=???，
， 因此，只要确定了i M 的值，就确定了整个表达式，i M 的计算方法如下： 令：
111
11111116()6(,,)i i i i i i i i i i i i i i i i i i i i i h h h h h h y y y y d f x x x h h h h μλμ++++--+++?
===-?++??
--?=-=?+?
， 则i M 满足如下n-1个方程：
1121,2,,1i i i i i i M M M d i n μλ-+++==???-，
对于第一种边界条件下有
???
?
???
-=
+-=+---00011
0111
)'],([62]),['(62h f x x M M h x x f f M M n n n n n n
如果令,])
,['(6,1,)'],[(6,11
1000100---==-=
=n n n n n n h x x f f d h f x x f d μλ那么解就可以为
??
?????
? ??=???????? ?????????? ?
?----n n n n n
n n d d d d M M M M 11011011
1
10
22
22
μλμλμλ
3 算例结果
s(x)可以表示为：
其中p 为4?n 的矩阵。
当把区间5等分时，输入如下：
图 1
矩阵p 输出如下
图 2
图形如下：
图 3 5等分图像
1);
-n 0,1(k , 1)]x(k [x(k),t x(k)),-(t *p(k,4)t)-1)(x(k *p(k,3)x(k))^3-(t *p(k,2)t)^3-1)(x(k *p(k,1)sx(k) =+∈+++++=
其在不同的区间的函数可以表示如下：
当n=10即区间10等分时，输入如下：
图 4 得到的矩阵p 如下：
图 5
图形如下：
图 6 10等分图像
?????
????∈∈∈∈--∈=[0.6,1]
x 0.6),-(x *0.2087+1)-(x *0.0549+0.6)^3-(x *0.7032-x)^3-(1*1.9057[0.2,0.6]x 0.2),-(x *0.0549-x)-(0.6*1.511+0.2)^3-(x *1.9057+0.6)^3-(x *1.6311[-0.2,0.2]
x 0.2),+(x *1.511+x)-(0.2*1.511+0.2)^3+(x *1.6311-0.2)^3-(x *1.6311][-0.6,-0.2x 0.6),+(x *1.511+0.2)+(x *0.0549+0.6)^3+(x *1.6311-0.2)^3+(x *1.9057-]6.0,1[x 1),+(x *0.0549-0.6)+(x *0.2087-1)^3+(x *1.9057+0.6)^3+(x *0.7032)(x
s
当把区间20等分时，输入如下：
图7
得到的矩阵p如下：
图8
得出的图像如下：
图9 20等分图像
运行gtu.m，输入如下：
图10
运行结果为图像，图形如下：
图11 13、20等分和原图的图像
二总结
拿到题目时，我首先先弄清三次样条插值函数的基本原理，因为只有这样，在以后的编程中才会更懂的如何编写程序，才会更不会混淆题目的目的。而且，不能为了做题而做题，在做题时还要应用到其它知识。
三次样条插值在不懂的边界条件出来的结果也不一样，表达形式也不一样，特别是计算过程，条件不一样，编写的程序也自然不同，题目给出的是第一边界条件，在第一边界条件中，首先要给出始末两点的导数值，而且要给出n个点的x值和y的对应的值，题目要求，所以有在程序输入x,y和始末两点的导数值。
这课设过程中，我学到了很多，知道了自己的很多不足，如对于某些函数不能巧妙的应用，编写的程序很粗糙，这次课设，我们通过自己的思考与实践，终于完成了。刚开始我对于三次样条插值很不了解，现在，我已经对于三次样条插值有了一定的了解，特别是其中的原理，但公式还是没背下来，不过，原理最重要，这次课设，更好的让我们掌握了Matlab和数值分析，我了解到理论与实践是分不开的，为了更好的掌握一个知识，我们必须通过不断的实践。要学好一门计算机语言，既要掌握其最基本的语言结构，而且要特别的熟悉，在这基础上，还要特别是上机操作，要把理论应用于实际
完稿日期: 2010 年 12 月 31 日
参考文献
1. 程丽华，周玲丽. 数学软件教程. 广东：中山大学出版社，2008.7
2. 王能超，易大义，李庆扬. 数值分析. 北京：清华大学出版社，2008.12
附 录
源程序:
程序1
功能：求出矩阵p ，其中p 为用于：
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)); for k=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; for k=1:n-1 A(k,k+1)=r(k); A(k+1,k)=u(k); end m=A\d; ft=d(1);
1);
-n 0,1(k , 1)]x(k [x(k),t x(k)),-(t *p(k,4)t)-1)(x(k *p(k,3)x(k))^3-(t *p(k,2)t)^3-1)(x(k *p(k,1) =+∈+++++
syms t
for k=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);
for i=1:n-1 %出点xt对应的y值
for k=1:kmax
if x(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')
hold on
plot(x,y,'*')
hold off
程序2
功能：得出插值函数的kmax等分点xt和对应的插值函数值fx：
function [xt,fx]=scyt2(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));
for k=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;
for k=1:n-1
A(k,k+1)=r(k);
A(k+1,k)=u(k);
end
m=A\d;
ft=d(1);
syms t
for k=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=100;
xt=linspace(x(1),x(n),kmax);
for i=1:n-1 %出点xt对应的y值
for k=1:kmax
if x(i)<=xt(k)&xt(k)<=x(i+1)
fx(k)=subs(sx(i),xt(k));
end
end
end
程序3
功能：定义函数fx
function fx=myfx(x)
fx=1./(1+25*x.^2);
程序 4
功能：画出13等分、20等分和原图。x1=-1:0.1:1;
x2=-1:0.15:1;
x=-1:0.01:1;
syms t
y1=feval(@myfx,x1);
y2=feval(@myfx,x2);
y=feval(@myfx,x);
[xt1,fx1]=scyt2(x1,y1,50/26^2,-50/26^2); [xt2,fx2]=scyt2(x2,y2,50/26^2,-50/26^2); plot(xt1,fx1,'r');
xlabel('x');
ylabel('y');
title('y=1/(1+25*x^2)');
hold on
plot(x1,y1,'mo');
plot(x2,y2,'c*');
plot(xt2,fx2,'g--');
plot(x,y,'k');
legend('20等分图','20等分点','13等分点','13等分图','原图'); hold off
2 三次样条插值程序 三次样条插值利用方案二（求解固支样条或压紧样条） 按照要求要起点和终点的一阶导数值已知， 可得关于01,,.....,n M M M 的严格对角占优势的三对角方程组 然后利用三对角法（追赶法）解此线性方程组。 （1）编写M 文件，并保存文件名scfit.m % x,y 分别为n 个节点的横坐标和纵坐标值组成的向量 % dx0和dxn 分别为S 的导数在x0和xn 处的值，即m 0和m n n=length(x)-1; h=diff(x); d=diff(y)./h; a=h(2:n-1); b=2*(h(1:n-1)+h(2:n)); c=h(2:n); u=6*diff(d); b(1)=b(1)-h(1)/2; u(1)=u(1)-3*(d(1)-dx0); b(n-1)=b(n-1)-h(n)/2; u(n-1)=u(n-1)-3*(dxn-d(n)); %追赶法部分 for k=2:n-1 temp=a(k-1)/b(k-1); b(k)=b(k)-temp*c(k-1); u(k)=u(k)-temp*u(k-1); end m(n)=u(n-1)/b(n-1); for k=n-2:-1:1 m(k+1)=(u(k)-c(k)*m(k+2))/b(k); end %求S K1，S K2，S K3，S K4 m(1)=3*(d(1)-dx0)/h(1)-m(2)/2; m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2; for k=0:n-1 00 ()S x m '=()n n S x m '=0011111111212212n n n n n n M d M d M d M d μλμλ----??????????????????????=??????????????????????????
样条插值函数及应用
摘要 样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。 在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。 插值法是数值逼近的重要方法之一，它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中,有许多插值函数的曲线要求具有较高的光滑性,在整个曲线中,曲线不但不能有拐点,而且曲率也不能有突变。因此，对于插值函数必须二次连续可微且不变号 ,这就需要用到三次样条插值。 关键词三次样条函数；插值法
目录 引言 0 第一章三次样条插值 (1) 1.1 样条插值函数简介 (1) 1.2 三次样条函数应用 (2) 第二章AMCM91A 估计水塔水流量 (4) 2.1 理论分析及计算 (5) 2.2运用MATLAB软件计算 (8) 参考文献 (13)
引言 样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。上世纪四十年代，在研究数据处理的问题中引出了样条函数，例如，在1946年Schoenberg将样条引入数学,即所谓的样条函数，直到五十年代，还多应用于统计数据的处理方面，从六十年代起，在航空、造船、汽车等行业中，开始大量采用样条函数。 在我国，从六十年代末开始，从船体数学放样到飞机外形设计，逐渐出现了一个使用样，逐渐出现了一个使用样条函数的热潮，并推广到数据处理的许多问题中。 在实际生活中有许多计算问题对插值函数的光滑性有较高的要求，例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数，用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时其函数值整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数值，不会发生“龙格现象”。 现在国内外学者对这方面的研究也越来越重视，根据我们的需要来解决不同的问题，而且函数的形式也在不断地改进，长期以来很多学者致力于样条插值的研究，对三次样条的研究已相当成熟。
沈阳航空航天大学 数学软件课程设计 （设计程序） 题目三次样条插值函数 班级 / 学号 学生姓名 指导教师
沈阳航空航天大学 课程设计任务书 课程名称数学软件课程设计 院（系）理学院专业信息与计算科学 班级学号姓名 课程设计题目三次样条插值函数 课程设计时间: 2010 年12月20日至2010 年12月31日 课程设计的内容及要求： 1.三次样条插值函数 给出函数在互异点处的值分别为。 (1)掌握求三次样条插值函数的基本原理； (2)编写程序求在第一边界条件下函数的三次样条插值函数； (3)在区间上取n=10，20，分别用等距节点对函数 作三次样条插值函数，利用(1)的结果画出插值函数的图形，并在该图形界面中同时画出的图形。 [要求] 1.学习态度要认真，要积极参与课程设计，锻炼独立思考能力； 2.严格遵守上机时间安排； 3.按照MATLAB编程训练的任务要求来编写程序； 4.根据任务书来完成课程设计论文； 5.报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”； 6.报告上交时间：课程设计结束时上交报告；
7.严谨抄袭行为。 指导教师年月日负责教师年月日学生签字年月日
沈阳航空航天大学 课程设计成绩评定单 课程名称数学软件课程设计 院（系）理学院专业信息与计算科学课程设计题目三次样条插值函数 学号姓名 指导教师评语： 课程设计成绩 指导教师签字 年月日
目录 一正文 (1) 1问题分析 (1) 1.1 题目 (1) 1.2 分析 (1) 2 研究方法原理 (1) 2.1 求三次样条插值多项式，算法组织 (1) 3 算例结果 (3) 二总结 (7) 参考文献 (8) 附录 (9) 源程序: (9) 程序1 (9) 程序2 (10) 程序3 (12) 程序 4 (12)
对样条函数及其插值问题的一点认识 样条函数是计算数学以及计算机辅助设计几何设计的重要工具。1946年，I. J. Schoenberg 著名的关于一元样条函数的奠定性论文“Contribution to the problem of application of equidistant data by analytic functions ”发表，建立了一元样条函数的理论基础。自此以后，关于样条函数的研究工作逐渐深入。随着电子计算机技术的不断进步，样条函数的理论以及应用研究得到迅速的发展和广泛的应用。经过数学工作者的努力，已经形成了较为系统的理论体系。 所谓（多项式）样条函数，乃指具有一定光滑性的分段（分片）多项式。一元n 次且n －1阶连续可微的样条函数具有如下的表示式： 1()()()()N n n j j j s x p x c x x x +==+--∞<<+∞∑[] 011,00,01,,...,,(1),...,(),,...,,n n n n N n N N u un u u u u x x x x x S x x x x ++++ +≥??=??
告: Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值(表格查找。该命令对数据点之间计算内插值。它找出一元函数f(x在中间点的数值。其中函数f(x由所给数据决定。x:原始数据点 Y:原始数据点 xi:插值点 Yi:插值点 格式 (1yi = interp1(x,Y,xi 返回插值向量yi,每一元素对应于参量xi,同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵,则按Y 的每列计算。yi 是阶数为length(xi*size(Y,2的输出矩阵。 (2yi = interp1(Y,xi 假定x=1:N,其中N 为向量Y 的长度,或者为矩阵Y 的行数。 (3yi = interp1(x,Y,xi,method 用指定的算法计算插值: ’nearest’:最近邻点插值,直接完成计算; ’linear’:线性插值(缺省方式,直接完成计算;
’spline’:三次样条函数插值。对于该方法,命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函 数。命令spline 用它们执行三次样条函数插值; ’pchip’:分段三次Hermite 插值。对于该方法,命令interp1 调用函数p chip,用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形; ’cubic’:与’pchip’操作相同; ’v5cubic’:在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量,使用方法’nearest’、’linear’、’v5cubic’的插值算法,相应地将返回NaN。对其他的方法,interp1 将对超出的分量执行外插值算法。 (4yi = interp1(x,Y,xi,method,'extrap' 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5yi = interp1(x,Y,xi,method,extrapval 确定超出x 范围的xi 中的分量的外插值extrapval,其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x; 2.>>xx = 0:.25:10; yy = interp1(x,y,xx; 3.>>plot(x,y,'kd',xx,yy 复制代码 例2 1.>> year = 1900:10:2010;
例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表： 且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s 本算法求解出的三次样条插值函数将写成三弯矩方程的形式： ) ()6()() 6()(6)(6)(211123 13 1j j j j j j j j j j j j j j j j x x h h M y x x h h M y x x h M x x h M x s -- + -- + -+ -= +++++其中，方程中的系数 j j h M 6， j j h M 61+， j j j j h h M y )6(2- ， j j j j h h M y ) 6(211++- 将由Matlab 代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。 以下为Matlab 代码： %============================= % 本段代码解决作业题的例1 %============================= clear all clc % 自变量x 与因变量y ，两个边界条件的取值 IndVar = [0, 1, 2, 3]; DepVar = [0, 0.5, 2, 1.5];
LeftBoun = 0.2; RightBoun = -1; % 区间长度向量，其各元素为自变量各段的长度h = zeros(1, length(IndVar) - 1); for i = 1 : length(IndVar) - 1 h(i) = IndVar(i + 1) - IndVar(i); end % 为向量μ赋值 mu = zeros(1, length(h)); for i = 1 : length(mu) - 1 mu(i) = h(i) / (h(i) + h(i + 1)); end mu(i + 1) = 1; % 为向量λ赋值 lambda = zeros(1, length(h)); lambda(1) = 1; for i = 2 : length(lambda) lambda(i) = h(i) / (h(i - 1) + h(i));
数值计算方法作业 实验4.3 三次样条差值函数 实验目的： 掌握三次样条插值函数的三弯矩方法。 实验函数: dt e x f x t ? ∞ -- = 2 221)(π 实验内容： （1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值； （3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线 比较插值结果。 实验4.5 三次样条差值函数的收敛性 实验目的： 多项式插值不一定是收敛的，即插值的节点多，效果不一定好。对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。 实验内容： 按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。 实验要求： （1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情 况，分析所得结果并与拉格朗日插值多项式比较； （2） 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子，考
虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一 算法描述： 拉格朗日插值： 错误!未找到引用源。 其中错误!未找到引用源。是拉格朗日基函数，其表达式为：() ∏ ≠=--=n i j j j i j i x x x x x l 0) ()( 牛顿插值： ) )...()(](,...,,[.... ))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N 其中????? ?? ?? ?????? --=--= --= -)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i j i j i j i 三样条插值： 所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a CENTRAL SOUTH UNIVERSITY 数值分析实验报告 三次样条插值方法的应用 一、问题背景 分段低次插值函数往往具有很好的收敛性，计算过程简单，稳定性好，并且易于在在电子计算机上实现，但其光滑性较差，对于像高速飞机的机翼形线船体放样等型值线往往要求具有二阶光滑度，即有二阶连续导数，早期工程师制图时，把富有弹性的细长木条（即所谓的样条）用压铁固定在样点上，在其他地方让他自由弯曲，然后沿木条画下曲线，称为样条曲线。样条曲线实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续，从数学上加以概括就得到数学样条这一概念。下面我们讨论最常用的三次样条函数及其应用。 二、数学模型 样条函数可以给出光滑的插值曲线（面），因此在数值逼近、常微分方程和偏微分方程的数值解及科学和工程的计算中起着重要的作用。 设区间[]b ,a 上给定有关划分b x x n =<<<= 10x a ，S 为[]b ,a 上满足下面条件的函数。 ● )(b a C S ,2∈； ● S 在每个子区间[]1,+i i x x 上是三次多项式。 则称S 为关于划分的三次样条函数。常用的三次样条函数的边界条件有三种类型： ● Ⅰ型 ()()n n n f x S f x S ''0'',==。 ● Ⅱ型 ()()n n n f x S f x S ''''0'''',==，其特殊情况为()()0''''==n n x S x S 。 ● Ⅲ型 ()() 3,2,1,0,0==j x S x S n j j ，此条件称为周期样条函数。 鉴于Ⅱ型三次样条插值函数在实际应用中的重要地位，在此主要对它进行详细介绍。 三、算法及流程 按照传统的编程方法，可将公式直接转换为MATLAB 可是别的语言即可；另一种是运用矩阵运算，发挥MATLAB 在矩阵运算上的优势。两种方法都可以方便地得到结果。方法二更直观，但计算系数时要特别注意。这里计算的是方法一的程序，采用的是Ⅱ型边界条件，取名为spline2.m 。 Matlab 代码如下： function s=spline2(x0,y0,y21,y2n,x) %s=spline2(x0,y0,y21,y2n,x) %x0,y0 are existed points,x are insert points,y21,y2n are the second 学习报告—— 三次样条函数插值问题的讨论 班级：数学二班 学号：152111033 姓名：刘楠楠 样条函数： 由一些按照某种光滑条件分段拼接起来的多项式组成的函数；最常用的样条函数为三次样条函数，即由三次多项式组成，满足处处有二阶连续导数。 一、三次样条函数的定义： 对插值区间[,]a b 进行划分，设节点011n n a x x x x b -=<< <<=，若 函数2()[,]s x c a b ∈在每个小区间1[,]i i x x +上是三次多项式，则称其为三次样条函数。如果同时满足()()i i s x f x = (0,1,2)i n =，则称()s x 为()f x 在 [,]a b 上的三次样条函数。 二、三次样条函数的确定： 由定义可设：101212 1(),[,] (),[,]()(),[,] n n n s x x x x s x x x x s x s x x x x -∈??∈?=???∈?其中()k s x 为1[,]k k x x -上的三次 多项式，且满足11(),()k k k k k k s x y s x y --== (1,2,,k n = 由2()[,]s x C a b ∈可得：''''''()(),()(),k k k k s x s x s x s x -+-+== 有''1()(),k k k k s x s x -++= ''''1()(),(1 ,2,,1)k k k k s x s x k n -+ +==-， 已知每个()k s x 均为三次多项式，有四个待定系数，所以共有4n 个待定系数，需要4n 个方程才能求解。前面已经得到22(1)42n n n +-=-个方程，因此要唯一确定三次插值函数，还要附加2个条件，一般上，实际问题通常对样条函数在端点处的状态有要求，即所谓的边界条件。 1、第一类边界条件：给定函数在端点处的一阶导数，即 ''''00(),()n n s x f s x f == 2、第二类边界条件：给定函数在端点处的二阶导数，即 MATLAB 程序设计期中考查 在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题——三次样条插值的基本提法： 对插值区间[]b a ,进行划分：b x x x a n ≤ ＭATLＡＢ三次样条插值之三弯矩法 首先说这个程序并不完善，为了实现通用（１，２，…,ｎ）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。 首先，通过函数saｎwanｊ得到方程的系数矩阵，即追赶法方程的四个向量参数,接下来调 用追赶法(在inteｒｓaｎwj函数中），得到三次样条分段函数系数因子,然后进行多项式合并 得到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算 改点的值。附:追赶法程序cｈase ％%%％％%%%%%%%％% fｕncｔiｏｎ [nｅwv,w,newu，neｗd]=sanwj（ｘ,y,x０,y0,y1ａ,y1b）?％三弯矩样 条插值?％将插值点分两次输入，x0y０单独输入?% 边值条件a的二阶导数 y1a 和ｂ 的二阶导数ｙ1b，这里建议将y１ａ和y1b换成y2a和ｙ2b,以便于和三转角代码相区别 ?n=length(x);m=leｎgth(y); ｉf m~=n?errｏｒ(＇x or y 输入有误,再来'）; end?v=ｏneｓ(ｎ-１，1);u=ｏnes(n－１,1）;d=zeros（n－１,１）;?w=２*ｏ nes(n+１)；?h0=x（１）-x0;?h=zerｏs(n-１,1); for k＝1:n-１?ｈ(k)＝x（k+1)－x（k);?enｄ ｖ(1)＝h0/(ｈ0+h(1）); u(1)＝1－v（1）; ｄ（１)＝6*((y（2）-y(1）)/h（1)-(y(1)-y０)/ｈ0）/(h0+ｈ(１）);?% ｆoｒ k＝2:ｎ-1?ｖ(ｋ）＝h(ｋ-１)/(ｈ(k-１)+h(ｋ)）;?u(k)=1－v(k);?d（k)＝ 6*((y(ｋ+１)-ｙ（ｋ)）/h（k）-(ｙ(k）-y(ｋ-1）)/ｈ(k－1))／(h（k-1)+h(k)); ｅnd newv=[v;1］;?newu=［1;u]； d0=6＊（(y(1)-y0)/h0-y1a)/h0; d(n)=6＊（y1ｂ-(ｙ(n）－ｙ(n-1))／ｈ(n-1))／h(ｎ-1）; newd=［d0;d]; %%%%%%％％%％％％ ｆunctｉon iｎｔersａnｗj（x，y,ｘ0,y0,y1a，y1b) %三弯矩样条插值?%第一部分?n=leｎgtｈ(ｘ）;ｍ＝length(y); iｆ m~＝n?eｒror（'ｘｏｒy 输入有误，再来'); enｄ?%重新定义h?h=zｅｒｏs(n,1）; h(1)=x（1）-ｘ０; fｏr k=2:n h（k)=ｘ(ｋ)－x(k-1)；?end ％sptep1调用三弯矩函数?[a,b,c,d]=sanwj(x，ｙ,ｘ0,ｙ0,y１ａ,y1b); 例：已知一组数据点，编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据 0.25 0.3 0.39 0.45 0.53 0.5000 0.5477 0.6245 0.6708 0.7280 求解，其中边界条件为. 1）三次样条插值自然边界条件源程序： function s=spline3(x,y,dy1,dyn) %x为节点，y为节点函数值，dy1，dyn分别为x=0.25,0.53处的二阶导 m=length(x);n=length(y); if m~=n error('x or y输入有误') return end h=zeros(1,n-1); h(n-1)=x(n)-x(n-1); for k=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); end 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; for i=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)); end for i=2:n-1; A(i,i-1)=v(i-1); A(i,i+1)=u(i-1); end A(n,n-1)=1; A(1,2)=1; A=A+2*eye(n); M=zhuigf(A,g); %调用函数，追赶法求M fprintf('三次样条（三对角）插值的函数表达式\n'); syms X; for k=1:n-1 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); end s=s.'; s=vpa(s,4); %画三次样条插值函数图像 for i=1:n-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); hold on End plot(x,y); grid on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %调用的函数： %追赶法 function M=zhuigf(A,g) n=length(A); L=eye(n); U=zeros(n); for i=1:n-1 U(i,i+1)=A(i,i+1); end U(1,1)=A(1,1); for i=2: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); end Y(1)=g(1); for i=2:n Y(i)=g(i)-L(i,i-1)*Y(i-1); end M(n)=Y(n)/U(n,n); for i=n-1:-1:1 M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i); 告： Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。x：原始数据点 Y：原始数据点 xi：插值点 Yi：插值点 格式 (1)yi = interp1(x,Y,xi) 返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵，则按Y 的每列计算。yi 是阶数为length(xi)*size(Y,2)的输出矩阵。 (2)yi = interp1(Y,xi) 假定x=1:N，其中N 为向量Y 的长度，或者为矩阵Y 的行数。 (3)yi = interp1(x,Y,xi,method) 用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。对于该方法，命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函 数。命令spline 用它们执行三次样条函数插值； ’pchip’：分段三次Hermite 插值。对于该方法，命令interp1 调用函数p chip，用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形； ’cubic’：与’pchip’操作相同； ’v5cubic’：在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1 将对超出的分量执行外插值算法。 (4)yi = interp1(x,Y,xi,method,'extrap') 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5)yi = interp1(x,Y,xi,method,extrapval) 确定超出x 范围的xi 中的分量的外插值extrapval，其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x); 2.>>xx = 0:.25:10; yy = interp1(x,y,xx); 3.>>plot(x,y,'kd',xx,yy) 复制代码 例2 1.>> year = 1900:10:2010; 2.>> product = [75.995 91.972 105.711 12 3.203 131.669 150.697 179.323 203.212 226.505 x0=[0 0.9211 1.8431 2.9497 3.8714 4.9781 5.9 7.0064 7.9286 8.9678 10.9542 12.0328 12.9544 13.8758 14.9822 15.9039 16.8261 17.9317 19.0375 19.9594 20.8392 22.9581 23.88 24.9869 25.9083]; >> >> y0=[14405 11180 10063 11012 8797 9992 8124 10160 8488 11018 19469 20196 18941 15903 18055 15646 13741 14962 16653 14496 14648 15225 15264 13708 9633]; >> x=0:0.1:25.9; >> y1=interp1(x0,y0,x,'spline'); >> pp1=csape(x0,y0); %样条插值工具箱函数 y2=ppval(pp1,x); %计算x对应的y值 pp2=csape(x0,y0,'second'); y3=ppval(pp2,x); xydata=[x',y1',y2',y3'] subplot(1,2,1) plot(x0,y0,'+',x,y1) title('Spline1') subplot(1,2,2) plot(x0,y0,'+',x,y2) title('Spline2') dx=diff(x); dy=diff(y2); dy_dx=dy./dx; dy_dx0=dy_dx(1) ytemp=y2(13<=x&x<=15); ymin=min(ytemp); xmin=x(y2==ymin); xymin_1315=[xmin,ymin] 实验四三次样条插值的应用 一、问题描述 The upper portion of this noble beast is to be approximated using clamped cubic spline interpolants. The curve is drawn on a grid from which the table is constructed. Use Algorithm 3.5 to construct the three clamped cubic splines. 二、模型建立 三次样条插值 给定一个列表显示的函数 yi=y(xi),i=0,1,2,...,N-1。特别注意在xj和xj+1之间的一个特殊的区间。该区间的线性插值公式为： (3.3.1)式和(3.3.2)式是拉格朗日插值公式(3.1.1)的特殊情况。 因为它是(分段)线性的，(3.3.1)式在每一区间内的二阶导数为零，在横坐标为xj处的二阶导数不定义或无限。三次样条插值的目的就是要得到一个内插公式，不论在区间内亦或其边界上，其一阶导数平滑，二阶导数连续。 做一个与事实相反的个假设，除yi的列表值之外，我们还有函数二阶导数y"的列表值，即一系列的yi"值，则在每个区间内，可以在(3.3.1)式的右边加上一个三次多项式，其二阶导数从左边的yj"值线性变化到右边的yj+1"值，这么做便得到了所需的连续二阶导数。如果还将三次多项式构造在xj和xj+1处为零，则不会破坏在终点xj和xj+1处与列表函数值yj和yj+1的一致性。 进行一些辅助计算便可知，仅有一种办法才能进行这种构造，即用 注意,(3.3.3)式和(3.3.4)式对自变量x的依赖，是完全通过A和B对x的线性依赖，以及C和D(通过A和B)对x的三次依赖而实现。可以很容易地验证，y"事实上是该插值多项式的二阶导数。使用ABCD的定义对x求(3.3.3)式的导数，计算dA/dx dB/dx dC/dx dD/dx,结果为一阶导数 三次样条插值C/C++程序（自己整理的） 具体推导看书<<数值分析>> code： #include freopen("threeInsert.in", "r", stdin); scanf("%d", &n); for (i = 0; i <= n; i++) scanf("%lf%lf", &x[i], &y[i]); // scanf("%lf%lf", &m[0], &m[n]); for (i = 0; i <= n - 1; i++) h[i] = x[i + 1] - x[i]; //第一种边界条件 //alph[0] = 0; alph[n] = 1; beta[0] = 2 * m[0]; beta[n] = 2 * m[n]; //第二种边界条件 alph[0] = 1; alph[n] = 0; beta[0] = 3 * (y[1] - y[0]) / h[0]; beta[n] = 3 * (y[n] - y[n - 1] / h[n - 1]); for (i = 1; i <= n - 1; i++) { alph[i] = h[i - 1] / (h[i - 1] + h[i]); beta[i] = 3 * ((1 - alph[i]) * (y[i] - y[i - 1]) / h[i - 1] + alph[i] * (y[i + 1] - y[i]) / h[i]); } a[0] = - alph[0] / 2; b[0] = beta[0] / 2; for (i = 1; i <= n; i++) { a[i] = - alph[i] / (2 + (1 - alph[i]) * a[i - 1]); b[i] = (beta[i] - (1 - alph[i]) * b[i - 1]) / (2 + (1 - alph[i]) * a[i - 1]); } m[n + 1] = 0; for (i = n; i >= 0; i--) { m[i] = a[i] * m[i + 1] + b[i]; } scanf("%lf", &xx); for (i = 0; i < n; i++) { if (xx >= x[i] && xx <= x[i + 1]) break ; } printf("%lf\n", sunc(xx, i)); 一．介绍 二．程序框图 开始 输入未知数X及 (xi,yi),i=0,1,…,n 计算步长H[i] 计算λ、μ、d 根据边界条件，求 解相应的方程得到 M1,…, Mn 将M代入原方程， 得到分段函数 结束 三．源码 syms h n=9;%插入节点数，可以根据题目更改 h=2/(n+1); u=0.5; v=0.5; f=inline('1/(1+25*x.^2)');%输入函数，这个也可以根据题目更改g=inline('3/h*((c-b)/h-(b-a)/h)','a','b','c','h'); for i=1:n+2 x(1)=-1; x(i+1)=x(i)+2/(n+1); y(i)=f(x(i)); end for i=1:n d(i)=g(y(i),y(i+1),y(i+2),h); end A=zeros(n,n); for i=1:n A(i,i)=2; end for i=1:n-1 A(i,i+1)=u; A(i+1,i)=v; end B=zeros(n,1); for i=1:n B(i,1)=d(i) end C=inv(A)*B for i=1:n M(i)=C(i,1); end x=(-1:h/50:1); k=1./(1+25*x.^2); cs=spline(x,k); plot(x,k,'r.'); hold on; ezplot('1/(1+25*x^2)',[-1 1]); title('三次样条插值曲线和f(x)曲线') 四．结果 系数矩阵 弯矩M 分段函数不同次幂x对应的系数 三次样条插值函数与原函数 function [a b f]=spline3(x,y,flag,vl,vr) %三次样条插值函数 %(x,y)为插值节点，xx为插值点； %flag表端点边界条件类型： %flag=1:第一类边界条件(端点一阶导数给定)； %flag=2:第二类边界条件(端点二阶导数给定)； %flag=3:第三类边界条件； %vl,vr表左右端点处的在边界条件值。 %样条函数为：Si(x)=yi+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3 %b,c,d分别为各子区间上的系数值 %yy表插值点处的函数值. if length(x)~=length(y) error('输入数据应成对!'); end n=length(x); a=zeros(n-1,1); b=zeros(n-1,1); dx=a;dy=a; A=zeros(n);B=zeros(n,1); for i=1:n-1 dx(i)=x(i+1)-x(i); dy(i)=y(i+1)-y(i); end for i=2:n-1 A(i,i-1)=dx(i-1)/(dx(i-1)+dx(i)); A(i,i)=2; A(i,i+1)=dx(i)/(dx(i-1)+dx(i)); B(i,1)=6*(dy(i)/dx(i)-dy(i-1)/dx(i-1))/(dx(i)+dx(i-1)); end %---------------------------------% %端点一阶导数条件% if flag==1 A(1,1)=2; A(1,2)=1; A(n,n-1)=1; A(n,n)=2; B(1,1)=6*(dy(1)/dx(1)-vl)/dx(1); B(n,1)=6*(vr-dy(n-1)/dx(n-1))/dx(n-1); c=A\B; end %---------------% %端点二阶导数条件% if flag==2 A(1,1)=2; A(n,n)=2; B(1,1)=2*vl; B(n,1)=2*vr; c=A\B; end三次样条插值方法的应用
关于三次样条插值函数的学习报告(研究生)资料
(精选)三次样条插值的MATLAB实现
MATLAB三次样条插值之三弯矩法
三次样条插值自然边界条件
Matlab中插值函数汇总和使用说明
三次样条插值函数matlab程序绝不坑爹
实验四三次样条插值Word版
三次样条插值的C程序(很全啊)
三次样条插值函数
三次样条程序
相关主题
文本预览