remberg.cpp

龙贝格求积公式的具体实现

// Remberg.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include<math.h>

float f(float x)             //被积函数
{
	float a1;
	a1=x*x*x;
	a1=sqrt(a1);
	return a1;
}

float squ(int y)             //求4的y次方函数
{
	int b;
	float a2=1;
	for(b=1;b<=y;b++)
		a2=a2*4;
	return a2;
}


void main()
{
	float T[20][20];         //存放T表 
	int i,j,k=0,a3;
	float h,a,b,d,m,n;
	printf("*******************************************************************\n");
	printf("*                         龙贝格求积公式                          *\n");
	printf("*******************************************************************\n");
	printf("输入积分下限:");
	scanf("%f",&a);
	printf("输入积分上限:");
	scanf("%f",&b);
	printf("输入精度:");
	scanf("%f",&d);
	h=b-a;
	T[0][0]=h/2*(f(a)+f(b));

	for(i=1;;i++)
	{
		m=1.0;n=0;
		while(m<=2/h-1)
		{
			n+=f(m*h/2);
			m=m+2.0;
		}
		T[0][i]=0.5*T[0][i-1]+h/2*n;
		for(a3=1;a3<=i;a3++)T[a3][i]=(1+1/(squ(a3)-1))*T[a3-1][i]-(1/(squ(a3)-1))*T[a3-1][i-1];
		if(T[i-1][i-1]-T[i][i]<=d)break;
		else h=h/2;
	}
	printf("\n-------------------------------------------------------------------\n");
	printf(" k  ");
	for(j=0;j<i;j++)
		printf("T%d          ",j);
	printf("\n-------------------------------------------------------------------\n");
	for(j=0;j<i;j++)
	{
		printf(" %d  ",j);
		for(k=0;k<=j;k++)
			printf("%f   ",T[k][j]);
		printf("\n");
	}
	printf("-------------------------------------------------------------------\n");

	printf("积分I=%f\n\n",T[i-1][i-1]);
	printf("    ╲︿⌒.︵︵,_.\n");
	printf("    （          ）\n");
	printf("     ︶~︶ˋ~︺︶ ╲  .．...... BY 云\n\n\n");
}