romberg.cpp

数值算法 龙贝格积分 C++实现 用迭代的方法求定积分

/* Copyright -- HIT Forum 
 * Romberg integral 
 */
#include <cstdio>
#include <vector>
#include <cmath>
using namespace std;
const double EPS = 1e-10;

double f(double s){
}

double Romberg(double a, double b) 
{ 
   vector<double> R; 
   int k = -1; 
   double r = 0.5 * (b - a) * (f(a) + f(b)); 
   R.push_back(r); 
   do 
   { 
      k += 1; 
      r = 0.0; 
      for (int i = 0; i < pow(2.0, k); i++) 
      { 
         r += f(a + (b - a) * (i + 0.5) / pow(2.0, k)); 
      } 
      r *= (b - a) / pow(2.0, k + 1); 
      r += 0.5 * R[k]; 
      R.push_back(r); 
      for (int m = 0; m <= k; m++) 
      { 
         R[k - m] = (pow(4.0, m + 1) * R[k + 1 - m] - R[k - m]) / 
            (pow(4.0, m + 1) - 1); 
      } 
   } while (fabs(R[0] - R[1]) > EPS); 
   return R[0]; 
} 

/*
f() 为要积分的函数
EPS 为控制的精度
a,b 分别为下界，和上界
*/