hypercycle.cpp

龙格库塔法求解二维反应扩散方程

/****************************************************
A program for hypercycle with Runge-Kutta. 2008.3.7 
****************************************************/
#define S  7               //the types of species.
#define m  147             //rectangular grid   
#define n  147
#include <stdlib.h> 
#include <conio.h> 
#include <stdio.h> 
#include <math.h>
#include "debug.h"
#define OFF


float X[S][m][n],XX[S][m][n],PXX[S][m][n];


//the parameters for the mothod of Runge-Kutta.
float K1[S][m][n], K2[S][m][n], K3[S][m][n],K4[S][m][n];
float kc[S][S];   //the rate of replication of species i catalyzed by species j.
 

//other parameters for the model
float sigma[S], rou[S],d[S]; 

float f(int sp, float sg, float r, float kk[S], float x[S]);
void extract1(int xi, int yj, float xtp[S]);
void extract2(int sp, float ktp[S]);
void extract3(int sp, int d, int xi, int yj, float xtp[S]);
void increase(int sp, int xi, int yi, float ta,float kx[S],float xtp[S]);
double gaussrand();


int main()
{
	float init=0,fint=200,tao=0.1,h=1,t;
	float xtemp[S],ktemp[S],k1[S],k2[S],k3[S],k4[S],temp;  
	int i,j,k,ss;

	//the file pointer for different species.
	FILE* fp[S];
	
	//the file pointer for array X.
	FILE* fpp;

	const char *file;
	char buff[]={'s','p','0','.','d','a','t','\0'};
	char sp[3]={'0','0','\0'};
	for(i=0;i<S;i++)
	{
		sprintf(sp, "%1d",i);
		buff[2]=sp[0];
		//buff[3]=sp[1];
		file=buff;
		if((fp[i]=fopen(file,"w"))==NULL)
		{
			printf("Creating files error!\n");
			exit(0);
		}
	}
	//initial values
	for(i=0;i<S;i++)
	{
		sigma[i]=1;
		rou[i]=2;
		d[i]=1;
	}

	for(i=0;i<S;i++)
		for(j=0;j<S;j++)
			kc[i][j]=0;
	//kc[1][0]=500; kc[2][1]=500; kc[0][2]=500;
	kc[0][5]=500; kc[1][0]=500; kc[2][1]=500; kc[3][2]=500;
	kc[4][3]=500; kc[5][4]=500; kc[6][5]=0;
	
	//initial array X
	if((fpp=fopen("randx0.dat","r"))!=NULL)
	{
		for(k=0;k<S;k++)
			for(i=0;i<m;i++)
				for(j=0;j<m;j++)
					fscanf(fpp,"%fg",&X[k][i][j]);				
	}
	else
	{
		//randx.dat store the value of the array X.
		fpp=fopen("randx.dat","w");  
		for(k=0;k<S;k++)
			for(i=0;i<m;i++)
			{
				for(j=0;j<m;j++)
				{
					X[k][i][j]=rand()/(RAND_MAX+1.0); //0--1均匀分布
					//X[k][i][j]=(gaussrand()/40.0+0.1)/2; //0--1正态分布
					fprintf(fpp,"%.4f ",X[k][i][j]);
				}
				fprintf(fpp,"\n");
			}
		fclose(fpp);
	}

	//all the concentration at different site should be unity,or normalization.
	for(i=0;i<m;i++)
		for(j=0;j<m;j++)
		{
			temp=0.0;
			for(k=0;k<S;k++)
				temp+=X[k][i][j];
			if(temp>1.0)
				for(k=0;k<S;k++)
					X[k][i][j] /= temp;
		}

	//initial value 
	t=init;
	ss = -1;
	//caculate
	for(;;)
	{
		//fist caculate the term of diffusion,
		//caculate the center.
		for(k=0;k<S;k++)
			for(i=1;i<m-1;i++)
				for(j=1;j<n-1;j++)
					PXX[k][i][j]=X[k][i+1][j]+X[k][i-1][j]+X[k][i][j+1]+ \
					        X[k][i][j-1]-4*X[k][i][j];
		
		//caculate the left and right boundary except two tips.
		for(k=0;k<S;k++)
			for(i=1;i<m-1;i++)
			{
				PXX[k][i][0]=X[k][i+1][0]+X[k][i-1][0]+2*X[k][i][1]- \
					4*X[k][i][0];
				PXX[k][i][n-1]=X[k][i+1][n-1]+X[k][i-1][n-1]+2*X[k][i][n-2]- \
					4*X[k][i][n-1];				
			}
		
		//caculate the top and bottom boundary except two tips.
		for(k=0;k<S;k++)
			for(j=1;j<n-1;j++)
			{
				PXX[k][0][j]=X[k][0][j+1]+X[k][0][j-1]+2*X[k][1][j]- \
					4*X[k][0][j];
				PXX[k][m-1][j]=X[k][m-1][j+1]+X[k][m-1][j-1]+2*X[k][m-2][j]- \
					4*X[k][m-1][j];
			}

		//caculate the four tops.
		for(k=0;k<S;k++)
		{
			PXX[k][0][0]=2*X[k][1][0]+2*X[k][0][1]-4*X[k][0][0];
			PXX[k][0][n-1]=2*X[k][1][n-1]+2*X[k][0][n-2]-4*X[k][0][n-1];
			PXX[k][m-1][0]=2*X[k][m-2][0]+2*X[k][m-1][1]-4*X[k][m-1][0];
			PXX[k][m-1][n-1]=2*X[k][m-2][n-1]+2*X[k][m-1][n-2]-4*X[k][m-1][n-1];
		}

#ifdef ON
		for(k=0;k<S;k++)
			for(i=0;i<m;i++)
			{	for(j=0;j<n;j++)
					TRACE("%.4f ",PXX[k][i][j]);
				TRACE("\n");
			}
#endif
		//caculate the function with Runge-Kutta.
		for(k=0;k<S;k++)
			for(i=0;i<m;i++)
				for(j=0;j<n;j++)
				{

#ifdef ON
					//debug
					TRACE("X[%d][%d][%d]=%.4f\n",k,i,j,X[k][i][j]);
#endif
					extract1(i,j,xtemp); //pick out X[k][i][j] at site(i,j) to xtemp[k];
					
#ifdef ON
					for(int ki=0;ki<S;ki++)//debug
						TRACE("xtemp[%d]=%.4f ",ki,xtemp[ki]);
					TRACE("\n");
#endif

					extract2(k,ktemp);   //pick out kc[S][S] for species k to ktemp[k];
					
#ifdef ON
					for(ki=0 ; ki<S; ki++)//debug
						TRACE("ktemp[%d]=%.4f ",ki,ktemp[ki]);
					TRACE("\n");
#endif
					//caculate the K1[k][i][j] for species k at site(i,j)
					K1[k][i][j]=f(k,sigma[k],rou[k],ktemp,xtemp);
#ifdef ON
					//debug
					TRACE("K1[%d][%d][%d]=%.4f\n",k,i,j,K1[k][i][j]);
#endif

					extract3(k,1,i,j,k1);//pick out K1[k][i][j] at site(i,j) to k1[k];
#ifdef ON					
					//debug
					TRACE("k1[%d]=%.4f\n",k,k1[k]);
#endif
					//caculate X[k][i][j]+k1[k][i][j]*tao/2
					increase(k,i,j,tao,k1,xtemp);
#ifdef ON
					//debug
					TRACE("after increase,xtemp[%d]=%.4f\n",k,xtemp[k]);
#endif
					//caculate K2[k][i][j] for species k at site(i,j)
					K2[k][i][j]=f(k,sigma[k],rou[k],ktemp,xtemp);
#ifdef ON
					//debug
					TRACE("K2[%d][%d][%d]=%.4f\n",k,i,j,K2[k][i][j]);
#endif

					extract3(k,2,i,j,k2);//pick out K2[k][i][j] at site(i,j) to k2[k];
#ifdef ON
					//debug
					TRACE("k2[%d]=%.4f\n",k,k2[k]);
#endif
					//caculate X[k][i][j]+k2[k][i][j]*tao/2
					increase(k,i,j,tao,k2,xtemp);
#ifdef ON
					//debug
					TRACE("after increase,xtemp[%d]=%.4f\n",k,xtemp[k]);
#endif
					//caculate K3[k][i][j] for species k at site(i,j)
					K3[k][i][j]=f(k,sigma[k],rou[k],ktemp,xtemp);
#ifdef ON
					//debug
					TRACE("K3[%d][%d][%d]=%.4f\n",k,i,j,K3[k][i][j]);
#endif
					extract3(k,3,i,j,k3);//pick out K3[k][i][j] at site(i,j) to k3[k];
#ifdef ON
					//debug
					TRACE("k3[%d]=%.4f\n",k,k3[k]);
#endif
					//caculate X[k][i][j]+k3[k][i][j]*tao
					increase(k,i,j,2*tao,k3,xtemp);
#ifdef ON
					//debug
					TRACE("after increase,xtemp[%d]=%.4f\n",k,xtemp[k]);
#endif
					//caculate K4[k][i][j] for species k at site(i,j)
					K4[k][i][j]=f(k,sigma[k],rou[k],ktemp,xtemp);
#ifdef ON
					//debug
					TRACE("K4[%d][%d][%d]=%.4f\n",k,i,j,K4[k][i][j]);
#endif
					XX[k][i][j]=X[k][i][j]+tao*((K1[k][i][j]+2*K2[k][i][j]+2*K3[k][i][j]+ \
						K4[k][i][j])/6+PXX[k][i][j]);
#ifdef ON
					//debug
					TRACE("XX[%d][%d][%d]=%.4f\n",k,i,j,XX[k][i][j]);
#endif
					if(XX[k][i][j]<0) 
						XX[k][i][j]=0.0;
/****************************/
				}
		
		//count the cycle
		ss = ss+1;
		if(ss>=100)
		{
			printf("progress:%.0f/%.0f\r",t,fint);
			ss=0;
		}
		t=t+tao;
		if(t>fint)
		{
			for(k=0;k<S;k++)
				for(i=0;i<m;i++)
				{
					for(j=0;j<n;j++)
					{
						fprintf(fp[k],"%04f ",X[k][i][j]);
					}
					fprintf(fp[k],"\n");
				}
			break;
		}
		for(k=0;k<S;k++)
			for(i=0;i<m;i++)
				for(j=0;j<n;j++)
					X[k][i][j]=XX[k][i][j];
	}
	for(k=0;k<S;k++)
		fclose(fp[k]);
	return 0;
}

void increase(int sp, int xi, int yj, float ta,float kx[],float xtp[S])
{
	xtp[sp]=X[sp][xi][yj]+ta*kx[sp]/2;
	return;
}

float f(int sp,float sg,float r, float k[S],float x[S])
{
	int i,j;
	float result=0,temp1=0,temp2=0;
	result=(-1)*sg*x[sp];
	for(i=0;i<S;i++)
		temp1+=x[i];
	temp1=1-temp1;
	for(j=0;j<S;j++)
	{
		if(j!=sp)
			temp2+=k[j]*x[j];
	}
	temp2=r+temp2;
	result+=temp1*temp2*x[sp];
	return result;
}

void extract1(int xi, int yj, float xtp[S])
{
	int i;
	for(i=0;i<S;i++)
		xtp[i]=X[i][xi][yj];		
	return;
}

void extract2(int sp, float ktp[S])
{
	int i;
	for(i=0;i<S;i++)
		ktp[i]=kc[sp][i];
	return;
}

void extract3(int sp,int d,int xi, int yj, float xtp[S])
{
	switch (d)
	{
		case 1:	xtp[sp]=K1[sp][xi][yj];break;
		case 2:	xtp[sp]=K2[sp][xi][yj];break;
		case 3:	xtp[sp]=K3[sp][xi][yj];break;
		case 4:	xtp[sp]=K4[sp][xi][yj];
	}
	return;
}

double gaussrand()
{
	static double V1, V2, s;
	static int phase = 0;
	double X;
	
	if(phase == 0) {
		do {
			double U1 = (double)rand() / RAND_MAX;
			double U2 = (double)rand() / RAND_MAX;
			
			V1 = 2 * U1 - 1;
			V2 = 2 * U2 - 1;
			s = V1 * V1 + V2 * V2;
		} while(s >= 1.0 || s == 0.0);
		
		X = V1 * sqrt(-2 * log(s) / s);
	} else
		X = V2 * sqrt(-2 * log(s) / s);
	
	phase = 1 - phase;
	
	return X;
}