不等距节点的表格函数的数值积分问题

实验类别：综合性专 业：弹药工程与爆炸技术班 级：08010643学 号：36姓 名：周渝捷中北大学理学院实验一 不等距节点的表格函数的数值积分问题【实验内容】1.通过画图准确地判断出拟合多项式的次数。 2.熟练掌握最小二乘法。 3.掌握梯形公式，辛普生公式，龙贝格公式 4.了解各种数值积分方法的异同点，并能融会贯通。 5.掌握龙贝格方法。【实验方法与步骤】最小二乘法： #include #include #include #includefloat average(int n,float *x) int i; float av; av=0; for(i=0;in;i+) av+=*(x+i); av=av/n; return(av); /平方和 float spfh(int n,float *x) int i; float a,b; a=0; for(i=0;in;i+) a+=(*(x+i)*(*(x+i); return(a); /和平方 float shpf(int n,float *x) int i; float a,b; a=0; for(i=0;in;i+) a=a+*(x+i); b=a*a/n; return(b); /两数先相乘，再相加 float dcj(int n,float *x,float *y) int i; float a; a=0; for(i=0;in;i+) a+=(*(x+i)*(*(y+i); return(a); /两数先相加，再相乘 float djc(int n,float *x,float *y) int i; float a=0,b=0; for(i=0;in;i+) a=a+*(x+i); b=b+*(y+i); a=a*b/n; return(a); /系数a float xsa(int n,float *x,float *y) float a,b,c,d,e; a=spfh(n,x); b=shpf(n,x); c=dcj(n,x,y); d=djc(n,x,y); e=(c-d)/(a-b); /printf(%f %f %f %f,a,b,c,d); return(e); float he(int n,float *y) int i; float a; a=0; for(i=0;in;i+) a=a+*(y+i); return(a); float xsb(int n,float *x,float *y,float a) float b,c,d; b=he(n,y); c=he(n,x); d=(b-a*c)/n; return(d); void main() int n,i; float *x,*y,a,b; printf(请输入将要输入的有效数值组数n的值); scanf(%d,&n); printf(n=%dn,n); x=(float*)calloc(n,sizeof(float); if(x=NULL) printf(内存分配失败); exit(1); y=(float*)calloc(n,sizeof(float); if(y=NULL) printf(内存分配失败); exit(1); printf(请输入x的值); for(i=0;in;i+)scanf(%f,x+i); printf(请输入y的值,请注意与x的值一一对应); for(i=0;in;i+)scanf(%f,y+i); /x0=0.1;x1=0.3;x2=0.4;x3=0.55;x4=0.7;x5=0.8;x6=0.95; /y0=15;y1=18;y2=19;y3=21;y4=22.6;y5=23.8;y6=26; for(i=0;in;i+) printf(x%d=%3.2f ,i,*(x+i); printf(y%d=%3.2fn,i,*(y+i); a=xsa(n,x,y); b=xsb(n,x,y,a); printf(经最小二乘法拟合得到的一元线性方程为:n); printf(f(x)=%3.2fx+%3.2fn,a,b梯形公式与辛普生公式：#include #includeusing namespace std;int main ()double x8;double y8;int xuan1;cout你是要使用默认数据还是自行输入数据：endl;cout1,默认数据； 2，自行输入endl;cout*xuan1;if(xuan1=1)y0=3.82;y1=4.455;y2=5.09;y3=5.725;y4=6.36;y5=6.995;y6=7.63;y7=8.265;y8=8.9;for(int i=1;i=9;i+)xi-1=1+(i-1.0);elsecout您可以输入9组数值endl;for(int f=1;f=9;f+)cout请输入xf-1xf-1;cout*;cout请输入yf-1yf-1;coutendl;while(1)/程序数据输入完毕 ，操作开始int xuan2;cout请选择计算公式：endl;cout1,梯形公式 2，辛浦生公式endl;cout*xuan2;if(xuan2=1)/选择了梯形公式，计算开始double h;h=(x8-x0)/16;double t;double temp=0;for(int e=1;e=7;e+)temp=temp+ye; t=h/2*(y0+y8+2*temp);cout采用梯形公式计算结果为：setiosflags(ios:fixed)setprecision(8)tendl;cout请选择下一步行动：endl;cout1,继续采用其他的运算方法运算 2，退出endl;cout*xuan3;if(xuan3=1);else break; else/选择了辛浦生计算方法，开始计算double h=(x8-x0)/16;double s2=0;/奇数和，乘以4double s1=0;/偶数和，乘以2for(int r=1;r=7;r=r+2)s2=s2+yr;for(int t=2;t=6;t=t+2)s1=s1+yt;double s;s=h/3*(y0+y8+s2*4+s1*2);cout采用辛浦生公式计算结果为：setiosflags(ios:fixed)setprecision(8)sendl;cout请选择下一步行动：endl;cout1,继续采用其他的运算方法运算 2，退出xuan4;if(xuan4=1);else break;return 0;龙贝格公式：#include#includeusing namespace std;#define f(x) (1.27*x+2.55) /举例函数#define epsilon 0.0001 /精度#define MAXREPT 10 /迭代次数,到最后仍达不到精度要求,则输出T(m=10).double Romberg(double aa, double bb) /aa,bb 积分上下限 int m, n;/m控制迭代次数, 而n控制复化梯形积分的分点数. n=2m double h, x; double s, q; double ep; /精度要求 double *y = new doubleMAXREPT;/为节省空间,只需一维数组 /每次循环依次存储Romberg计算表的每行元素,以供计算下一行,算完后更新 double p;/p总是指示待计算元素的前一个元素(同一行) /迭代初值 h = bb - aa; y0 = h*(f(aa) + f(bb)/2.0; m = 1; n = 1; ep = epsilon + 1.0; /迭代计算 while (ep = epsilon) & (m MAXREPT) /复化积分公式求T2n(Romberg计算表中的第一列),n初始为1,以后倍增 p = 0.0; for (int i=0; in; i+)/求Hn x = aa + (i+0.5)*h; p = p + f(x); p = (y0 + h*p)/2.0;/求T2n = 1/2(Tn+Hn),用p指示 /求第m行元素,根据Romberg计算表本行的前一个元素(p指示), /和上一行左上角元素(yk-1指示)求得. s = 1.0; for (int k=1; k=m; k+) s = 4.0*s; q = (s*p - yk-1)/(s - 1.0); yk-1 = p; p =