lce.c

经典混沌系统Lorenz系统的Lyapunov指数计算C程序

//Lorenz系统的Lyapunov指数谱
//success
#include "stdio.h"
#include "math.h"

long   N=40000;
double h=0.02;

void dverk(y)
double y[12];
{
	int    i,j;
	double y1[12],y2[12],k[8][12],b[8][12];
	double c[8]={41.0,0.0,216.0,27.0,272.0,27.0,216.0,41.0};
	for(j=0;j<12;j++)  y1[j]=y[j];
	for(i=0;i<8;i++)
	{
		for(j=0;j<12;j++)
		{
			b[0][j]=0.0;
			b[1][j]=k[0][j]/9.0;
			b[2][j]=(k[0][j]+3.0*k[1][j])/24.0;
			b[3][j]=(k[0][j]-3.0*k[1][j]+4.0*k[2][j])/6.0;
			b[4][j]=(-5.0*k[0][j]+27.0*k[1][j]-24.0*k[2][j]+6.0*k[3][j])/8.0;
			b[5][j]=(221.0*k[0][j]-981.0*k[1][j]+867.0*k[2][j]-102.0*k[3][j]+k[4][j])/9.0;
			b[6][j]=(-183.0*k[0][j]+678.0*k[1][j]-472.0*k[2][j]-66.0*k[3][j]+80.0*k[4][j]+3.0*k[5][j])/48.0;
			b[7][j]=(716.0*k[0][j]-2079.0*k[1][j]+1002.0*k[2][j]+834.0*k[3][j]-454.0*k[4][j]-9.0*k[5][j]+72.0*k[6][j])/82.0;
			y2[j]=y1[j]+b[i][j]*h;
		}
		k[i][0]=16.0*(y2[1]-y2[0]);
		k[i][1]=-y2[0]*y2[2]+45.92*y2[0]-y2[1];
		k[i][2]=y2[0]*y2[1]-4.0*y2[2];
		for(j=0;j<3;j++)
		{
			k[i][3+j]=16.0*(y2[6+j]-y2[3+j]);
			k[i][6+j]=(45.92-y2[2])*y2[3+j]-y2[6+j]-y2[0]*y2[9+j];
			k[i][9+j]=y2[1]*y2[3+j]+y2[0]*y2[6+j]-4.0*y2[9+j];
		}
		for(j=0;j<12;j++)  y[j]+=h*c[i]*k[i][j]/840.0;
	}
}

void gs(y,znorm)
double y[12];
double znorm[3];
{
	int    i,j,k;
	double gsc[2];
	znorm[0]=0.0;
	for(i=1;i<=3;i++) znorm[0]+=y[3*i]*y[3*i];
	znorm[0]=sqrt(znorm[0]);
	for(i=1;i<=3;i++) y[3*i]/=znorm[0];
	for(i=1;i<3;i++)
	{
		for(k=0;k<=i-1;k++)
		{
			gsc[k]=0.0;
			for(j=1;j<=3;j++) gsc[k]+=y[3*j+i]*y[3*j+k];
		}
		for(k=1;k<=3;k++)
		{
			for(j=0;j<=i-1;j++) y[3*k+i]-=gsc[j]*y[3*k+j];
		}
		znorm[i]=0.0;
		for(k=1;k<=3;k++) znorm[i]+=y[3*k+i]*y[3*k+i];
		znorm[i]=sqrt(znorm[i]);
		for(k=1;k<=3;k++) y[3*k+i]/=znorm[i];
	}
}

void main()
{
	long   i;
	int    j;
	double x[12]={9.23667,14.04743,29.92179,1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};
	double z[3],cum[3]={0.0,0.0,0.0};
	for(i=0;i<N;i++)
	{
		dverk(x);
		gs(x,z);
		for(j=0;j<3;j++) cum[j]+=log(z[j]);
	}
	for(j=0;j<3;j++)
	{
		cum[j]/=(N*h*log(2.0));
		printf("%f\n",cum[j]);
	}
}