111.cpp

费米函数的拟合算法

#include<iomanip.h>
#include<iostream.h>
#include<fstream.h>
#include<math.h>




double cmp(double x,double y,double p )              // xy分别为往前走和往后走后所得的结果   p为步长

{
	if(x<y) 
		p=-p;
	else p=p;
	return p;
}

double add(double d, double a, double b, double c,double m[],double n[])                    //累加求和
{
	double t=0;
	int i=0;
	while(m[i]!=0)
	{
		t+=(d+a/(exp ((m[i]-b)/c)+1)-n[i])*(d+a/(exp ((m[i]-b)/c)+1)-n[i]);
		i++;
	}
	return t;
}



double AA(double D,double a,double b,double c,double p,double m[],double n[],double Ta,double Tb)                  //递归求出适当的系数值和步长p的值
{
	                                                                    // Tb用来记录走步之后的拟合度 再把它同原来的拟合度进行比较
	D=D+p;
	Tb=add(D,a,b,c,m,n);
	if(Ta>=Tb)
	{
		Ta=Tb;
		return AA(D,a,b,c,p,m,n,Ta,Tb);
	
	}
	else return D;

	
}

double BB(double D,double a,double b,double c,double p,double m[],double n[],double Ta,double Tb)                  //递归求出适当的系数值和步长p的值
{
	                                                                    // Tb用来记录走步之后的拟合度 再把它同原来的拟合度进行比较
	a=a+p;
	Tb=add(D,a,b,c,m,n);
	if(Ta>=Tb)
	{
		Ta=Tb;
		return BB(D,a,b,c,p,m,n,Ta,Tb);
	
	}
	else return a;

	
}

double CC(double D,double a,double b,double c,double p,double m[],double n[],double Ta,double Tb)                  //递归求出适当的系数值和步长p的值
{
	                                                                    // Tb用来记录走步之后的拟合度 再把它同原来的拟合度进行比较
	b=b+p;
	Tb=add(D,a,b,c,m,n);
	if(Ta>=Tb)
	{
		Ta=Tb;
		return CC(D,a,b,c,p,m,n,Ta,Tb);
	
	}
	else return b;

	
}

double DD(double D,double a,double b,double c,double p,double m[],double n[],double Ta,double Tb)                  //递归求出适当的系数值和步长p的值
{
	                                                                    // Tb用来记录走步之后的拟合度 再把它同原来的拟合度进行比较
	c=c+p;
	Tb=add(D,a,b,c,m,n);
	if(Ta>=Tb)
	{
		Ta=Tb;
		return DD(D,a,b,c,p,m,n,Ta,Tb);
	
	}
	else return c;

	
}




double climbD(double D,double a,double b,double c,double p,double m[],double n[])          //此处固定abc来确定系数D的值
{
	double Ta ;double Tb ;double Tc;
	for(int i=0;i<=30;i++)
	{
		Tb=add(D-p, a ,b,c,m,n);
		Tc=add(D+p,a,b,c,m,n);
		if(Ta<Tb&&Ta<Tc)Ta=Ta;
		else
		{
			p=cmp(Tb,Tc,p);
			if(p>0)  Ta=Tc;
			else Ta=Tb;
		}
		D=AA(D,a,b,c,p,m,n,Ta,0);
		p=p/2;
	}
	
	return D;
}


double climba(double D,double a,double b,double c,double p,double m[],double n[])             // 此处用来确定系数a的值
{
	double Ta ;double Tb;double Tc;
		
	for(int i=0;i<=30;i++)
	{
		Tb=add(D,a-p ,b,c,m,n);
		Tc=add(D,a+p,b,c,m,n);
		if(Ta<Tb&&Ta<Tc)Ta=Ta;
		else
		{
			p=cmp(Tb,Tc,p);
			if(p>0)  Ta=Tc;
			else Ta=Tb;
		}
		a=BB(D,a,b,c,p,m,n,Ta,0);
		p=p/2;
	}
	return a;
}

double climbb(double D,double a,double b,double c,double p,double m[],double n[])              //此处用来确定系数b的值
{
	double Ta ;double Tb ;double Tc;
		for(int i=0;i<=50;i++)
		{
			Tb=add(D, a-p ,b,c,m,n);
			Tc=add(D,a+p,b,c,m,n);
			if(Ta<Tb&&Ta<Tc)Ta=Ta;
			else
			{
				p=cmp(Tb,Tc,p);
				if(p>0)  Ta=Tc;
				else Ta=Tb;
			}
			b=CC(D,a,b,c,p,m,n,Ta,0);
			p=p/2;
		}
		return b;
}

double climbc(double D,double a,double b,double c,double p,double m[],double n[])              //此处用来确定系数c的值
{
	double Ta;double Tb ;double Tc;
		for(int i=0;i<=50;i++)
		{
			Tb=add(D, a-p ,b,c,m,n);
			Tc=add(D,a+p,b,c,m,n);
			if(Ta<Tb&&Ta<Tc)Ta=Ta;
			else
			{
				p=cmp(Tb,Tc,p);
				if(p>0)  Ta=Tc;
				else Ta=Tb;
			}
			c=DD(D,a,b,c,p,m,n,Ta,0);
			p=p/2;
		}
		
		return c;
}


void main()
{
	ifstream f1("x.txt");
	double e[29*80+1];
	int n=0;
	while(f1>>e[n])
	{
		n++;
	}
 ifstream f2("y.txt");
	double f[29*80];
	int m=0;
	while(f2>>f[m])
	{
		m++;
	}
	double D,a,b,c,p;
	D=f[0];a=0.75*(f[29*80-1]-f[0]);b=1.68965517241377;
	c=-0.15;p=0.75;

	double x,y;
	x=add(D,a+0.1,b,c,e,f);
	y=add(D,a-0.1,b,c,e,f);
	cout<<setprecision(20)<<x<<' ';
	cout<<setprecision(20)<<y<<' ';

	D=climbD( D, a, b, c, p,e,f) ;
	cout<<setprecision(20)<<D<<endl;   

	p=0.1;


	a=climba( D, a, b, c, p,e,f) ;
	cout<<setprecision(20)<<a<<endl; 
	

	p=0.01;


	c=climbc( D, a, b, c, p,e,f) ;
	cout<<setprecision(20)<<c<<endl; 
	
}