数值分析报告.doc

实验报告一、对列表函数的计算【开发语言及实现平台或实验环境】C+/C#1.实验目的（1）掌握拉格郎日插值多项式的用法，适用范围及精确度。 （2）掌握牛顿插值多项式的用法，适用范围及精确度。（3）掌握埃特金插值多项式的用法，适用范围及精确度。2.算法思想(1)拉格朗日插值算法已知函数y=f(x)在n+1个不同的点x0 ,x1 ,x2上的函数值分别为y0 ,y1 ,yn ,求一个次数不超过n的多项式Pn (x),使其满足:Pn (xi )=yi , (i=0,1,n),即n+1个不同的点可以唯一决定一个n次多项式。插值基函数过n+1个不同的点分别决定n+1个n次插值基函数l0 (x),l1 (x),ln (X)每个插值基本多项式li (x)满足：a. li (x)是n次多项式;b. li (xi )=1,而在其它n个li (xk )=0 ,(ki)。由于li (xk )=0 ,(ki), 故有因子: (x-x0 )(x-xi-1 )(x-xi+1 )(x-xn )因其已经是n次多项式，故而仅相差一个常数因子。令:li (x)=a(x-x0 )(x-xi-1 )(x-xi+1 )(x-xn )由li (xi )=1,可以定出a, 进而得到： n次拉格朗日型插值多项式Pn (x)Pn (x)是n+1个n次插值基本多项式l0 (x),l1 (x),ln (X)的线性组合,相应的组合系数是y0 ,y1 ,yn。即:Pn (x)=y0 l0 (x)+y1 l1 (x)+yn ln (x),从而Pn (x)是一个次数不超过n的多项式,且满足Pn (xi )=yi , (i=0,1,2,n).(2)牛顿插值算法牛顿插值是基于下面这些的公式：fx0,x1,.xk=(fx1,.xk-fx0,.xk-1)/(xk-x0)fx=f(x)f(x)=fx0+fx0,x1(x-x0)+fx0,x1,x2(x-x0)(x-x1)+.fx0,.xn(x-x0).(x-xn-1)+Rn(x)前两个是均差的递推关系式，而后一个就是牛顿插值公式，其中N(x)=f(x)-Rn(x)，即目标多项式，Rn(x)是n阶插值余项，我们就是用N(x)去近似f(x)。可以构造这样一个均方差表：xk f(xk) 一阶均差 二阶均差 .x0 f(x0)x1 f(x1) fx0,x1x2 f(x2) fx1,x2 fx0,x1,x2.如果有n个点插值，表会有(n*n)/2+n个表项，如果直接编程会有O(n*n)的空间复杂度，编程时做个简单的改进，不难发现在这个表中只有部分数据有用，对角线（斜行）它们是目标值，用来表示多项式的，左边的两纵行（实际上只需要x一行）以及最底下的一行，表示当前插值的状态。经过改进后只需要O(n)的空间复杂度。两个过程：新增加一个点时的更新。只须更新最底下一行数据，其递推关系由均差公式给出，最后算出高一队的均差值，需时O(n)插入点完成后如何计算多项式在另外给定点的值N(x)。由牛顿插值公式，最终的表达式为：N(x)=fx0+fx0,x1(x-x0)+fx0,x1,x2(x-x0)(x-x1)+.fx0,.xn(x-x0).(x-xn-1)如果直接将它展开，再算实在麻烦，实际上大可不必这样做，还记得多项式求值的秦九韶算法吗？将多项式叠起来，从内层括号往外一层层拨开，n次多项多的计算，只需要做n次乘法，同样的思想，将上式改写成：N(x)=fx0+(x-x0)fx0,x1+(x-x1)fx0,x1,x2+(x-x2).fx0,.xn-1+(x-xn-1)fx0,.xn.就可以同样简单的计算了，时间复杂度O(n)综合起来的性能：对于n个点的插值，产生多项式的时间复杂度是O(n*n)，最终进行一个点的计算的时间复杂度是O(n)。(3)埃特金插值算法埃特金插值又称分段插值。分段插值是将被插值函数逐步多项式化。分段插值的处理过程分两步，将区间分成几个子段，并在每个子段上构造插值多项式装配在一起，作为整个区间的插值函数。在分化的每个节点给出数据，连接相邻节点得一折线，该折线函数可以视作插值问题的解。 举例说明：例 已知列表函数 1 2 3 4 0 5 6 3用Aitken算法求的近似值。解 ，用列表方法求解。故3.插值算法流程图（1）拉格朗日插值法（2）牛顿插值法（3）埃特金插值法4.实验内容X01234F（x）121782257求任意结点x对应的函数值。5.插值算法代码#include#include#include#includetypedef struct data float x; float y;Data;/变量x和函数值y的结构Data d20;/最多二十组数据float f(int s,int t)/牛顿插值法，用以返回插商 if(t=s+1) return (dt.y-ds.y)/(dt.x-ds.x); else return (f(s+1,t)-f(s,t-1)/(dt.x-ds.x); float Newton(float x,int count) int n; while(1) coutn; if(n=count-1)/插值次数不得大于count1次 break; /初始化t，y，yt。 float t=1.0; double y=d0.y; float yt=0.0;/计算y值 for(int j=1;j=n;j+) t=(x-dj-1.x)*t; yt=f(0,j)*t; /coutf(0,j)endl; y=y+yt; return y;float lagrange(float x,int count) double y=0.0; for(int k=0;kcount;k+)/这儿默认为count1次插值 float p=1.0;/初始化p for(int j=0;jcount;j+) /计算p的值 if(k=j)continue;/判断是否为同一个数 p=p*(x-dj.x)/(dk.x-dj.x); y=y+p*dk.y;/求和 return y;/返回y的值float aitken(float x0,int count)double y=0.0;int k=1;while(k!=count)for(int i=k;icount;i+)di.y=(x0-dk-1.x)/(di.x-dk-1.x)*di.y+(x0-di.x)/(dk-1.x-di.x)*dk-1.y;k+;y=dcount-1.y;return y;void main() float x,y; int count; while(1) cout请输入xi,yi的组数，不得超过20组:count; if(count=20) break;/检查输入的是否合法 /获得各组数据 for(int i=0;icount;i+) cout-endl; cout请输入第i+1di.x; cout请输入第i+1di.y; cout-endl; coutx; while(1) int choice; cout请您选择使用哪种插值法计算：endl; cout (0):退出endl; cout (1):Lagrangeendl; cout (2):Newtonendl; cout (3):Aitkenendl; coutchoice;/取得用户的选择项 if(choice!=0&choice!=1&choice!=2&choice!=3) cout输入错误!endl; if(choice=0)break; if(choice=1)cout你选择了拉格朗日插值计算方法，其结果为：;y=lagrange(x,count);coutx , yendl;/输出最终结果/break;/调用相应的处理函数 if(choice=2)cout你选择了牛顿插值计算方法endl;y=Newton(x,count);cout其结果为：;coutx , yendl;/输出最终结果/break;/调用相应的处理函数 if(choice=3) cout你选择了埃特金插值计算方法，其结果为：;y=aitken(x,count);coutx , yendl;/输出最终结果 6.输入输出结果二、数值积分【开发语言及实现平台或实验环境】C+/C#1.实验目的（1）掌握复化梯形公式求积分的用法，适用范围及精确度。（2）掌握复化辛普生公式求积分的用法，适用范围及精确度。（3）掌握插值型求积分的用法，适用范围及精确度。2.算法思想（1）复化梯形公式将积分区间a,b等分成n个子区间:其中：然后，在每个子区间上使用梯形求积公式计算积分值，再将这些积分值求和就是总的积分值。即：（2）复化辛普生公式将积分区间a,b等分成n个子区间:其中：然后，在每个子区间上使用Simposon求积公式计算积分值，再将这些积分值求和就是总的积分值。即：（3）牛顿-柯特斯公式 对求积节点作适当的限制和选择，可以简化问题的复杂性和提高公式的性能。设将积分区间a,b划分为n等分，步长h=(b-a)/n，选取等距节点x k=a+kh构造出的差值型求积公式 称为牛顿柯特斯公式，式中称为柯特斯系数。按（1.6）式，引进变换x=a+th，则有 , （其中）3.实验内容求函数1/x在区间a,b上的定积分4.定积分算法代码#include #include double f(double x)double y;y=1/x;return y;/-梯形公式-double t2ntn(double a,double b)int n=1;double hn=b-a;double tn=0.0;double t2n=(f(a)+f(b)*hn/2.0;while(fabs(t2n-tn)1e-4)tn=t2n;t2n=0.0;for(int k=0;k=n-1;k+)t2n+=f(a+hn/2.0+k*hn);t2n=tn/2.0+t2n*hn/2.0;n=n*2;hn=hn/2.0;return t2n;/-/-辛普森公式-double sn(double a,double b,int n)/将区间分为n等份的辛普生公式double h=(b-a)/n;double sk=0.0,sk2=0.0;double xk;double xk2;for(int k=1;k=n-1;k+)xk=a+k*h;sk+=f(xk);for(k=0;k=n-1;k+)xk2=a+k*h+0.5*h;sk2+=f(xk2);return h/6.0*(f(a)+f(b)+2.0*sk+4.0*sk2);double s2n(double a,double b,int n)/将区间分为2n等份的辛普生公式double h=(b-a)/(2*n);double sk=0.0,sk2=0.0;double xk;double xk2;for(int k=1;k=(2*n-1);k+)xk=a+k*h;sk+=f(xk);for(k=0;k=(2*n-1);k+)xk2=a+k*h+0.5*h;sk2+=f(xk2);return h/6.0*(f(a)+f(b)+2.0*sk+4.0*sk2);double s2nsn(double m,double n)/for(int c=1;c+;)if(fabs(s2n(m,n,c)-sn(m,n,c)1e-5)break;return s2n(m,n,c);/-/-柯特斯公式-/设置全局数组牛顿 科特斯公式系数表 double C78=1.0/2,1.0/2,1.0/6,4.0/6,1.0/6,1.0/8,3.0/8,3.0/8,1.0/8, 7.0/90,16.0/45,2.0/15,16.0/45,7.0/90, 19.0/288,25.0/96,25.0/144,25.0/144,25.0/96,19.0/188, 41.0/840,9.0/35,9.0/280,34.0/105,9.0/280,9.0/35,41.0/840,751.0/17280,3577.0/17280,1323.0/17280,2989.0/17280,2989.0/17280,1323.0/17280,3577.0/17280,751.0/17280; double cotes1(double a,double b,int n)double cn=0.0;for(int j=0;j=n;j+) cn=cn+Cn-1j/(j*(b-a)/n)+a);return cn;double cotes2(double a,double b,int n)double cn=0.0;for(int j=0;j=n+1;j+) cn=cn+Cnj/(j*(b-a)/(n+1)+a);return cn;double Cotes(double m,double n)while(1)for(int c=1;c+;)if(fabs(cotes2(m,n,c)-cotes1(m,n,c)1e-5)break;return cotes2(m,n,c);/-void main()double a,b;cout求函数1/x的积分endl;cout请输入所求函数的积分区间endl;couta;coutb;while(1)int choice;cout请选择所使用的方法endl;cout 退出endl;cout 复化梯形公式endl;cout复化辛卜生公