qsimp.c

来自「一些常用数值算法的 C程序源代码 对初级学数值计算很有帮助」· C语言 代码 · 共 26 行

#include <math.h>
#define EPS 1.0e-6
#define JMAX 20

float qsimp(float (*func)(float), float a, float b)
{
	float trapzd(float (*func)(float), float a, float b, int n);
	void jswlerror(char error_text[]);
	int j;
	float s,st,ost,os;

	ost = os = -1.0e30;
	for (j=1;j<=JMAX;j++) {
		st=trapzd(func,a,b,j);
		s=(4.0*st-ost)/3.0;
		if (fabs(s-os) < EPS*fabs(os)) return s;
		if (s == 0.0 && os == 0.0 && j > 6) return s;
		os=s;
		ost=st;
	}
	jswlerror("Too many steps in routine qsimp");
	return 0.0;
}
#undef EPS
#undef JMAX