龙贝格积分计算定积分速度很快,和复化梯形积分方法做过对比,不同的被积函数f(x)效率有所不同,但是至少龙贝格积分的速度会是复化梯形的2倍,多的时候有时候有几十倍。
另外网上大部分代码为里查德森外推,虽然理查德外推和龙贝格很相似,但是龙贝格算法的过程是
复化梯形积分->辛普森积分->牛顿科特斯积分->龙贝格积分。
只比复化梯形积分多了3步计算,里查德森外推会多出k步。
代码:
#include <stdio.h>
#include <algorithm>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include <iomanip>
#define EPS 1e-8
#define PI acos(-1.0)
using namespace std;
//使用龙贝格积分计算定积分
//采用C++11的新特性lambda函数作为构造函数传值,这样可以很方便的把f(x)的形式做替换
//只有带C++11的编译器才能正常编译输出
//如果不是C++11,还是老老实实的写一个double f(double x){...}或者写一个class也行
class Calculus{
public:
Calculus(function<double(double x)> &fun,double a,double b,double e);
double CalcIntegration();
private:
function<double(double x)> f;
double low,high,eps;
};
Calculus::Calculus(function<double(double x)> &fun,double a,double b,double e)
{
this->f=fun;
this->low=a;
this->high=b;
this->eps=e;
}
//网上大部分的龙贝格积分代码都是错的,龙贝格积分只应该比复化梯形积分多3次计算
//多k次计算的那个叫里查德森外推
//T-复化梯形积分,S-辛普森积分,C-牛顿科特斯积分,R-龙贝格积分
double Calculus::CalcIntegration()
{
double T1,T2,S1,S2,C1,C2,R1,R2,h;
int n=1,k,m=1;
h=high-low;
T1=0.5*h*(f(high)+f(low));
while(m<=10)
{
double sum=0.0;
for(k=0;k<n;k++)
{
double x=low+(k+0.5)*h;
sum=sum+f(x);
}
n=n*2,h=h*0.5;
T2=0.5*(T1+2*h*sum),S1=(4*T2-T1)/3;
if(m==1)
{
T1=T2,++m;
continue;
}
S2=(4*T2-T1)/3,C1=(16*S2-S1)/15;
if(m==2)
{
S1=S2,T1=T2,++m;
continue;
}
S2=(4*T2-T1)/3,C2=(16*S2-S1)/15;
if(m==3)
{
R1=(64*C2-C1)/63;
C1=C2,S1=S2,T1=T2,++m;
continue;
}
if(m>=4)
{
R2=(64*C2-C1)/63;
cout<<setprecision(16)<<R1<<" "<<R2<<endl;
if(fabs(R2-R1)<eps)
break;
R1=R2,C1=C2,S1=S2,T1=T2,++m;
}
}
cout<<m<<endl;
return R2;
}
int main()
{
//用积分算圆面积,和用PI*R^2公式计算,做对比
function<double(double x)>f=[](double x){ return sqrt(10*x-x*x); };
Calculus romberg(f,0,10,EPS);
double ans=romberg.CalcIntegration();
cout<<setprecision(16)<<2*ans<<" "<<25*PI<<endl;
return 0;
}