sw.cpp

本程序模拟计算复杂网络中的混沌同步控制

/// hybrid  BBV  model programme
#include <stdio.h>
#include<stdlib.h>
#include <math.h>
#include <iostream>
using namespace std;

#define N 100                   //节点数
#define K 6				   //邻居数
double pp[18]={0.001 ,0.004 ,0.008 ,0.01, 0.04,0.08 ,0.1 ,0.4, 0.6 ,0.7 , 0.9 ,1}; 
//double pp[18]={0.5};
int Pw[10000][10000];
int Pw1[10000][10000];
long int i,j;
long int z,xx,t,tt;
long int bb;



//%%%%%%%%%%%%%%%%%%%%long ge tu ta   %%%%%%%%%%
double x[N][3];				   //各节点的状态
double e[N];				   //各节点的耦合输入
double d1;					   //耦合输入变量
double a1=20,a2=28,a3=8/3;
double c;         //耦合强度
double w=6.283185307179586476925286766559;
double h=0.001;
double X;


double f1(double x1,double x2,double x3)
//{ return(x2+c*d1);}
{return(a1*(x2-x1));}

double f2(double x1,double x2,double x3)
//{ return(-(a+b*cos(x3))*(a+b*cos(x3))*x1+kk/x1+1/pow(x1,3.0));}
{return(a2*x1-x2-x1*x3+c*X*d1);}

double f3(double x1,double x2,double x3)
//{ return(w);}
{return(x1*x2-a3*x3);}


void lungekuta()						  //积分
{
	int i;
	double k01,k02,k03,k11,k12,k13,k21,k22,k23,k31,k32,k33;
	for(i=0;i<N;i++)
	{
		d1=e[i];
		k01=f1(x[i][0],x[i][1],x[i][2]);
		k02=f2(x[i][0],x[i][1],x[i][2]);
		k03=f3(x[i][0],x[i][1],x[i][2]);
		k11=f1(x[i][0]+h*k01/2.0,x[i][1]+h*k02/2.0,x[i][2]+h*k03/2.0);
		k12=f2(x[i][0]+h*k01/2.0,x[i][1]+h*k02/2.0,x[i][2]+h*k03/2.0);
		k13=f3(x[i][0]+h*k01/2.0,x[i][1]+h*k02/2.0,x[i][2]+h*k03/2.0);
		k21=f1(x[i][0]+h*k11/2.0,x[i][1]+h*k12/2.0,x[i][2]+h*k13/2.0);
		k22=f2(x[i][0]+h*k11/2.0,x[i][1]+h*k12/2.0,x[i][2]+h*k13/2.0);
		k23=f3(x[i][0]+h*k11/2.0,x[i][1]+h*k12/2.0,x[i][2]+h*k13/2.0);
		k31=f1(x[i][0]+h*k21,x[i][1]+h*k22,x[i][2]+h*k23);
		k32=f2(x[i][0]+h*k21,x[i][1]+h*k22,x[i][2]+h*k23);
		k33=f3(x[i][0]+h*k21,x[i][1]+h*k22,x[i][2]+h*k23);
		x[i][0]=x[i][0]+h*(k01+2*k11+2*k21+k31)/6.0;
		x[i][1]=x[i][1]+h*(k02+2*k12+2*k22+k32)/6.0;
		x[i][2]=x[i][2]+h*(k03+2*k13+2*k23+k33)/6.0;
	}
}



//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


void main(void)
{


for (z=5;z<6;z++)
{


int t=0;
int tt=0;


   for(j=1;j<=K/2;j++)
	{
		for(i=0;i<N-j;i++)
		{
			Pw[i][i+j]=1;    
			Pw[i+j][i]=1;    
		}
		for(i=N-j;i<N;i++)
		{
			Pw[i][(i+j)%N]=1;    
			Pw[(i+j)%N][i]=1;    
		}
	}


  for(i=0;i<N;i++)
  {
	for (j=0;j<N;j++)
	{
		Pw1[i][j]=Pw[i][j];
	}
  }


   double p;
   p=pp[z];
    
   for (i=0;i<N;i++)
   {
      for (j=i+1;j<N;j++)
	  {
	    if (Pw1[i][j]==1)
		{
		   double r0;
           r0=double(rand())/32767;
		  // cout<<"r0=  "<<r0<<endl;

         if (r0<=0.9)
		 {		
			 tt=tt+1;
			 Pw[i][j]=0;
             Pw[j][i]=0;

          for (xx=0;xx<100000;xx++)
		  {
			   bb=long int(rand()%(N-1))+1;
				
			   if (bb!=i)
			   {
					 if (bb!=j)
					 {
						 if (Pw1[i][bb]==0) 	
						 {
							if (Pw[i][bb]==0)	
							{
						
						//	cout<<"bb=  "<<bb<<endl;
								break;
							}
							
						 }
					 }
					 
			   }
		  }  
	  Pw[i][bb]=1;
	  Pw[bb][i]=1;
		    
		 }
		}
	  }
   }




 /*  for(i=0;i<N;i++)
   {
	for (j=0;j<N;j++)
	{
		if (Pw[i][j]==1)
		{
			t=t+1;
		}
			cout<<Pw[i][j];

	}
	cout<<endl;

   }

cout<<endl;
	cout<<"t=   "<<t;
cout<<endl;
	cout<<"tt=  "<<tt;  */

///  龙格库塔法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


	FILE   *fp;
	fp=fopen("dmax.dat","w");
    FILE   *fp1;
	fp1=fopen("xmax.dat","w");

	int    b,kkk,a[N][N],k;
	double ddmax;
    double dmax[30000];
	double xmax[30000][3];
    
    for(i=0;i<N;i++)
	{
		for (j=0;j<N;j++)
		{
             a[i][j]=Pw[i][j];
		}
	}
	

	for(i=0;i<N;i++)								  //对角元素
	{   
		b=0;
		for(j=0;j<i;j++) 
		{
			if (a[i][j]!=0)  b=b+1;
		}
        for(j=i+1;j<N;j++) 
		{
		    if (a[i][j]!=0)  b=b+1;
		}
		a[i][i]=-b;
	} 



												
	for(c=39;c<39.01;c+=0.1)             	 //耦合强度：1.0--2.0,每次增加0.1
	{
		for(i=0;i<N;i++)							 //各节点的初始状态
		{ 
			x[i][0]=double(rand())/32768;
			x[i][1]=double(rand())/32768;
			x[i][2]=double(rand())/32768;
			
		
		}
			//for (j=0;j<N;j++)
			 //cout<<x[j][1]<<endl;



		for(kkk=0;kkk<30000;kkk++)					   //时步
		{
			cout<<kkk<<endl;
			for(i=0;i<N;i++)
			{
				e[i]=0.0;
				for(j=0;j<N;j++) 
				e[i]+=a[i][j]*x[j][1];		   //第i个节点的耦合输入
			}



   	for (i=0;i<1;i++)
	{
	
   	 double V1=0, V2=0, S=0;
     int phase = 0;
     double X;
     if ( phase == 0 ) 
	{
        do 
		{
		 //   srand((unsigned)time(0));
            double U1 = (double)rand() /32767;
            double U2 = (double)rand() / 32767;
            
            V1 = 2 * U1 - 1;
            V2 = 2 * U2 - 1;
            S = V1 * V1 + V2 * V2;
        } 
		while(S >= 1 || S == 0);
        X = V1 * sqrt(-2 * log(S) / S);
    } else
        X = V2 * sqrt(-2 * log(S) / S);
        
    phase = 1 - phase;
    cout<<X<<endl;
	}


			lungekuta();
			xmax[kkk][0]=0.0;
            xmax[kkk][1]=0.0;
			xmax[kkk][2]=0.0;
			xmax[kkk][0]=x[2][0]-x[N-N/2][0];
            xmax[kkk][1]=x[2][1]-x[N-N/2][1];
            xmax[kkk][2]=x[2][2]-x[N-N/2][2];

		for(i=0;i<N;i++)
		{
			for(j=i+1;j<N;j++)
			{
				ddmax=0;
				dmax[kkk]=0;
				for(k=0;k<3;k++) ddmax+=(x[i][k]-x[j][k])*(x[i][k]-x[j][k]);	  
				ddmax=sqrt(ddmax);
				if(ddmax>dmax[kkk]) dmax[kkk]=ddmax;					 //所有节点对的轨道距离中最大者
			} 
		}
	//	fprintf(fp,"%f\n",c); 
		fprintf(fp,"%f\n",dmax[kkk]); 
		fprintf(fp1,"  %f  ",xmax[kkk][0]); 
		fprintf(fp1,"  %f  ",xmax[kkk][1]); 
		fprintf(fp1,"  %f \n ",xmax[kkk][2]); 
       // cout<<dmax[k]<<endl;
	

		}
    
	}  ///  for c

	fclose(fp);
    fclose(fp1);
}  // for z

for (j=0;j<N;j++)
{
	cout<<x[j][0];
    cout<<"  "<<x[j][1];
    cout<<"  " <<  x[j][2]<<endl;
}  

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

}  // void