longbeige.cpp

数值计算方法中求积分的程序代码

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "conio.h"
double s(double z)
{
	double rr;
	if((z>=0)&&(z<1900))
	{
		rr=exp(log(0.5)-0.523*z/1000);
	}
	else if((z>=1900)&&(z<2100))
	{
		rr=0.185;
	}
	else
	{
		rr=exp(log(0.41)-0.379*z/1000);
	}
	return rr;

}

double m(double z)
{
	double rr;
	if((z>=0)&&(z<1900))
	{
		rr=exp(log(0.55)-0.603*z/1000);
	}
	else if((z>=1900)&&(z<2100))
	{
		rr=0.175;
	}
	else
	{
		rr=exp(log(0.416)-0.412*z/1000);
	}
	return rr;

}

double LongBeiGe(double a,double b,double min,int kmin,int kmax,double (* fx)())
{
	double T[4][2],temp0=0,temp1=0,temp2=0,temp3=0;
	T[0][0]=(b-a)*((*fx)(a)+(*fx)(b))/2;
	T[0][1]=T[0][0]/2+(b-a)/2*(*fx)((a+b)/2);
	T[1][0]=(T[0][1]*4-T[0][0])/3;
	temp0=T[0][1]/2+(b-a)/4*((*fx)((3*a+b)/4)+(*fx)((a+3*b)/4));
	T[0][0]=T[0][1];
	T[0][1]=temp0;
	T[1][1]=(T[0][1]*4-T[0][0])/3;
	T[2][0]=(T[1][1]*16-T[1][0])/15;
	temp0=T[0][1]/2+(b-a)/8*((*fx)(a+(b-a)/8)+(*fx)(a+3*(b-a)/8)+(*fx)(a+5*(b-a)/8)+(*fx)(a+7*(b-a)/8));
    T[0][0]=T[0][1];
	T[0][1]=temp0;
    temp1=(T[0][1]*4-T[0][0])/3;
    T[1][0]=T[1][1];
	T[1][1]=temp1;
    T[2][1]=(T[1][1]*16-T[1][0])/15;
	T[3][0]=(T[2][1]*64-T[2][0])/63;
	T[3][1]=T[3][0];
	double temp=0;
	for(int l=3;l<kmax;l++)
	{
		if(kmin<4)
		{
		    printf("kmin错误，最少分解次数至少为4\n");
			exit(-1);
		}
		temp=0;
		for(int i=1;i<=pow(2,l-1);i++)
		{
			temp+=(*fx)(a+(2*i-1)*(b-a)/pow(2,l));
		}
		temp0=T[0][1]/2+(b-a)/pow(2,l)*temp;
		T[0][0]=T[0][1];
	    T[0][1]=temp0;
		temp1=(T[0][1]*4-T[0][0])/3;
        T[1][0]=T[1][1];
	    T[1][1]=temp1;
		temp2=(T[1][1]*16-T[1][0])/15;
        T[2][0]=T[2][1];
	    T[2][1]=temp2;
		temp3=(T[2][1]*64-T[2][0])/63;
        T[3][0]=T[3][1];
	    T[3][1]=temp3;
		if((fabs(T[3][1]-T[3][0])<min)&&(l>=kmin))  break;
	}
    if(l>=kmax) 
	{
		printf("已达最大分解次数，但还没有满足精度要求\n");
		exit(-1);
	}
	return T[3][1];
}

double hs(double z1,double z2,float ps,float pm,float pl)
{
	double hs=0,ss=0,mm=0,ll=0;
	ss=LongBeiGe(z1,z2,1.0e-6,4,1000,s);
	mm=LongBeiGe(z1,z2,1.0e-6,4,1000,m);
	ll=mm;
	hs=z2-z1-(ps*ss+pm*mm+pl*ll);
	return hs;

}
double z2(double z1,double hs,float ps,float pm,float pl)
{
	double ss=0,mm=0,ll=0,z20=0,z21=0;
	double z2temp=0;
	z20=z1+hs;
	for(register int i=0;i<=1000;i++)
	{
		ss=LongBeiGe(z1,z20,1.0e-6,4,1000,s);
	    mm=LongBeiGe(z1,z20,1.0e-6,4,1000,m);
	    ll=mm;
		z21=hs+z1+(ps*ss+pm*mm+pl*ll);
		if(fabs(z21-z20)<0.001)
			break;
		z20=z21;
	}
	return z21;
}

typedef struct DiCeng
{
	char code[7];
	float depth;
	float time;
	float ps;
	float pm;
	float pl;
	float strip;
	float striptime;
	double hs;
}DiCeng;

void main()
{


}