数值积分之Simpson公式与梯形公式

Simpson（辛普森）公式和梯形公式是求数值积分中很重要的两个公式，可以帮助我们使用计算机求解数值积分，而在使用过程中也有多种方式，比如复合公式和变步长公式。这里分别给出其简单实现（C++版）： 

1、复合公式：

 1 #include<iostream>
 2 #include<iomanip>
 3 #include <cmath>
 4 using namespace std; 
 5 
 6 double Simpson(double a,double b,double n);
 7 double Compound_Trapezoid(double a,double b,double n);
 8 
 9 int main()
10 {
11     int n;
12     double a, b;
13     cout << "区间数n:";
14     cin >> n;
15     cout << "区间端点a:";
16     cin >> a;
17     cout<<"区间端点b:";
18     cin >> b;
19     cout<<setprecision(20)<<Simpson(a,b,n)<<endl;
20     cout<<setprecision(20)<<Compound_Trapezoid(a,b,n)<<endl;
21     getchar();
22     getchar();
23     return 0;
24 }
25 
26 /* 
27  * Simpson算法 
28  */
29 double Simpson(double a,double b,double n)
30 {
31     double h=(b-a)/n;
32     double Sn=exp(a)-exp(b);
33     for (double x=a+h;x<=b;x+=h)
34     {
35         Sn+=4*exp(x-h/2)+2*exp(x);
36     }
37     Sn *= h/6;
38     return Sn;
39 }
40 
41 /*
42  * 复合梯形算法
43  */
44 double Compound_Trapezoid(double a,double b,double n)
45 {
46     double h=(b-a)/n;
47     double Sn=exp(a)+exp(b);
48     for(double x=a+h;x<b;x+=h)
49     {
50         Sn += 2 * exp(x);
51     }
52     Sn *= h/2;
53     return Sn;
54 }

2、变步长公式

  1 /*
  2  * e^x，1/x求1到3的积分
  3  * 精确到1E-5
  4  */
  5 #include<iostream>
  6 #include<iomanip>
  7 #include<cmath>
  8 
  9 using namespace std;
 10 
 11 
 12 //变步长梯形法
 13 double ex_Variable_step_size_trape(double ,double ,double);
 14 double x_Variable_step_size_trape(double ,double ,double);
 15 //变步长Simpson法
 16 double ex_Variable_step_size_Simpson(double ,double ,double);
 17 double x_Variable_step_size_Simpson(double ,double ,double);
 18 //Romberg法
 19 //double Romberg();
 20 
 21 int main()
 22 {
 23     //左端点a，右端点b，允许误差E 
 24     double a,b,E;
 25     cout << "请输入左端点a：";
 26     cin >> a;
 27     cout << "请输右端点b：";
 28     cin >> b;
 29     cout << "请输入允许误差E：";
 30     cin >> E;
 31     cout << "变步长梯形(e^x)：" << setiosflags(ios::fixed) << setprecision(5) << ex_Variable_step_size_trape(a,b,E) << endl;
 32     cout << "变步长Simpson(e^x)：" << setiosflags(ios::fixed) << setprecision(5) << ex_Variable_step_size_Simpson(a,b,E) << endl;
 33     cout << "变步长梯形(1/x)：" << setiosflags(ios::fixed) << setprecision(5) << x_Variable_step_size_trape(a,b,E) << endl;
 34     cout << "变步长Simpson(1/x)：" << setiosflags(ios::fixed) << setprecision(5) << x_Variable_step_size_Simpson(a,b,E) << endl;
 35     getchar();
 36     getchar();
 37     return 0;
 38 } 
 39 
 40 double ex_Variable_step_size_trape(double a,double b,double E)
 41 {
 42        double h = b - a, e = 0 ,T2 = 0;
 43        double T1 = h/2 * (exp(a) + exp(b));
 44        do
 45        {
 46               double S = 0, x = a + h/2;
 47               do 
 48               {
 49                      S += exp(x);
 50                      x += h;
 51               }while(x < b);
 52               T2 = T1/2 + h/2 * S;
 53               e = (T2 > T1)?(T2 - T1):(T1 - T2);
 54               h = h/2;
 55               T1 = T2;
 56        }while(e > E);
 57        return T2;
 58 }
 59 
 60 double x_Variable_step_size_trape(double a,double b,double E)
 61 {
 62        double h = b - a, e = 0 ,T2 = 0;
 63        double T1 = h/2 * (1/a + 1/b);
 64        do
 65        {
 66               double S = 0, x = a + h/2;
 67               do 
 68               {
 69                      S += 1/x;
 70                      x += h;
 71               }while(x < b);
 72               T2 = T1/2 + h/2 * S;
 73               e = (T2 > T1)?(T2 - T1):(T1 - T2);
 74               h = h/2;
 75               T1 = T2;
 76        }while(e > E);
 77        return T2;
 78 }
 79 
 80 
 81 double ex_Variable_step_size_Simpson(double a,double b,double E)
 82 {
 83        double h = b - a, e = 0 ,T2 = 0;
 84        double T1 = h/6 * (exp(a) - exp(b));
 85        do
 86        {
 87               double S = 0, x = a + h/2;
 88               do 
 89               {
 90                      S += 2 * exp(x);
 91                      x += h/2;
 92                      S += 1 * exp(x);
 93                      x += h/2;
 94               }while(x <= b);
 95               T2 = T1/2 + h/6 * S;
 96               e = (T2 > T1)?(T2 - T1):(T1 - T2);
 97               h = h/2;
 98               T1 = T2;
 99        }while(e > E);
100        return T2;
101 }
102 
103 double x_Variable_step_size_Simpson(double a,double b,double E)
104 {
105        double h = b - a, e = 0 ,T2 = 0;
106        double T1 = h/6 * (1/a - 1/b);
107        do
108        {
109               double S = 0, x = a + h/2;
110               do 
111               {
112                      S += 2 * 1/x;
113                      x += h/2;
114                      S += 1 * 1/x;
115                      x += h/2;
116               }while(x <= b);
117               T2 = T1/2 + h/6 * S;
118               e = (T2 > T1)?(T2 - T1):(T1 - T2);
119               h = h/2;
120               T1 = T2;
121        }while(e > E);
122        return T2;
123 }

 

作者：耑新新，发布于  博客园

转载请注明出处，欢迎邮件交流：zhuanxinxin@foxmail.com

转载于:https://www.cnblogs.com/Amedeo/p/7886011.html