三次样条插值函数.doc

沈阳航空航天大学数学软件课程设计（设计程序）题目 三次样条插值函数 班级 / 学号 学 生 姓 名 指 导 教 师 沈阳航空航天大学课 程 设 计 任 务 书课 程 名 称 数学软件课程设计 院（系） 理学院 专业 信息与计算科学 班级 学号 姓名 课程设计题目 三次样条插值函数 课程设计时间: 2010 年 12月 20日至 2010 年 12月 31日课程设计的内容及要求：1. 三次样条插值函数给出函数在互异点处的值分别为。(1) 掌握求三次样条插值函数的基本原理；(2) 编写程序求在第一边界条件下函数的三次样条插值函数；(3) 在区间上取n=10，20，分别用等距节点对函数作三次样条插值函数，利用(1)的结果画出插值函数的图形，并在该图形界面中同时画出的图形。要求1. 学习态度要认真，要积极参与课程设计，锻炼独立思考能力；2. 严格遵守上机时间安排；3. 按照MATLAB编程训练的任务要求来编写程序；4. 根据任务书来完成课程设计论文；5. 报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”；6. 报告上交时间：课程设计结束时上交报告；7. 严谨抄袭行为。指导教师 年 月 日负责教师 年 月 日学生签字 年 月 日沈阳航空航天大学课 程 设 计 成 绩 评 定 单课 程 名 称 数学软件课程设计 院（系） 理学院 专业 信息与计算科学 课程设计题目 三次样条插值函数 学号 姓名 指导教师评语：课程设计成绩 指导教师签字 年 月 日数学软件课程设计 沈阳航空航天大学 课程设计用纸 目 录目 录一 正文11问题分析11.1 题目11.2 分析12 研究方法原理12.1 求三次样条插值多项式，算法组织13 算例结果2二 总结7参考文献8附 录9源程序:9程序19程序210程序312程序 412第 7 页数学软件课程设计 沈阳航空航天大学 课程设计用纸 正 文一 正文1问题分析1.1 题目三次样条插值函数给出函数在互异点处的值分别为。（1）掌握求三次样条插值函数的基本原理；（2）编写程序求在第一边界条件下函数的三次样条插值函数；（3）在区间上取n=10，20，分别用等距节点对函数 作三次样条插值函数，利用(1)的结果画出插值函数图形，并在该图形界面中同时画出的图形。1.2 分析一般认为插值次数n越高，的精度就越高，但实际并非如此，20世纪初龙格(Runge)就发现了这一现象，因此就提出了分段低次插值分段线性插值有一致收敛性，但光滑性差，而三次样条插值具有二介光滑度，三次样条插值首先要给定n个点和对应的函数值，还要给出边界条件如第一边界条件，第二边界条件，而题目要求是在给定第一边界条件下的三次样条插值。2 研究方法原理2.1 求三次样条插值多项式，算法组织所谓三次样条插值多项式是一种分段函数，它在节点分成的每个小区间上是3次多项式，其在此区间上的表达式如下：因此，只要确定了的值，就确定了整个表达式，的计算方法如下：令：则满足如下n-1个方程：对于第一种边界条件下有如果令那么解就可以为 3 算例结果s(x)可以表示为：其中p为的矩阵。当把区间5等分时，输入如下：图 1矩阵p输出如下图 2 图形如下：图 3 5等分图像其在不同的区间的函数可以表示如下： 当n=10即区间10等分时，输入如下： 图 4得到的矩阵p如下：图 5图形如下：图 6 10等分图像当把区间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第 13 页数学软件课程设计 沈阳航空航天大学 课程设计用纸 附 录附 录源程序:程序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);endA=eye(n,n)*2;for k=1:n-1 A(k,k+1)=r(k); A(k+1,k)=u(k);endm=Ad;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=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); 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-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);endA=eye(n,n)*2;for k=1:n-1 A(k,k+1)=r(k); A(k+1,k)=u(k);endm=Ad;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:kmax if 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等分和原图。x1=-1:0.1:1; x2=-1:0.15:1; x=-1:0.01:1;syms ty1=feval(myfx,x1);y2=feval(myfx,x2);y=feval(myfx,x);xt1,fx1=scyt2(x1,y1