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