沈阳航空航天大学数学软件课程设计（设计程序）题目三次样条插值函数班级 / 学号学生姓名指导教师沈阳航空航天大学课程设计任务书课程名称数学软件课程设计院（系）理学院专业信息与计算科学班级学号姓名课程设计题目三次样条插值函数课程设计时间: 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 -=<<⋅⋅⋅<<=分成的每个小区间1[,]i i x x -上是3次多项式，其在此区间上的表达式如下：22331111111()[()()]()()666[,]1,2,,.i i i i i i i i i i i i i i ii 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 的计算方法如下： 令：11111111116()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 μλ-+++==⋅⋅⋅-，对于第一种边界条件下有⎪⎪⎩⎪⎪⎨⎧-=+-=+---000110111)'],([62]),['(62h f x x M M h x x f f M M n n n n n n如果令,]),['(6,1,)'],[(6,111000100---==-==n n n n n n h x x f f d h f x x f d μλ那么解就可以为⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n nn n d d d d M M M M 110110111102222μλμλμλ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)(xs当把区间20等分时，输入如下：图7得到的矩阵p如下：图8得出的图像如下：图9 20等分图像运行gtu.m，输入如下：图10运行结果为图像，图形如下：图11 13、20等分和原图的图像二总结拿到题目时，我首先先弄清三次样条插值函数的基本原理，因为只有这样，在以后的编程中才会更懂的如何编写程序，才会更不会混淆题目的目的。

而且，不能为了做题而做题，在做题时还要应用到其它知识。

三次样条插值在不懂的边界条件出来的结果也不一样，表达形式也不一样，特别是计算过程，条件不一样，编写的程序也自然不同，题目给出的是第一边界条件，在第一边界条件中，首先要给出始末两点的导数值，而且要给出n个点的x值和y的对应的值，题目要求，所以有在程序输入x,y和始末两点的导数值。

这课设过程中，我学到了很多，知道了自己的很多不足，如对于某些函数不能巧妙的应用，编写的程序很粗糙，这次课设，我们通过自己的思考与实践，终于完成了。

刚开始我对于三次样条插值很不了解，现在，我已经对于三次样条插值有了一定的了解，特别是其中的原理，但公式还是没背下来，不过，原理最重要，这次课设，更好的让我们掌握了Matlab和数值分析，我了解到理论与实践是分不开的，为了更好的掌握一个知识，我们必须通过不断的实践。

要学好一门计算机语言，既要掌握其最基本的语言结构，而且要特别的熟悉，在这基础上，还要特别是上机操作，要把理论应用于实际完稿日期: 2010 年 12 月 31 日参考文献1. 程丽华，周玲丽. 数学软件教程. 广东：中山大学出版社，2008.72. 王能超，易大义，李庆扬. 数值分析. 北京：清华大学出版社，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-1u(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)); endA=eye(n,n)*2;for k=1:n-1A(k,k+1)=r(k);A(k+1,k)=u(k);endm=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 tfor 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)); endkmax=1000;xt=linspace(x(1),x(n),kmax);for i=1:n-1 %出点xt对应的y值for k=1:kmaxif x(i)<=xt(k)&xt(k)<=x(i+1)fx(k)=subs(sx(i),xt(k));endendendplot(xt,fx,'r'); xlabel('x');ylabel('y'); title('f');text(x(fix(n/2)),y(fix(n/2)),'f')hold onplot(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-1u(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)); endA=eye(n,n)*2;for k=1:n-1A(k,k+1)=r(k);A(k+1,k)=u(k);endm=A\d;ft=d(1);syms tfor 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)); endkmax=100;xt=linspace(x(1),x(n),kmax);for i=1:n-1 %出点xt对应的y值for k=1:kmaxif x(i)<=xt(k)&xt(k)<=x(i+1)fx(k)=subs(sx(i),xt(k));endendend程序3功能：定义函数fxfunction fx=myfx(x)fx=1./(1+25*x.^2);程序 4功能：画出13等分、20等分和原图。