qromb.cpp

龙贝格积分

#include "stdafx.h"
#include <math.h>
#include "stdio.h"
#include "Polint.h"
#include "Trapzd.h"

#define EPS 1.0e-6
#define JMAX 200
#define JMAXP JMAX+1
#define K 5

double qromb(double (*func)(double),double a,double b)
{
	double ss,dss;
	double s[JMAXP+1],h[JMAXP+1];
	int j;
	
	h[1]=1.0;
	for (j=1;j<=JMAX;j++) {
		s[j]=trapzd(func,a,b,j);
		if (j >= K) {
			polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss);
			if (fabs(dss) < EPS*fabs(ss)) return ss;
		}
		s[j+1]=s[j];
		h[j+1]=0.25*h[j];
	}
	printf("Too many steps in routine QROMB");
	return 0;
}

#undef EPS
#undef JMAX
#undef JMAXP
#undef K