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 }
原文链接: https://www.cnblogs.com/Arthurian/p/7886011.html
欢迎关注
微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍
原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/264035
非原创文章文中已经注明原地址,如有侵权,联系删除
关注公众号【高性能架构探索】,第一时间获取最新文章
转载文章受原作者版权保护。转载请注明原作者出处!