zoukankan      html  css  js  c++  java
  • 数值积分之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@aliyun.com

  • 相关阅读:
    Overloaded的方法是否可以改变返回值的类型
    parseXXX的用法
    java的类型转换问题。int a = 123456;short b = (short)a;System.out.println(b);为什么结果是-7616?
    UVA 10405 Longest Common Subsequence(简单DP)
    POJ 1001 Exponentiation(大数处理)
    POJ 2318 TOYS(计算几何)(二分)
    POJ 1265 Area (计算几何)(Pick定理)
    POJ 3371 Flesch Reading Ease (模拟题)
    POJ 3687 Labeling Balls(拓扑序列)
    POJ 1094 Sorting It All Out(拓扑序列)
  • 原文地址:https://www.cnblogs.com/Arthurian/p/7886011.html
Copyright © 2011-2022 走看看