昆明理工大学数值分析上机报告2

1、函数插值方法(课题五)、问题提出对于给定的一元函数 y = f(x)的n+1个节点值 Yj = f(Xj )j=(0,1，n)。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。(2)数据如下：(Dxj0.40.550.650.800.951.05Yj0.4107 50.5781 50.6967 50.901.001.25382求五次Lagrange多项式L5(x),和分段三次插值多项式，计算xj1234567Yj0.3680.1350.0500.0180.0070.0020.001f (0.596 )f(0.99 )的值。试构造Lagrange多项式L6(x),计算

2、f(1.8 )的值结果 f 1.80.165299 f 6.150.00213348、wl1、利用Lagrange插值公式nLn X = 'k =0nni =0x - xiYk xk -xi编写出插值多项式程序;2、给出插值多项式或分段三次插值多项式的表达式;3、根据节点选取原则，对问题(2)用三点插值或二点插值，其结果如何;4、对此插值问题用Newton插值多项式其结果如何。2 .作业环境(包括选用的程序语言、运行环境)本题中的插值多项式程序采用的编程语言为C+,因此运行环境可以在装有Microsoft VC+的windows XP或2000的系统下运行程序。3 .数学(理论背景)描

3、述在生产实践和科学研究所遇到的大量函数中， 相当一部分是通过测量或实验得到的。虽然其函数关系y=f(x)在某个区间a, b上是客观存在的，但是却 不知道具体的解析表达式，只能通过观察、测量或实验得到函数在区间a, b 上一些离散点上的函数值、导数值等，因此，希望对这样的函数用一个比较简单的函数表达式来近似地给出整体上的描述。还有些函数，虽然有明确的解析表达式，但却过于复杂而不便于进行理论分析和数值计算，同样希望构造一个既能反映函数的特性又便于计算的简单函数，近似代替原来的函数。插值法就是寻求近似函数的方法之一。在用插值法寻求近似函数的过程中，根据所讨论问题的特点，对简单函数的类型可有不同的 选

4、取，如多项式、有理式、三角函数等，其中多项式结构 简单，并有良好的性质，便于数值计算和理论分析，因此被广泛采用。设函数y=f(x)在区间a,b上有定义，丫。，丫1,yn且已知函数在区间a,b上n+1个互异点*。，卬.方上的函数值，若存在一个简单函数尸p(x),使其经过 y=f(x)上的这 n+1 个已知点(Xo,y0),(Xl，yi),；(xn,yn),即p(Xi)= yi,i=0,1,n那么，函数p(x)称为插值函数，点Xo，X1，.Xn称为插节点，点(Xo，yo),(X1,y1), ； (Xn，yn)称为插值点，包含插值节点的区间a，b称为 插值区间，求p (X)的方法称为插值法，f(x)

5、称为被插函数。若p(x)是次数不超过 n的多项式，用Pn(x)表示，即pn(x)= a。a1x a2x2 anxn则称pn(x)为n次插值多项式，相应的插值法称为多项式插值；若 P(x)为分段 多项式，称为分段插值，多项式插值和分段插值称为代数插值。4 .数值计算公式1. Lagrange插值公式：Ln(X )= g 口x -xii =0 xk Xi京 Jyk2.分段三次插值公式nIh(x)f- i (x ) 1fi i =0i (x),其中 Ih(xi) =fi,I h(x) = f i(i =01An) ; %(x), Pi(x)如下:X _x-xi -Xi 1 x -x：1 2-xi为x

6、 j - x _ xi (i = 0略去)二 i(x)=（x）=x -xi 1x-x1Vx - Xx -x 1x -x1 +2-(xx(xx、x -X V0, x在其余地方x _ x _ x 1(i = 0略去),x.-x - xi (i = 0略去)，为-x - x, 1(i = n略去)5 .算法程序流程与程序结构（程序中的函数调用关系）参数说明:双精度实型一维数组，长度为no存放n个不等距结点的值（从小到大）双精度实型一维数组，长度为n0存放n个不等距结点上的函数值。整型变量。给定不等距结点的个数。双精度实型变量。指定插值点的值。整型变量。插值时启始结点的位置。整型变量。插值时最后结点的

7、位置。i,j整型变量。数组下标双精度实型变量，存放运算过程中间值z 双精度实型变量，初始值为0.0。存放最终插值h 双精度实型变量，给定n个结点的步长。算法描述：1)全区不等值插值算法输入数据，判断n值。如果<1返回z；如果=1则返回y0;如果=2则进行两点线性插值z=(y1*(t-x0)-y0*(t-x0-h)/hReturn(z)；若n>=3,则进行全区间插值For (i=0;i<n;i+) s=1.0For (j=0;j<n;j+)If (j!=i) s=s*(t-xj)/(xi-xj);z=z+s*yi;最后返回z值。2)全区间等值插值算法输入数据，判断n值。如

8、果<1返回z；如果=1则返回y0;如果=2则进行两点线性插值z=(y1*(t-x0)-y0*(t-x0-h)/hReturn(z)；若n>=3,则进行全区间插值For (i=0; i<n;i+) s=1.0; xi=x0+i*h;For (j=0;j<n;j+)If (j!=i) xj=x0+j*h; s=s*(t-xj)/(xi-xj); z=z+s*yi;最后返回z值。3)三点等值插值算法输入数据，判断n值。如果<1返回z；如果=1则返回y0;如果=2则进行两点线性插值z=(y1*(t-x0)-y0*(t-x0-h)/hReturn(z)；如果n>=3,

9、则判断t值如果t<=x0+h ,则取前面三个结点k=0;m=2;进行三点二次抛物 插值；如果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+i*h;For (j=k;j<=m;j+)If (j!=i) xj=x0+j*h; s=s*(t-

10、xj)/(xi-xj);z=z+s*yi;最后返回z值。4)函数调用关系如下:主函数 main ()输入数据；全区不等插值函数信函数elg3 ();输出；)6.实验数据和实验结果对于(1)中：nlgr ()或全区等值插值函数elgr ()或三点等值插(打印或用屏幕图形拷屏表示，可加为附页)5L5(x)，5 x- xi L)yk- y0l0(x)y1l1(x)y2l2(x)y3l3(x)y4l4(x)y5l5(x)其中 y0=0.41075, yi=0.57815, y2=0.69675, y3=0.90, y4=1.00, y5=1.25382;一、 (x - 0.5 5()-(0x6-5 )

11、 ( x 0.8 0 )x： - 0 .9 5 ) (1.0 5 )l0(x)=(0.4 -0.55)(0.4 -0.65)(0.4 -0.80)(0.4 -0.95)(0.4 -1.05)li(x)="x)l3 (x)l4(x)(x 一 0.4)(x -0.65)( x 一 0.80)( x 一0.95)( x 一 1.05)(0.55 -0.4)(0.55 -0.65)(0.55 -0.80)(0.55 -0.95)(0.55 -1.05)(x -0.4)(x -0.55)(x -0.80)(x-0.95)(x -1.05)(0.65 -0.4)(0.65 -0.55)(0.65

12、 -0.80)(0.65 -0.95)(0.65 -1.05)(x 0.4)(x0.55)(x0.65)(x0.95)(x1.05)(0.8 -0.4)(0.8 -0.55)(0.8 -0.65)(0.8 -0.95)(0.8 -1.05)(x0.4)(x 0.55)(x0.65)(x0.80)(x1.05)(0.95 -0.4)(0.95 -0.55)(0.95 -0.65)(0.95 -0.8)(0.95 -1.05)l5(x)(x -0.4)(x -0.55)(x-0.65)(x-0.80)(x -0.95)(1.05 -0.4)(1.05-0.55)(1.05 -0.65)(1.05

13、-0.8)(1.05 -0.95)对于(2)中:L6(x)八6 x - xi (:k =0i =0 xk 一 xii ；k)Yk=y°l0(x) y/Kx) y2"x) y，") y4x) ys/x) yeU)其中 y0=0.368, y1=0.135, y2=0.050, y3=0.018, y4=0.007, y5=0.002, y0=0.001;l0(x)=(x -2)(x -3)(x -4)(x -5)(x-6)(x-7)(1 -2)(1 -3)(1 -4)(1 -5)(1 -6)(1 -7)l1(x)-(x -1)(x -3)(x -4)(x -5)(x

14、 -6)(x -7)(2 -1)(2 -3)(2 -4)(2 -5)(2 -6)(2 -7)"x)(x -1)(x -2)(x -4)(x -5)(x-6)(x-7)(3 -1)(3 -2)(3 -4)(3 -5)(3-6)(3 -7)l3 (x)(x -1)(x-2)(x -3)(x -5)(x-6)(x -7)(4 -1)(4 -2)(4 -3)(4 -5)(4 -6)(4 -7)L(x)(x-1)(x-2)(x-3)(x-4)(x-6)(x-7)(5-1)(5-2)(5 -3)(5-4)(5-6)(5 -7)l5(x)(x -1)(x-2)(x -3)(x -4)(x-5)(x

15、-7)(6 -1)(6 -2)(6 -3)(6 -4)(6 -5)(6 -7)l6(x)(x -1)(x -2)(x -3)(x -4)(x-5)(x -6)(7 -1)(7 -2)(7 -3)(7 -4)(7 -5)(7 -6)上机实验截图：1)全区间不等值插值结果:2）全区间等值插值结果:3）三点等值二次插值结果:7 .讨论(包括题目要求的讨论和方法的适用性讨论)1)三点插值与二点插值比较根据结点选取原则，对问题( 2)用三点插值的结果如上6中所示为1.6976e-001;二用两点插值的结果为0.1816;与题中所给的准确值相比，可以 看出明显是三点插值更接近准确值。但两点插值比三点插值算

16、法更简单些。2)拉格朗日插值与牛顿插值的比较两者都是通过给定n+1个互异的插值节点，让你求一条 n次代数曲线近似 地表示待插值的函数曲线Lagrange插值代数和Newton法插值都属于代数插值 的范畴.区别：Lagrange插值法是通过构造n+1个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.x n这样的公式，代进去就可以得到。Lagrange插值法在求每个基本多项式的时 候要用到所有那些结点，因此如果需要

17、再多加进去一个结点的话， 需要重新求 出基本多项式才可，而这需要大量的工程。而对 Newton插值，如果再加进去 一个结点是不是只要在它后面再加上一个(x-x0)(x-x1).(x-x(n-1)(x-xn)fx0,x1.xn,x(n+1)就行了。并且，它每一项的系数就是各阶差商，比拉格朗日插值计算量小，且便于编程。所以，拉格朗日插值公式和牛顿插值公式本质上是同一个多项式，只是形式不同而已，对此问题如用newton插值多项式计算其结果会相同。8 .附一源程序(打包邮件)源程序：1)全区间不等值插值程序：#include"math.h"#include"stdio.h

18、"double 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(n<1) return(z);if(n=1) z=y0;return(z);if(n=2) z=double(y0*(t-x1)-y1*(t-x0)/(x0-x1);return(z);第10页z=0.0;for(i=0;i<n;i+)(s=1.0;for(j=0;j<n;j+)if(j!=i) s=s*(t-xj)/(xi-xj);z=z+s*yi;return(z);voi

19、d 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)=%en",t,z);t=0.99;z=nlg3(x,y,6,t);printf("x=%6.3f,f(x)=%en",t,z);printf("n");2)

20、全区间等值插值程序：#include "math.h"#include "stdio.h"double elg3(double x0,double h,int n,double y,double t)第11页int i,j;double z,s,xi,xj;z=0.0;if(n<1) 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;i<n;i+)s=1.0;xi=x0+i*h;for(j=0;j<n;j+)if(j!

21、=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;第12页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);printf("n");3)三点等值二次插值程序：#include "math.h"#include "stdio.h"double elg3(double x0,doub