s3.cpp

数值分析的实验,可工大家参考的,能通过,实验的代码

#include<iostream>
#include<math.h>
#define E 2.71828
#define ERRORLIMIT 0.000001
#define F f2
#define INIT f2init
#define PI 3.1415927
int k;
double a,b,h;
double *x,*T;
using namespace std;
double f1(double x)
{
	return x*x*exp(-1.*x*x);
}
double f2(double x)
{
	return 1./tan(x);
}
double f3(double x)
{
	return 1./(x*x-1.);
}
double f4(double x)
{
	return pow(x,3./2.);
}
void f1init()
{
	a=0;
	b=2;
}
void f2init()
{
	a=PI/2.;
	b=(3./4.)*PI;
}
void f3init()
{
	a=2;
	b=3;
}
void f4init()
{
	a=0;
	b=1;
}
void Long()
{
	int i;
	double Ef,error(1);
	double *xtemp,*Ttemp;
	INIT();
	k=0;
	h=b-a;
	T=new double;
	T[0]=(h/2.)*(F(a)+F(b));
	x=new double[2];
	x[0]=a;
	x[1]=b;
	cout<<k<<' '<<T[0]<<endl;
	while(error>=ERRORLIMIT)
	{
		k++;
		xtemp=new double[(int)pow(2,k)];
		for(i=0;i<(int)pow(2,k);i++)
			xtemp[i]=x[i];
		x=new double[(int)pow(2,k)+1];
		for(i=0;i<(int)pow(2,k)+1;i+=2)
			x[i]=xtemp[i/2];
		Ef=0;
		for(i=1;i<(int)pow(2,k);i+=2)
		{
			x[i]=(x[i-1]+x[i+1])/2.;
			Ef+=F(x[i]);
		}
		Ttemp=new double[k];
		for(i=0;i<k;i++)
			Ttemp[i]=T[i];
		T=new double[k+1];
		T[0]=Ttemp[0]/2.+(h/2.)*Ef;
		h/=2.;
		for(i=1;i<k+1;i++)
			T[i]=(pow(4.,i)/(pow(4.,i)-1.))*T[i-1]-(1./(pow(4.,i)-1.))*Ttemp[i-1];
		error=fabs(T[k]-Ttemp[k-1]);
		cout<<k<<' ';
		for(i=0;i<k+1;i++)
			cout<<T[i]<<' ';
		cout<<endl;
	}
	cout<<"I="<<T[k]<<endl;
}
int main()
{
	Long();
	return 0;
}