lyapunovexponent.m

计算微分方程系统的Lyapunv exponent,在非线性动力学和有关混沌的研究中非常有用。

#include <iostream.h>
#include <fstream.h>
#include <math.h>
#define M 1000
#define N1 50
#define N2 50
#define PI 3.14159265

double *f(double t,double x[4])
{
double a=0.13,e=1.0,w=1.0,y=0.18,d=0.02;
double *p,m[4];
m[0]=x[1];
m[1]=x[0]-a*x[1]-e*x[0]*(x[0]*x[0]+x[2]*x[2])+y*cos(w*t);
m[2]=x[3];
m[3]=x[2]-a*x[3]-e*x[2]*(x[0]*x[0]+x[2]*x[2])+y*cos(w*t);
p=m;
return p;
}

double *Runge_Kutta(double t0,double t1,double x00,double x01,double x02,double x03,double h)
{
//Runge_Kutta法解微分方程组
//x'=f(x) length(x)=4
double t,x[4],*p;
x[0]=x00;
x[1]=x01;
x[2]=x02;
x[3]=x03;
double *k1,*k2,*k3,*k4;
for(t=t0;t<t1;t=t+h)
{
k1=f(t,x);
x[0]+=h*k1[0]/2;
x[1]+=h*k1[1]/2;
x[2]+=h*k1[2]/2;
x[3]+=h*k1[3]/2;
k2=f(t+h/2,x);
x[0]+=h*k2[0]/2;
x[1]+=h*k2[1]/2;
x[2]+=h*k2[2]/2;
x[3]+=h*k2[3]/2;
k3=f(t+h/2,x);
x[0]+=h*k3[0];
x[1]+=h*k3[1];
x[2]+=h*k3[2];
x[3]+=h*k3[3];
k4=f(t+h,x);
x[0]+=h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6;
x[1]+=h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6;
x[2]+=h*(k1[2]+2*k2[2]+2*k3[2]+k4[2])/6;
x[3]+=h*(k1[3]+2*k2[3]+2*k3[3]+k4[3])/6; 
}
p=x;
return p;
}
double Lyapunov_Exponent(double x[2],double z[2])
{
double delt=2*PI;
double d0=sqrt((x[0]-z[0])*(x[0]-z[0])+(x[1]-z[1])*(x[1]-z[1]));
double y[2],d[M];
for(int i=0;i<M;i++)
d=0;
for(i=0;i<M;i++)
{
x[0]=Runge_Kutta(i*delt,(i+1)*delt,x[0],x[1],x[0],x[1],0.001)[0];
x[1]=Runge_Kutta(i*delt,(i+1)*delt,x[0],x[1],x[0],x[1],0.001)[1];
y[0]=Runge_Kutta(i*delt,(i+1)*delt,z[0],z[1],z[0],z[1],0.001)[0];
y[1]=Runge_Kutta(i*delt,(i+1)*delt,z[0],z[1],z[0],z[1],0.001)[1];
d=sqrt((x[0]-y[0])*(x[0]-y[0])+(x[1]-y[1])*(x[1]-y[1]));
z[0]=x[0]+(d0/d)*(x[0]-y[0]);
z[1]=x[1]+(d0/d)*(x[1]-y[1]);
}
double le=0;
for(i=0;i<M;i++)
le+=log(d/d0);
le/=(M*delt);
return le;
}
void main()
{
float n=2,s=-2,w=-2,e=2;
float h1=(e-w)/N1,h2=(n-s)/N2;
double le[N1][N2];
double x[2],z[2];
for(int i=0;i<N1;i++)
for(int j=0;j<N2;j++)
le[j]=0;

for(i=0;i<N1;i++)
for(int j=0;j<N2;j++)
{
x[0]=i*h1+w;
x[1]=i*h2+s;
z[0]=(i+1)*h1+w;
z[1]=(i+1)*h2+s; 
le[j]=Lyapunov_Exponent(x,z);
cout<<"i="<<i<<"\tj="<<j<<"\n";
}
ofstream file;
file.open("Lyapunov_Exponent.txt");
if(!file)
{
cout<<"can't open file Lyapunov_Exponent.txt\n";
}
for(i=0;i<N1;i++)
{
for(int j=0;j<N2;j++)
{
cout<<le[j]<<"\t";
file<<le[j];
if(j==N2-1)
file<<";";
else
file<<",";
}
cout<<"\n";
}

}