henon.cpp

将henon映射由牛顿算法进行控制,求出其周期2的点,这个例子很有代表性

#include <stdio.h>
#include <stdlib.h>
#include "math.h"

#define A 1.4
#define B 0.3
#define pnumber 1000
#ifndef EPSIOLON
#define EPSILON 0.000000001
#endif

#define DIM 2

double map(double xn[])
{
	double x,y;
	x=A-xn[0]*xn[0]+B*xn[1];
	y=xn[0];

	xn[0]=x;
	xn[1]=y;

	return *xn;
}

double onemap(double xn[])
{
	double x,y;
	x=A-xn[1]*xn[1]+B*xn[0];
	y=A-xn[0]*xn[0]+B*xn[1];

	xn[0]=x;
	xn[1]=y;

	return *xn;
}
double df(double xn[],double df[],int i)
{
	if(i=2)
	{
	df[0]=1-B;
	df[1]=2*xn[1];
	df[2]=2*xn[0];
	df[3]=1-B;
	}
	return *df;
}
int gjcpeli( int process, double a[DIM][DIM], double xx[DIM] )
{
	int k,i,j,i0;
double pelement;
if( process == 1 )
    printf("The process of elimination\n");
for(k=0;k<DIM;k++)
  {
  pelement=fabs(a[k][k]);  i0=k;
  for(i=k+1;i<DIM;i++)
    {
    if( fabs(a[i][k]) > pelement )
      { pelement=fabs(a[i][k]); i0=i; }
  }
  if( i0 != k )
    {
    for(j=k;j<DIM;j++)
      { pelement=a[k][j]; a[k][j]=a[i0][j]; a[i0][j]=pelement; }
    pelement=xx[k]; xx[k]=xx[i0]; xx[i0]=pelement;
    }
  if( fabs(a[k][k]) < EPSILON ) return(1);
  for(j=k+1;j<DIM; j++ ) a[k][j]=a[k][j]/a[k][k];
  xx[k]=xx[k]/a[k][k];
  a[k][k]=1.0;
  for(i=0;i<DIM;i++)
    {
    if( i!=k )
     {
     xx[i]=xx[i]-a[i][k]*xx[k];
     for(j=k+1;j<DIM;j++) a[i][j]=a[i][j]-a[i][k]*a[k][j];
     a[i][k]=0.0;
     }
    }
  if( process == 1 )
    {
    for(i=0;i<DIM;i++)
      {
      for(j=0;j<DIM;j++) printf("%10.4f",a[i][j]);
      printf("   |%10.4f\n",xx[i]);
      }
    printf("\n");
    }
  }
	return 0;
}

double delt(double dfx[],double xn[],double kf[],double delta[])
{
int i,j;
double a[DIM][DIM];	
double b[DIM];
for(i=0;i<DIM;i++)
{
	for(j=0;j<DIM;j++)
	{
		a[i][j]=dfx[i*DIM+j];
	}

}
for(i=0;i<DIM;i++)
{
	 b[i]=xn[i]-kf[i];
}

if(gjcpeli( 1, a, b ) == 1)
  {
  printf(" The linear system has'nt solution!\n");
  printf(" Strike any key to exit!\n");  exit(1);
  }
printf("Gauss-Jordan column pivot\n");
for(i=0;i<DIM;i++)  printf("    %10.6f\n",b[i]);
for(i=0;i<DIM;i++)
{
	 delta[i]=b[i];
}
 return 0;
}

void main()
{
	FILE   *fp1;
	if((fp1=fopen("point.dat","wr"))==NULL)
	{
	 printf("Can't open reading file 1 !");
	 exit(0);
	}

	double x,y;
	x=1.4; y=-0.6;
//	x=random()*1.5;
	double xn[DIM];
	double dfx[4];
	double kf[DIM];
	double delta[DIM];

	for(int i=0; i<4; i++)
	{
		dfx[i]=0;
	}
	for( i=0; i<DIM; i++)
	{
		delta[i]=0;
	}

	xn[0]=x;
	xn[1]=y;

	for(i=0; i<pnumber; i++)
	{
		kf[0]=xn[0];
		kf[1]=xn[1];

		df(xn,dfx,2);

		onemap(xn);

		//df(xn,dfx,2);

		delt(dfx, xn, kf, delta);

		if(fabs(delta[0])<0.000001 && fabs(delta[1])<0.000001)
		{
			i=pnumber;
			printf("    %10.6f, %10.6f\n",xn[0],xn[1]);

		}
		else
		{
		xn[0]=xn[0]+delta[0];
		xn[1]=xn[1]+delta[1];
		if (i==pnumber)
		{
			printf("lose");
		}
		}
	//	fprintf(fp1,"%lf %lf\n",xn[0],xn[1]);
		//	     printf("okkkkkkkkkkkkkkkkkk");

	}
//	double z=xn[0],u=xn[1];
	

	fclose(fp1);
}