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@aliyun.com