dfuntcion.cpp

数值计算中

// DFuntcion.cpp: implementation of the CDFuntcion class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "d.h"
#include "DFuntcion.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CDFuntcion::CDFuntcion()//构造函数
{

}

CDFuntcion::~CDFuntcion()//析构函数
{

}

double CDFuntcion::Func(double x)//运行主函数
{
		return sin(x)/x;
}

double CDFuntcion::Guass(double a,double b,double eps)//高斯迭代方法
{
       int n,k,temp=0;
	   double t1,t2;
       double s1,s2,ep,p,x,fa,fb,h;
	   eps=0.000001;
       n=1;
	   if(a+b==0)
	   { 
		   temp=1;
		   a=0;
	   }
	   h=b-a;
	   fa=a*a;
	   fb=b*b;
	   t1=h*(fa+fb)/2.0;
	   s1=t1;
	   ep=eps+1;
	   while(ep>=eps)
	   {
		   p=0;
		   for(k=0;k<=n-1;k++)
		   {
			   x=a+(k+0.5)*h;
               p=p+sin(x)/x;
		   }
		   t2=(t1+h*p)/2.0;
		   s2=(4*t2-t1)/3.0;
		   ep=fabs(s2-s1);
		   t1=t2;s1=s2;n=n+n;h=h/2;
       }
	   if(temp==1)
		   s1=2*s1;
	return(s1);
}

double CDFuntcion::tixing(double a,double b,double eps)//梯形迭代方法
{
	int i,j=0,n=1,temp=0;
	double c,Tn1=0,Tn2=0;
    double fa,fb,h,sum=0;
	eps=0.000001;
	c=0.001;
     if(a+b==0)
	   { 
		   temp=1;
		   a=0;
	   }
	fa=sin(a)/(c+a);
  	fb=sin(b)/b;
	if(b==0)
		fb=1;
	do
	{
		Tn1=Tn2;
		h=(b-a)/n;
		for(i=1;i<n;i++)
		{
			sum=sum+sin(a+i*h)/(a+i*h);
		}
			Tn2=(h/2)*(fa+fb+2*sum);
			n=2*n;
			sum=0;
			j=j+1;
			
	}
	while(fabs(Tn2-Tn1)>eps);
	if(temp==1)
		Tn2=2*Tn2;
	return(Tn2);

}

double CDFuntcion::Simpson(double a,double b,double eps)//辛普森迭代方法
{
	   int m,i,j,temp=0;
       double s,ep,p,x,aa,bb,h,w,g;
	   static double t[5]={-0.9061798459,-0.538493101,0.0,0.538493101,
		   0.9061798459};
	   static double c[5]={0.2369268851,0.4786286705,0.5688888889,
		   0.4786286705,0.2369268851};
	   eps=0.000001;
	   m=1;
       if(a+b==0)
	   { 
		   temp=1;
		   a=0;
	   }
	   h=b-a;
	   s=fabs(0.001*h);
	   p=1.0e+35;
	   ep=eps+1.0;
	   while((ep>=eps)&&(fabs(h)>s))
	   {
		   g=0;
		   for(i=1;i<=m;i++)
		   {
			   aa=a+(i-1)*h;
			   bb=a+i*h;
			   w=0;

			   for(j=0;j<=4;j++)
			   {
				   x=((bb-aa)*t[j]+(bb+aa))/2.0;
				   w=w+sin(x)/x*c[j];
			   }
			   g=g+w;
		   }
		   g=g*h/2.0;
		   ep=fabs(g-p)/(1.0+fabs(g));
		   p=g;
		   m=m+1;
		   h=(b-a)/m;
	   }
	   if(temp==1)
		   g=2*g;
	return (g);

}