龙格贝积分算法.c

在C语言环境下的一个龙格贝积分算法例子

//龙格贝积分算法//
//对于不同f(x),修改程序#define f(x)项,本例f(x)=sin(x),区间[a,b]=[1,2],初始h=(b-a)/n=(b-a)/20
#include "stdio.h"
#include "math.h"
#define f(x) (sin(x))
#define N_H 20
#define MAXREPT 10
#define a 1.0
#define b 2.0
#define epsilon 0.00001
double computeT(double aa, double bb, long int n) //复化梯形公式
{
	int i; double sum,h=(bb-aa)/n;sum=0;
	for(i=1;i<n;i++)
		sum+=f(aa+i*h);
	sum+=(f(aa)+f(bb))/2;
	return (h*sum);
}
void main()
{
	int i;
	long int n=N_H,m=0;
	double T[2][MAXREPT+1];
	T[1][0]=computeT(a,b,n);
	n*=2;
	for(m=1;m<MAXREPT;m++)
	{
		for(i=0;i<m;i++)
		{T[0][i]=T[1][i];}
		T[1][0]=computeT(a,b,n);
		n*=2;

		for(i=1;i<=m;i++)
			T[1][i]=T[1][i-1]+(T[1][i-1]-T[0][i-1])/(pow(2,2*m)-1);
		if((T[0][m-1] -T[1][m])<epsilon)
		{
			printf("T=%lf\n",T[1][m]);
			return;
		}


	}
	printf("Return no solved...\n");
}
//-----------End of File------------