求NewtonCotes系数

求NewtonCotes系数

1. #include<iostream.h>
2. #include<math.h>
3. # define Precision 0.01//积分精度要求
4. #include <iomanip.h>
5. #include <windows.h>
6. #define jieshu  10   //NewtonCotes的阶数
9. double Factorial (long n)
10. {
11.     long s;
12.     s=1.0;
13.     while (n>0)
14.     {
15.         s=s*n;
16.         n--;
17.     }
19.     return s;
20. }
21. double function(int n,int k,double t)//被积函数
22. {
23.     double s,d;
24.     s=1.0;
25.     for(int j=0;j<=n;j++)
26.     {
27.         if(j!=k)
28.         {
29.             s=s*(t-j);
30.         }
31.         
32.     }
34.     return s;
35. }
36. double VariableStepSimpson(int n1,int k1,double a,double b,double f(int n1,int k1,double x))//变步长Simpson公式
37. {
38.     int n,k,n2,k2;
39.     n2=n1;
40.     k2=k1;
41.     double h,fa,fb,Tn,T2n,s,s1,s2,xk,p;
42.     n=1;
43.     h=b-a;
44.     fa=f(n2,k2,a);
45.     fb=f(n2,k2,b);
46.     Tn=h*(fa+fb)/2.0;//计算T`1`=(b-a)/[f(a)+f(b)];
47.     s1=Tn;
48.     p=Precision+1.0;
49.     while(p>=Precision)
50.     {
51.         s=0.0;
52.         for(k=0;k<n;k++)
53.         {
54.            xk=a+(k+0.501)*h;
56.            s=s+f(n2,k2,xk);
57.         }
58.         T2n=(Tn+h*s)/2.0;//计算T2n
59.         s2=(4.0*T2n-Tn)/3.0;//计算Sn
60.         p=fabs(s2-s1);
61.         Tn=T2n;
62.         s1=s2;
63.         n=n+n;
64.         h=h/2.0;
65.     }
66.     return s2;
67. }
69. main()
70. {
71.   HANDLE hOut;
72.   COORD Position;
73.   hOut = GetStdHandle(STD_OUTPUT_HANDLE);
74.   Position.X = 500;
75.   Position.Y = 500;
76.   SetConsoleScreenBufferSize(hOut,Position);//设置输出屏幕的大小
78.     double a[jieshu+1][jieshu+1];
79.     int n,k;
80.     double s,s1,t;
82.     for(n=1;n<=jieshu;n++)
83.     {
84.         for(k=0;k<=n;k++)
85.         {
86.             a[n][k]=10;
87.         }
88.     }                       //初始化数组
90.       
92.     for(n=1;n<=jieshu;n++)
93.     {
94.         for(k=0;k<=n;k++)
95.         {
96.             s=VariableStepSimpson(n,k,0,n,function);
97.             s=s*pow(-1,n-k);
98.             s1=(n*Factorial(k)*Factorial(n-k));
99.             s=s/s1;
100.             a[n][k]=s;
101.         }                                              //计算系数
102.     } 
105.     for(n=1;n<jieshu;n++)
106.     {
107.         for(k=0;k<=n/2;k++)
108.         {
109.             t=(a[n][k]+a[n][n-k])/2;
110.             a[n][k]=a[n][n-k]=t;    
111.         }
112.     }                                                 //利用对称相等的性质把对称的两个数相加除以2
115.    for(n=1;n<=jieshu;n++)
116.     {
117.         if(n<10)
118.         {
119.             cout<<"n="<<n<<"   ";
120.         }
121.         else
122.         {
123.             cout<<"n="<<n<<"  ";
124.         }
125.     
126.             for(k=0;k<=n;k++)
127.         {
128.             cout<<setiosflags(ios::fixed)<<setprecision(6)<<a[n][k]<<"  ";
129.         }
130.         cout<<endl;
131.     }
132.     
134.     
135. }

xk=a+(k+0.501)*h
#define jieshu  13    //NewtonCotes的阶数
# define Precision 0.001//积分精度要求
程序的运行时间特别长，主要和以上三个因素有关，在xk=a+(k+0.5001)*h中，本应该是乘以0.5的，但是这样运算时，在double function(int n,int k,double t)函数中会得到0，从而使部分运算结果为0（错误），所以就改为0.501了，阶数和积分精度也花了很长运行时间，根据运算的结果，我觉得课本上好像是正确的，-4540/28350=-0.16014109347442680776014109347443 ，-45440/28350= -1.6028218694885361552028218694885,
由于取了0.5001，最后算出的结果Cotes系数不是很对称，最后把对应的两个数相加除以二，这样得到的结果就更加精确，和真实值更加接近

在求积分的时候用了变步长Simpson公式

下面是两种情况下运行的结果，第二种花费时间长一些