《昆明理工大学数值分析上机报告5.doc》由会员分享,可在线阅读,更多相关《昆明理工大学数值分析上机报告5.doc(11页珍藏版)》请在三一办公上搜索。
1、数值分析实验报告姓 名: 学 号:2006202066专 业:材料学学 院:材料与冶金工程学院授课教师: 昆 明 理 工 大 学 06工科硕士 数值分析上机实验报告专业: 材料学 姓名:牛延龙 学号: 2006202066 任课教师: 作业完成实验室: 实验成绩:理论描述(10)数值公式(10)程序流程图和程序结构(20)数据和结果(20)讨论(20)源程序(20)总分(100)实验内容:1题目/要求:函数插值方法一、问题提出对于给定的一元函数 的n+1个节点值 。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下:(1)0.40.550.650.800.95
2、1.050.410750.578150.696750.901.001.25382求五次Lagrange多项式L,和分段三次插值多项式,计算 的值。(2) 12345670.3680.1350.0500.0180.0070.0020.001 试构造Lagrange多项式L,计算的值。 结果0.165299 0.00213348 二、要求1、 利用Lagrange插值公式 编写出插值多项式程序;2、 给出插值多项式或分段三次插值多项式的表达式;3、 根据节点选取原则,对问题(2)用三点插值或二点插值,其结果如何;4、 对此插值问题用Newton插值多项式其结果如何。2作业环境(包括选用的程序语言、
3、运行环境)本题中的插值多项式程序采用的编程语言为c+,因此运行环境可以在装有Microsoft VC+的windows XP 或2000的系统下运行程序。3数学(理论背景)描述在生产实践和科学研究所遇到的大量函数中,相当一部分是通过测量或实验得到的。虽然其函数关系y=f(x)在某个区间a,b上是客观存在的,但是却不知道具体的解析表达式,只能通过观察、测量或实验得到函数在区间a,b上一些离散点上的函数值、导数值等, 因此,希望对这样的函数用一个比较简单的函数表达式来近似地给出整体上的描述。还有些函数,虽然有明确的解析表达式,但却过于复杂而不便于进行理论分析和数值计算,同样希望构造一个既能反映函数
4、的特性又便于计算的简单函数,近似代替原来的函数。插值法就是寻求近似函数的方法之一。在用插值法寻求近似函数的过程中,根据所讨论问题的特点,对简单函数的类型可有不同的 选取,如多项式、有理式、三角函数等,其中多项式结构简单,并有良好的性质,便于数值计算和理论分析,因此被广泛采用。设函数y=f(x)在区间a,b上有定义且已知函数在区间a,b上n+1个互异点上的函数值,若存在一个简单函数y=p(x ),使其经过y=f(x)上的这n+1个已知点(),(),,(),即 p()= , i=0,1,n那么,函数p(x)称为插值函数,点称为插节点, 点(),(),,()称为插值点,包含插值节点的区间a,b称为插
5、值区间,求p (x)的方法称为插值法,f(x)称为被插函数。若p(x)是次数不超过n的多项式,用Pn(x)表示,即则称为n次插值多项式,相应的插值法称为多项式插值;若P(x)为分段多项式,称为分段插值,多项式插值和分段插值称为代数插值。4数值计算公式 1Lagrange插值公式: 2分段三次插值公式 其中;如下: 5算法程序流程与程序结构(程序中的函数调用关系)参数说明:X 双精度实型一维数组,长度为n。存放n个不等距结点的值(从小到大)。Y 双精度实型一维数组,长度为n。存放n个不等距结点上的函数值。n 整型变量。给定不等距结点的个数。t 双精度实型变量。指定插值点的值。k 整型变量。插值时
6、启始结点的位置。m 整型变量。插值时最后结点的位置。i,j 整型变量。数组下标。s 双精度实型变量,存放运算过程中间值。z 双精度实型变量,初始值为0.0。存放最终插值。h 双精度实型变量,给定n个结点的步长。算法描述:1) 全区不等值插值算法 输入数据,判断n值。如果=3,则进行全区间插值 For (i=0;in;i+) s=1.0 For (j=0;jn;j+) If (j!=i) s=s*(t-xj)/(xi-xj); z=z+s*yi; 最后返回z值。2) 全区间等值插值算法输入数据,判断n值。如果=3,则进行全区间插值 For (i=0; in;i+) s=1.0; xi=x0+i*
7、h; For (j=0;jn;j+) If (j!=i) xj=x0+j*h; s=s*(t-xj)/(xi-xj); z=z+s*yi; 最后返回z值。3) 三点等值插值算法输入数据,判断n值。如果=3,则判断t值 如果t=x0+(n-3)*h, 则取最后三个结点k=n-3;m=n-1;进行三点二次抛物插值;否则 取离t最近的中间三个结点进行插值 i=(int)(t-x0/h)+1; If (fabs(t-x0-i*h)=fabs(t-x0-(i-1)*h) k-i-2;m=I; Else k=i-1;m=m+1;三点二次插值程序 For (i=k; i=m;i+) s=1.0; xi=x0
8、+i*h; For (j=k;j=m;j+) If (j!=i) xj=x0+j*h; s=s*(t-xj)/(xi-xj); z=z+s*yi; 最后返回z值。4) 函数调用关系如下:主函数main() 输入数据; 全区不等插值函数nlgr()或全区等值插值函数 elgr()或三点等值插值函数elg3(); 输出;6实验数据和实验结果(打印或用屏幕图形拷屏表示,可加为附页)对于(1)中:其中y0=0.41075, y1=0.57815, y2=0.69675, y3=0.90, y4=1.00, y5=1.25382; 对于(2)中:其中y0=0.368, y1=0.135, y2=0.05
9、0, y3=0.018, y4=0.007, y5=0.002, y0=0.001;7讨论(包括题目要求的讨论和方法的适用性讨论)1)三点插值与二点插值比较根据结点选取原则,对问题(2)用三点插值的结果如上6中所示为1.6976e-001;二用两点插值的结果为0.1816;与题中所给的准确值相比,可以看出明显是三点插值更接近准确值。但两点插值比三点插值算法更简单些。2)拉格朗日插值与牛顿插值的比较 两者都是通过给定n+1个互异的插值节点,让你求一条n次代数曲线近似地表示待插值的函数曲线Lagrange插值代数和Newton法插值都属于代数插值的范畴区别:Lagrange插值法是通过构造n+1个
10、n次基本多项式,然后线性组合而得到的,而Newton法插值是通过求各阶差商,递推得到的一个 f(x)=f(x0)+(x-x0)fx0,x1+(x-x0)(x-x1)fx0,x1,x2.+(x-x0).(x-x(n-1)fx0,x1.xn这样的公式,代进去就可以得到。Lagrange插值法在求每个基本多项式的时候要用到所有那些结点,因此如果需要再多加进去一个结点的话,需要重新求出基本多项式才可,而这需要大量的工程。而对Newton插值,如果再加进去一个结点是不是只要在它后面再加上一个(x-x0)(x-x1).(x-x(n-1)(x-xn)fx0,x1.xn,x(n+1)就行了。并且,它每一项的系
11、数就是各阶差商,比拉格朗日插值计算量小,且便于编程。 所以,拉格朗日插值公式和牛顿插值公式本质上是同一个多项式,只是形式不同而已,对此问题如用newton插值多项式计算其结果会相同。8附源程序源程序:1) 全区间不等值插值程序:#includemath.h#includestdio.hdouble nlg3(double x,double y,int n,double t)/int n;/double t,x,y;int i,j;double z, s; z=0.0; if(n1) return(z);if(n=1) z=y0;return(z);if(n=2)z=double(y0*(t-x
12、1)-y1*(t-x0)/(x0-x1);return(z);z=0.0;for(i=0;in;i+)s=1.0;for(j=0;jn;j+)if(j!=i) s=s*(t-xj)/(xi-xj);z=z+s*yi;return(z);void main()double t,z;static double x7=0.4,0.55,0.65,0.80,0.95,1.05;static double y7=0.41075,0.57815,0.69675,0.90,1.00,1.25382;printf(n);t=0.596;z=nlg3(x,y,6,t);printf(x=%6.3f, f(x)=
13、%en,t,z);t=0.99;z=nlg3(x,y,6,t);printf(x=%6.3f, f(x)=%en,t,z);printf(n);2) 全区间等值插值程序:#include math.h#include stdio.hdouble elg3(double x0,double h,int n,double y,double t)int i,j;double z,s,xi,xj;z=0.0;if(n1) return(z);if(n=1) z=y0;return(z);if(n=2)z=(y1*(t-x0)-y0*(t-x0-h)/h;return(z);for(i=0;in;i+)
14、s=1.0;xi=x0+i*h;for(j=0;jn;j+)if(j!=i)xj=x0+j*h;s=s*(t-xj)/(xi-xj);z=z+s*yi;return(z);void main()double t,z,x0,h;static double y7=0.368,0.135,0.050,0.018,0.007,0.002,0.001;x0=1.0;h=1.0;printf(n);t=1.8; z=elg3(x0,h,7,y,t);printf(x=%f, f(x)=%en,t,z);t=6.15; z=elg3(x0,h,7,y,t);printf(x=%f, f(x)=%en,t,z
15、);printf(n);3) 三点等值二次插值程序:#include math.h#include stdio.hdouble elg3(double x0,double h,int n,double y,double t)int i,j,k,m;double z,s,xi,xj;z=0.0;if(n1) return(z);if(n=1) z=y0;return(z);if(n=2)z=(y1*(t-x0)-y0*(t-x0-h)/h;return(z);if(t=x0+(n-3)*h) k=n-3;m=n-1;else i=(int)(t-x0/h)+1;if(fabs(t-x0-i*h)
16、=fabs(t-x0-(i-1)*h) k=i-2;m=i;else k=i-1;m=i+1;z=0.0;for(i=k;i=m;i+)s=1.0;xi=x0+i*h;for(j=k;j=m;j+)if(j!=i)xj=x0+j*h;s=s*(t-xj)/(xi-xj);z=z+s*yi;return(z);void main()double t,z,x0,h;static double y7=0.368,0.135,0.050,0.018,0.007,0.002,0.001;x0=1.0;h=1.0;printf(n);t=1.8; z=elg3(x0,h,7,y,t);printf(x=%6.3f, f(x)=%en,t,z);t=6.15; z=elg3(x0,h,7,y,t);printf(x=%6.3f, f(x)=%en,t,z);printf(n);第11页