计算方法上机实习报告5 - 华中科技大学

MACROBUTTONMTEditEquationSection2SEQMTEqn\r\hSEQMTSec\r1\hSEQMTChap\r1\h计算方法上机实习报告5提出问题给出概率积分的数据表：0.460.470.480.490.48465550.49375420.50274980.5116683试用二次插值计算。已知的函数表1.51.61.70.997490.999570.99166试构造出差商表，利用二次Newton插值公式计算sin1.609（保留五位有效数字）。求不高于4次的多项式，使它满足并写出余项表达式。用最小二乘法求一个形如的经验公式，使与下列数据相拟合x1925313841y19.032.349.073.397.8分析问题（1）题目给出概率积分的四个插值节点，要求用二次插值计算，而我们知道二次插值只需要三个插值节点即可，在该题中我们尝试取前三个点作为插值节点。 该题的重点在求插值基函数而它们由公式 给出。而给出插值公式。（2）Newton基本插值公式 它的各项系数就是函数的各阶差商，每增加一个插值节点，只需要在原来的基础上多计算一项，这一性质被称作承袭性。（3）该题需要确定一个4次插值多项式，就是要获得其各项系数的解。4次多项式有5个系数，而题目正好给予了5个条件，这样我们会获得关于系数的非奇异五元一次方程组，在运用第三章的选列主元消元法求出各系数。（4）本题要求用最小二乘法求经验公式。实际上也就是确定拟合曲线的各项系数，关于系数的线性方程组如下：三.解决问题 (1)C语言代码如下：#include"stdio.h"#include"iostream.h"#include"math.h"#defineN3voidGet_l(intk,doublex,doublexk[N],doubleL[N]);voidmain(){ inti; doublel[N],L[N],x,Lx=0; doublexk[N]={0.46,0.47,0.48}; doublef[N]={0.4846555,0.4937542,0.5027498}; printf("输入需插值节点x:\n"); scanf("%lf",&x); for(i=0;i<N;i++) { Get_l(i,x,xk,l); L[i]=l[i]*f[i]; Lx+=L[i]; } cout.precision(7); cout<<"f(x)="<<Lx<<endl;}/*获得基函数lk(x)*/voidGet_l(inti,doublex,doublexk[N],doublel[N]){ intk; doubley1=1,y2=1; for(k=0;k<N;k++) if(k!=i) { y1*=(x-xk[k]); y2*=(xk[i]-xk[k]); } l[i]=y1/y2;}结果如下：（2）C语言代码如下:#include"stdio.h"#include"math.h"#defineN3#defineM4voidPut_out(inti,doublea[N][M]);doubleGet_value(doublex,doublea[N][M]);voidmain(){ inti,j; doublea[N][M],sinx,x1; doublex[N]={1.5,1.6,1.7}; doublef[N]={0.99749,0.99957,0.99166}; for(i=0;i<N;i++) { a[i][0]=x[i]; a[i][1]=f[i]; } for(i=0;i<N;i++) Put_out(i,a); printf("k\tf(xk)\t\tf(x0,xk)\tf(x0,x1,xk)\tf(x0,x1,x2,xk)\n"); for(i=0;i<N;i++) { printf("%d\t",i); for(j=0;j<i+2;j++) { printf("%lf\t",a[i][j]); } printf("\n"); } printf("Putinx:\n"); scanf("%lf",&x1); sinx=Get_value(x1,a); printf("sin(%lf)=%lf\n",x1,sinx);}/*输出差商表*/voidPut_out(inti,doublea[N][M]){ intj; for(j=2;j<M;j++) { if(j<i) a[i][j]=(a[i][j-1]-a[0][j-1])/(a[i][0]-a[j-1][0]); if(j<i+2) a[i][j]=(a[i][j-1]-a[j-2][j-1])/(a[i][0]-a[j-2][0]); else a[i][j]=0; }}以上代码仅是主函数和输出差商部分，求值用的是选列主元消元法，前面有所涉及，在这里不再赘述。输出差商表如下：k01.500000.9974911.600000.999570.0208021.700000.99166-0.02915-0.49950运行结果如图：（3）C语言代码如下：#include"stdio.h"#include"math.h"#defineN5voidGet_a(doublec[2][N],double[N][N]);voidfound(doublea[N][N],doubleb[N],intj);voidGauss(doublea[N][N],doubleb[N],intj);voidGet_root(doublea[N][N],doubleb[N],doublex[N]);voidmain(){ intj; doublea[N][N],x[N]; doublec[2][N]={{1,0,1,0,1},{1,1,2,2,3}}; doubleb[N]={-2,4,0,0,2}; Get_a(c,a); for(j=0;j<N-1;j++) { found(a,b,j); /*寻找第j列主元*/ Gauss(a,b,j); /*对第j列进行消元*/ } Get_root(a,b,x); /*求解*/ printf("输出各项系数(升序):\n"); for(j=0;j<N;j++) printf("%lf\t",x[j]); }voidGet_a(doublec[2][N],doublea[N][N]){ inti,j; for(i=0;i<N;i++) { if(c[0][i]==1) for(j=0;j<N;j++) a[i][j]=pow(c[1][i],double(j)); if(c[0][i]==0) for(j=0;j<N;j++) a[i][j]=j*pow(c[1][i],double(j-1)); }}这是主函数和获得关于系数的五元一次方程组部分，解系数采用选列主元消元法。运行结果如下：也就是说，也就是参考答案中的展开式。（4）C语言代码如下：#include"stdio.h"#include"math.h"#defineN5#defineM2voidGet_x(doublex[M][M],doubleflag[M],doublexy[N]);voidGet_d(doubled[M],doubleflag[M],doublexy[N],doubley[N]);voidfound(doublea[M][M],doubleb[M],intj);voidGauss(doublea[M][M],doubleb[M],intj);voidGet_root(doublea[M][M],doubleb[M],doublex[M]);voidmain(){ inti,j; doubleab[M]; doubleflag[M]={0,2}; doublex[M][M],d[M]; doublexy[N]={19,25,31,38,44}; doubley[N]={19,32.3,49,73.3,97.8}; Get_x(x,flag,xy); Get_d(d,flag,xy,y); for(j=0;j<M-1;j++) { found(x,d,j); /*寻找第j列主元*/ Gauss(x,d,j); /*对第j列进行消元*/ } Get_root(x,d,ab); /*求解*/ for(j=0;j<M;j++) { printf("第%d次项系数:%lf\n",int(flag[j]),ab[j]); }}voidGet_x(doublex[M][M],doubleflag[M],doublexy[N]){ inti,j,k; doubletemp; for(i=0;i<M;i++) for(j=0;j<M;j++) { temp=0; if(j<i) x[i][j]=x[j][i]; else { for(k=0;k<N;k++) { temp+=pow(xy[k],flag[i])*pow(xy[k],flag[j]); } x[i][j]=temp; } }}voidGet_d(doubled[M],doubleflag[M]